¿Adaptar la distribución empírica a las teóricas con Scipy (Python)?

139

INTRODUCCIÓN : Tengo una lista de más de 30,000 valores enteros que van de 0 a 47, inclusive, por ejemplo, [0,0,0,0,..,1,1,1,1,...,2,2,2,2,...,47,47,47,...]muestreados de alguna distribución continua. Los valores en la lista no están necesariamente en orden, pero el orden no importa para este problema.

PROBLEMA : Según mi distribución, me gustaría calcular el valor p (la probabilidad de ver valores mayores) para cualquier valor dado. Por ejemplo, como puede ver, el valor p para 0 se acercará a 1 y el valor p para números más altos tenderá a 0.

No sé si tengo razón, pero para determinar las probabilidades creo que necesito ajustar mis datos a una distribución teórica que sea la más adecuada para describir mis datos. Supongo que se necesita algún tipo de prueba de bondad de ajuste para determinar el mejor modelo.

¿Hay alguna manera de implementar dicho análisis en Python ( Scipyo Numpy)? ¿Podría presentar algún ejemplo?

¡Gracias!

s_sherly
fuente
2
¿Solo tiene valores empíricos discretos pero desea una distribución continua? ¿Entiendo eso correctamente?
Michael J. Barber
1
Parece absurdo. ¿Qué representan los números? Mediciones con precisión limitada?
Michael J. Barber
1
Michael, le expliqué lo que representan los números en mi pregunta anterior: stackoverflow.com/questions/6615489/…
s_sherly
66
Esos son los datos de conteo. No es una distribución continua.
Michael J. Barber
1
Compruebe la respuesta aceptada a esta pregunta stackoverflow.com/questions/48455018/…
Ahmad Suliman

Respuestas:

209

Ajuste de distribución con suma de error cuadrado (SSE)

Esta es una actualización y modificación de la respuesta de Saullo , que utiliza la lista completa de las scipy.statsdistribuciones actuales y devuelve la distribución con la menor SSE entre el histograma de la distribución y el histograma de los datos.

Ejemplo de montaje

Usando el conjunto de datos de El Niño destatsmodels , las distribuciones se ajustan y se determina el error. Se devuelve la distribución con el menor error.

Todas las distribuciones

Todas las distribuciones ajustadas

Distribución de mejor ajuste

Distribución de mejor ajuste

Código de ejemplo

%matplotlib inline

import warnings
import numpy as np
import pandas as pd
import scipy.stats as st
import statsmodels as sm
import matplotlib
import matplotlib.pyplot as plt

matplotlib.rcParams['figure.figsize'] = (16.0, 12.0)
matplotlib.style.use('ggplot')

# Create models from data
def best_fit_distribution(data, bins=200, ax=None):
    """Model data by finding best fit distribution to data"""
    # Get histogram of original data
    y, x = np.histogram(data, bins=bins, density=True)
    x = (x + np.roll(x, -1))[:-1] / 2.0

    # Distributions to check
    DISTRIBUTIONS = [        
        st.alpha,st.anglit,st.arcsine,st.beta,st.betaprime,st.bradford,st.burr,st.cauchy,st.chi,st.chi2,st.cosine,
        st.dgamma,st.dweibull,st.erlang,st.expon,st.exponnorm,st.exponweib,st.exponpow,st.f,st.fatiguelife,st.fisk,
        st.foldcauchy,st.foldnorm,st.frechet_r,st.frechet_l,st.genlogistic,st.genpareto,st.gennorm,st.genexpon,
        st.genextreme,st.gausshyper,st.gamma,st.gengamma,st.genhalflogistic,st.gilbrat,st.gompertz,st.gumbel_r,
        st.gumbel_l,st.halfcauchy,st.halflogistic,st.halfnorm,st.halfgennorm,st.hypsecant,st.invgamma,st.invgauss,
        st.invweibull,st.johnsonsb,st.johnsonsu,st.ksone,st.kstwobign,st.laplace,st.levy,st.levy_l,st.levy_stable,
        st.logistic,st.loggamma,st.loglaplace,st.lognorm,st.lomax,st.maxwell,st.mielke,st.nakagami,st.ncx2,st.ncf,
        st.nct,st.norm,st.pareto,st.pearson3,st.powerlaw,st.powerlognorm,st.powernorm,st.rdist,st.reciprocal,
        st.rayleigh,st.rice,st.recipinvgauss,st.semicircular,st.t,st.triang,st.truncexpon,st.truncnorm,st.tukeylambda,
        st.uniform,st.vonmises,st.vonmises_line,st.wald,st.weibull_min,st.weibull_max,st.wrapcauchy
    ]

    # Best holders
    best_distribution = st.norm
    best_params = (0.0, 1.0)
    best_sse = np.inf

    # Estimate distribution parameters from data
    for distribution in DISTRIBUTIONS:

        # Try to fit the distribution
        try:
            # Ignore warnings from data that can't be fit
            with warnings.catch_warnings():
                warnings.filterwarnings('ignore')

                # fit dist to data
                params = distribution.fit(data)

                # Separate parts of parameters
                arg = params[:-2]
                loc = params[-2]
                scale = params[-1]

                # Calculate fitted PDF and error with fit in distribution
                pdf = distribution.pdf(x, loc=loc, scale=scale, *arg)
                sse = np.sum(np.power(y - pdf, 2.0))

                # if axis pass in add to plot
                try:
                    if ax:
                        pd.Series(pdf, x).plot(ax=ax)
                    end
                except Exception:
                    pass

                # identify if this distribution is better
                if best_sse > sse > 0:
                    best_distribution = distribution
                    best_params = params
                    best_sse = sse

        except Exception:
            pass

    return (best_distribution.name, best_params)

