probabilidad - Ajustando la distribución empírica a las teóricas con Scipy(Python)?

poisson en python

Question

INTRODUCCIÓN: Tengo una lista de más de 30 000 valores que van de 0 a 47, por ejemplo [0,0,0,0, .., 1,1,1,1, ..., 2,2,2,2, ..., 47 etc.] que es la distribución continua.

PROBLEMA: De acuerdo con 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 a 1 y el valor p para números más altos tendrían tendencia a 0.

No sé si estoy en lo correcto, 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 forma de implementar dicho análisis en Python (Scipy o Numpy)? ¿Podría presentar algún ejemplo?

¡Gracias!

Answer 1

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 distribuciones scipy.stats 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 ajuste

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

Todas las distribuciones

La mejor distribución de 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 Propbability 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_fir_paramms = 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
pdf = make_pdf(best_dist, best_fir_paramms)

# 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_fir_paramms)])
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'')

Answer 2

AFAICU, tu distribución es discreta (y nada discreta). Por lo tanto, basta con contar las frecuencias de diferentes valores y normalizarlos 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 superiores a 1 es simplemente (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.

Answer 3

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 el número de ocurrencias 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.

Answer 4

Suena como un problema de estimación de densidad de probabilidad para mí.

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:])

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

Answer 5

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

1, el que te da la mayor probabilidad de registro.

2, el que te da los valores más pequeños de AIC, BIC o BICc (ver wiki: http://en.wikipedia.org/wiki/Akaike_information_criterion , básicamente se puede ver como log likelihood ajustado por el número de parámetros, como la 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 debe describir sus datos (según las teorías en su campo particular) y desea apegarse a eso, omitirá el paso de identificar la mejor distribución de ajuste.

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

Answer 6

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 método fit() . Verifique el código a continuación para obtener más detalles:

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), color=''w'')

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 ajustadas, 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'']