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!

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

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')

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:

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

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?

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.

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 .

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()

# 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()

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.

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:

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.

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.