def make_pdf(dist, params, size=10000):
    """Generate distributions's Probability Distribution Function """

    # Separate parts of parameters
    arg = params[:-2]
    loc = params[-2]
    scale = params[-1]

    # Get sane start and end points of distribution
    start = dist.ppf(0.01, *arg, loc=loc, scale=scale) if arg else dist.ppf(0.01, loc=loc, scale=scale)
    end = dist.ppf(0.99, *arg, loc=loc, scale=scale) if arg else dist.ppf(0.99, loc=loc, scale=scale)

    # Build PDF and turn into pandas Series
    x = np.linspace(start, end, size)
    y = dist.pdf(x, loc=loc, scale=scale, *arg)
    pdf = pd.Series(y, x)

    return pdf

# Load data from statsmodels datasets
data = pd.Series(sm.datasets.elnino.load_pandas().data.set_index('YEAR').values.ravel())

# Plot for comparison
plt.figure(figsize=(12,8))
ax = data.plot(kind='hist', bins=50, normed=True, alpha=0.5, color=plt.rcParams['axes.color_cycle'][1])
# Save plot limits
dataYLim = ax.get_ylim()

# Find best fit distribution
best_fit_name, best_fit_params = best_fit_distribution(data, 200, ax)
best_dist = getattr(st, best_fit_name)

# Update plots
ax.set_ylim(dataYLim)
ax.set_title(u'El Niño sea temp.\n All Fitted Distributions')
ax.set_xlabel(u'Temp (°C)')
ax.set_ylabel('Frequency')

# Make PDF with best params 
pdf = make_pdf(best_dist, best_fit_params)

# Display
plt.figure(figsize=(12,8))
ax = pdf.plot(lw=2, label='PDF', legend=True)
data.plot(kind='hist', bins=50, normed=True, alpha=0.5, label='Data', legend=True, ax=ax)

param_names = (best_dist.shapes + ', loc, scale').split(', ') if best_dist.shapes else ['loc', 'scale']
param_str = ', '.join(['{}={:0.2f}'.format(k,v) for k,v in zip(param_names, best_fit_params)])
dist_str = '{}({})'.format(best_fit_name, param_str)

ax.set_title(u'El Niño sea temp. with best fit distribution \n' + dist_str)
ax.set_xlabel(u'Temp. (°C)')
ax.set_ylabel('Frequency')
tmthydvnprt
fuente
2
Increíble. Considere usar en density=Truelugar de normed=Trueen np.histogram(). ^^
Peque
1
@tmthydvnprt Tal vez podría deshacer los cambios en los .plot()métodos para evitar futuras confusiones. ^^
Peque
10
Para obtener los nombres de distribución: from scipy.stats._continuous_distns import _distn_names. Luego puede usar algo como getattr(scipy.stats, distname)para cada uno distnameen _distn_names`. Útil porque las distribuciones se actualizan con diferentes versiones de SciPy.
Brad Solomon el
1
¿Puede explicar por qué este código verifica solo el mejor ajuste de las distribuciones continuas y no puede verificar las distribuciones discretas o multivariadas? Gracias.
Adam Schroeder
66
Muy genial. Tuve que actualizar el parámetro de color -ax = data.plot(kind='hist', bins=50, normed=True, alpha=0.5, color=list(matplotlib.rcParams['axes.prop_cycle'])[1]['color'])
basswaves
147

Hay 82 funciones de distribución implementadas en SciPy 0.12.0 . Puede probar cómo algunos de ellos se ajustan a sus datos utilizando su fit()método . Consulte el código a continuación para obtener más detalles:

ingrese la descripción de la imagen aquí

import matplotlib.pyplot as plt
import scipy
import scipy.stats
size = 30000
x = scipy.arange(size)
y = scipy.int_(scipy.round_(scipy.stats.vonmises.rvs(5,size=size)*47))
h = plt.hist(y, bins=range(48))

dist_names = ['gamma', 'beta', 'rayleigh', 'norm', 'pareto']

for dist_name in dist_names:
    dist = getattr(scipy.stats, dist_name)
    param = dist.fit(y)
    pdf_fitted = dist.pdf(x, *param[:-2], loc=param[-2], scale=param[-1]) * size
    plt.plot(pdf_fitted, label=dist_name)
    plt.xlim(0,47)
plt.legend(loc='upper right')
plt.show()

Referencias

- Distribuciones de ajuste, bondad de ajuste, valor p. ¿Es posible hacer esto con Scipy (Python)?

- Ajuste de distribución con Scipy

Y aquí una lista con los nombres de todas las funciones de distribución disponibles en Scipy 0.12.0 (VI):

dist_names = [ 'alpha', 'anglit', 'arcsine', 'beta', 'betaprime', 'bradford', 'burr', 'cauchy', 'chi', 'chi2', 'cosine', 'dgamma', 'dweibull', 'erlang', 'expon', 'exponweib', 'exponpow', 'f', 'fatiguelife', 'fisk', 'foldcauchy', 'foldnorm', 'frechet_r', 'frechet_l', 'genlogistic', 'genpareto', 'genexpon', 'genextreme', 'gausshyper', 'gamma', 'gengamma', 'genhalflogistic', 'gilbrat', 'gompertz', 'gumbel_r', 'gumbel_l', 'halfcauchy', 'halflogistic', 'halfnorm', 'hypsecant', 'invgamma', 'invgauss', 'invweibull', 'johnsonsb', 'johnsonsu', 'ksone', 'kstwobign', 'laplace', 'logistic', 'loggamma', 'loglaplace', 'lognorm', 'lomax', 'maxwell', 'mielke', 'nakagami', 'ncx2', 'ncf', 'nct', 'norm', 'pareto', 'pearson3', 'powerlaw', 'powerlognorm', 'powernorm', 'rdist', 'reciprocal', 'rayleigh', 'rice', 'recipinvgauss', 'semicircular', 't', 'triang', 'truncexpon', 'truncnorm', 'tukeylambda', 'uniform', 'vonmises', 'wald', 'weibull_min', 'weibull_max', 'wrapcauchy'] 
Saullo GP Castro
fuente
77
¿Qué pasa si normed = Trueal trazar el histograma? No te multiplicarías pdf_fittedpor el size, ¿verdad?
aloha
3
Vea esta respuesta si desea ver cómo se ven todas las distribuciones o para tener una idea de cómo acceder a todas ellas.
tmthydvnprt
@SaulloCastro ¿Qué representan los 3 valores en param, en la salida de dist.fit?
shaifali Gupta
2
Para obtener los nombres de distribución: from scipy.stats._continuous_distns import _distn_names. Luego puede usar algo como getattr(scipy.stats, distname)para cada uno distnameen _distn_names`. Útil porque las distribuciones se actualizan con diferentes versiones de SciPy.
Brad Solomon el
1
Quitaría color = 'w' del código; de lo contrario, el histograma no se muestra.
Eran
12

fit()El método mencionado por @Saullo Castro proporciona estimaciones de máxima verosimilitud (MLE). La mejor distribución para sus datos es la que le brinda la más alta puede determinarse de varias maneras diferentes: como

1, el que le brinda la mayor probabilidad de registro.

2, el que le proporciona los valores más pequeños de AIC, BIC o BICc (consulte wiki: http://en.wikipedia.org/wiki/Akaike_information_criterion , básicamente se puede ver como la probabilidad de registro ajustada por el número de parámetros, como distribución con más se espera que los parámetros se ajusten mejor)

3, el que maximiza la probabilidad posterior bayesiana. (ver wiki: http://en.wikipedia.org/wiki/Posterior_probability )

Por supuesto, si ya tiene una distribución que debería describir sus datos (en base a las teorías en su campo particular) y desea apegarse a eso, omitirá el paso de identificar la distribución de mejor ajuste.

scipyno viene con una función para calcular la probabilidad de registro (aunque se proporciona el método MLE), pero el código uno es fácil: consulte ¿Las funciones de densidad de probabilidad incorporadas de `scipy.stat.distributions` son más lentas que las proporcionadas por el usuario?

CT Zhu
fuente
1
¿Cómo aplicaría este método a una situación en la que los datos ya se han agrupado, es decir, ya es un histograma en lugar de generar un histograma a partir de los datos?
Pete
@pete, esa sería una situación de datos censurados por intervalos, hay un método de máxima probabilidad para ello, pero actualmente no está implementado enscipy
CT Zhu
No olvides la evidencia
jtlz2
5

AFAICU, su distribución es discreta (y nada más que discreta). Por lo tanto, solo contar las frecuencias de diferentes valores y normalizarlos debería ser suficiente para sus propósitos. Entonces, un ejemplo para demostrar esto:

In []: values= [0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 3, 3, 4]
In []: counts= asarray(bincount(values), dtype= float)
In []: cdf= counts.cumsum()/ counts.sum()

Por lo tanto, la probabilidad de ver valores más altos que 1simplemente (de acuerdo con la función de distribución acumulativa complementaria (ccdf) :

In []: 1- cdf[1]
Out[]: 0.40000000000000002

Tenga en cuenta que ccdf está estrechamente relacionado con la función de supervivencia (sf) , pero también se define con distribuciones discretas, mientras que sf se define solo para distribuciones contiguas.

comer
fuente
2

A mí me parece un problema de estimación de densidad de probabilidad.

from scipy.stats import gaussian_kde
occurences = [0,0,0,0,..,1,1,1,1,...,2,2,2,2,...,47]
values = range(0,48)
kde = gaussian_kde(map(float, occurences))
p = kde(values)
p = p/sum(p)
print "P(x>=1) = %f" % sum(p[1:])

Consulte también http://jpktd.blogspot.com/2009/03/using-gaussian-kernel-density.html .

emre
fuente
1
Para futuros lectores: esta solución (o al menos la idea) proporciona la respuesta más simple a las preguntas de los OP ('cuál es el valor p'); sería interesante saber cómo se compara esto con algunos de los métodos más involucrados que se ajustan Una distribución conocida.
Greg
¿Las regresiones de kernel gaussianas funcionan para todas las distribuciones?
@mikey Como regla general, no hay regresiones que funcionen para todas las distribuciones. Sin embargo
TheEnvironmentalist
2

Prueba la distfitbiblioteca.

pip install distfit

# Create 1000 random integers, value between [0-50]
X = np.random.randint(0, 50,1000)

# Retrieve P-value for y
y = [0,10,45,55,100]

# From the distfit library import the class distfit
from distfit import distfit

# Initialize.
# Set any properties here, such as alpha.
# The smoothing can be of use when working with integers. Otherwise your histogram
# may be jumping up-and-down, and getting the correct fit may be harder.
dist = distfit(alpha=0.05, smooth=10)

# Search for best theoretical fit on your empirical data
dist.fit_transform(X)

> [distfit] >fit..
> [distfit] >transform..
> [distfit] >[norm      ] [RSS: 0.0037894] [loc=23.535 scale=14.450] 
> [distfit] >[expon     ] [RSS: 0.0055534] [loc=0.000 scale=23.535] 
> [distfit] >[pareto    ] [RSS: 0.0056828] [loc=-384473077.778 scale=384473077.778] 
> [distfit] >[dweibull  ] [RSS: 0.0038202] [loc=24.535 scale=13.936] 
> [distfit] >[t         ] [RSS: 0.0037896] [loc=23.535 scale=14.450] 
> [distfit] >[genextreme] [RSS: 0.0036185] [loc=18.890 scale=14.506] 
> [distfit] >[gamma     ] [RSS: 0.0037600] [loc=-175.505 scale=1.044] 
> [distfit] >[lognorm   ] [RSS: 0.0642364] [loc=-0.000 scale=1.802] 
> [distfit] >[beta      ] [RSS: 0.0021885] [loc=-3.981 scale=52.981] 
> [distfit] >[uniform   ] [RSS: 0.0012349] [loc=0.000 scale=49.000] 

# Best fitted model
best_distr = dist.model
print(best_distr)

# Uniform shows best fit, with 95% CII (confidence intervals), and all other parameters
> {'distr': <scipy.stats._continuous_distns.uniform_gen at 0x16de3a53160>,
>  'params': (0.0, 49.0),
>  'name': 'uniform',
>  'RSS': 0.0012349021241149533,
>  'loc': 0.0,
>  'scale': 49.0,
>  'arg': (),
>  'CII_min_alpha': 2.45,
>  'CII_max_alpha': 46.55}

# Ranking distributions
dist.summary

# Plot the summary of fitted distributions
dist.plot_summary()

ingrese la descripción de la imagen aquí

# Make prediction on new datapoints based on the fit
dist.predict(y)

# Retrieve your pvalues with 
dist.y_pred
# array(['down', 'none', 'none', 'up', 'up'], dtype='<U4')
dist.y_proba
array([0.02040816, 0.02040816, 0.02040816, 0.        , 0.        ])

# Or in one dataframe
dist.df

# The plot function will now also include the predictions of y
dist.plot()

Mejor ajuste

Tenga en cuenta que en este caso, todos los puntos serán significativos debido a la distribución uniforme. Puede filtrar con dist.y_pred si es necesario.

erdogant
fuente
1

Con OpenTURNS , usaría los criterios BIC para seleccionar la mejor distribución que se ajuste a dichos datos. Esto se debe a que este criterio no da demasiada ventaja a las distribuciones que tienen más parámetros. De hecho, si una distribución tiene más parámetros, es más fácil para la distribución ajustada estar más cerca de los datos. Además, el Kolmogorov-Smirnov puede no tener sentido en este caso, porque un pequeño error en los valores medidos tendrá un gran impacto en el valor p.

Para ilustrar el proceso, cargo los datos de El-Nino, que contienen 732 mediciones de temperatura mensuales desde 1950 hasta 2010:

import statsmodels.api as sm
dta = sm.datasets.elnino.load_pandas().data
dta['YEAR'] = dta.YEAR.astype(int).astype(str)
dta = dta.set_index('YEAR').T.unstack()
data = dta.values

Es fácil obtener las 30 fábricas de distribuciones univariadas incorporadas con el GetContinuousUniVariateFactoriesmétodo estático. Una vez hecho, el BestModelBICmétodo estático devuelve el mejor modelo y la puntuación BIC correspondiente.

sample = ot.Sample(data, 1)
tested_factories = ot.DistributionFactory.GetContinuousUniVariateFactories()
best_model, best_bic = ot.FittingTest.BestModelBIC(sample,
                                                   tested_factories)
print("Best=",best_model)

que imprime:

Best= Beta(alpha = 1.64258, beta = 2.4348, a = 18.936, b = 29.254)

Para comparar gráficamente el ajuste al histograma, utilizo los drawPDFmétodos de la mejor distribución.

import openturns.viewer as otv
graph = ot.HistogramFactory().build(sample).drawPDF()
bestPDF = best_model.drawPDF()
bestPDF.setColors(["blue"])
graph.add(bestPDF)
graph.setTitle("Best BIC fit")
name = best_model.getImplementation().getClassName()
graph.setLegends(["Histogram",name])
graph.setXTitle("Temperature (°C)")
otv.View(graph)

Esto produce:

Beta ajustado a las temperaturas de El-Nino

Más detalles sobre este tema se presentan en el documento BestModelBIC . Sería posible incluir la distribución Scipy en SciPyDistribution o incluso con distribuciones ChaosPy con ChaosPyDistribution , pero supongo que el script actual cumple con la mayoría de los propósitos prácticos.

Michael Baudin
fuente
2
Probablemente deberías declarar un interés?
jtlz2
0

Perdóneme si no entiendo su necesidad, pero ¿qué pasa con el almacenamiento de sus datos en un diccionario donde las claves serían los números entre 0 y 47 y valora el número de apariciones de sus claves relacionadas en su lista original?
Por lo tanto, su probabilidad p (x) será la suma de todos los valores para claves mayores que x dividido por 30000.

PierrOz
fuente
En este caso, la p (x) será la misma (igual a 0) para cualquier valor mayor que 47. Necesito una distribución de probabilidad continua.
s_sherly
2
@s_sherly: probablemente sería algo bueno si pudiera editar y aclarar mejor su pregunta, ya que, de hecho, la "probabilidad de ver valores más grandes" , como usted dice, ES cero para los valores que están por encima del valor más alto en el grupo .
mac