INTRODUCTION : J'ai une liste de plus de 30 000 valeurs entières allant de 0 à 47, inclusivement, par exemple [0,0,0,0,..,1,1,1,1,...,2,2,2,2,...,47,47,47,...]échantillonnées à partir d'une distribution continue. Les valeurs de la liste ne sont pas nécessairement dans l'ordre, mais l'ordre n'a pas d'importance pour ce problème.

PROBLÈME : Sur la base de ma distribution, je voudrais calculer la valeur p (la probabilité de voir des valeurs plus élevées) pour une valeur donnée. Par exemple, comme vous pouvez le voir, la valeur p pour 0 approcherait 1 et la valeur p pour les nombres plus élevés tendrait vers 0.

Je ne sais pas si j'ai raison, mais pour déterminer les probabilités, je pense que j'ai besoin d'adapter mes données à une distribution théorique qui est la plus appropriée pour décrire mes données. Je suppose qu'une sorte de test d'ajustement est nécessaire pour déterminer le meilleur modèle.

Existe-t-il un moyen d'implémenter une telle analyse en Python ( ScipyouNumpy )? Pouvez-vous présenter des exemples?

Je vous remercie!

Ajustement de distribution avec somme des erreurs carrées (SSE)

Il s'agit d'une mise à jour et d'une modification de la réponse de Saullo , qui utilise la liste complète des scipy.statsdistributions actuelles et renvoie la distribution avec le moins SSE entre l'histogramme de la distribution et l'histogramme des données.

Exemple de montage

En utilisant l' ensemble de données El Niño destatsmodels , les distributions sont ajustées et l'erreur est déterminée. La distribution avec le moins d'erreur est renvoyée.

Toutes les distributions

Meilleure distribution d'ajustement

Exemple de code

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

Il y a 82 fonctions de distribution implémentées dans SciPy 0.12.0 . Vous pouvez tester l'adaptation de certains d'entre eux à vos données à l'aide de leur fit()méthode . Consultez le code ci-dessous pour plus de détails:

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

Références:

- Distributions d'ajustement, qualité d'ajustement, valeur p. Est-il possible de faire cela avec Scipy (Python)?

- Raccord de distribution avec Scipy

Et voici une liste avec les noms de toutes les fonctions de distribution disponibles dans 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()La méthode mentionnée par @Saullo Castro fournit des estimations du maximum de vraisemblance (MLE). La meilleure distribution pour vos données est celle qui vous donne la plus élevée peut être déterminée de plusieurs manières différentes: comme

1, celui qui vous donne la probabilité log la plus élevée.

2, celui qui vous donne les plus petites valeurs AIC, BIC ou BICc (voir wiki: http://en.wikipedia.org/wiki/Akaike_information_criterion , peut essentiellement être considéré comme une vraisemblance logarithmique ajustée pour le nombre de paramètres, comme une distribution avec plus les paramètres devraient mieux s'adapter)

3, celui qui maximise la probabilité bayésienne postérieure. (voir wiki: http://en.wikipedia.org/wiki/Posterior_probability )

Bien sûr, si vous avez déjà une distribution qui devrait décrire vos données (sur la base des théories de votre domaine particulier) et que vous souhaitez vous en tenir à cela, vous sauterez l'étape d'identification de la distribution la mieux adaptée.

scipyne vient pas avec une fonction pour calculer la vraisemblance logarithmique (bien que la méthode MLE soit fournie), mais le code en dur est facile: voir Les fonctions de densité de probabilité intégrées de `scipy.stat.distributions` sont-elles plus lentes que celles fournies par l'utilisateur?

AFAICU, votre distribution est discrète (et rien que discrète). Par conséquent, le simple fait de compter les fréquences de différentes valeurs et de les normaliser devrait suffire à vos fins. Donc, un exemple pour démontrer ceci:

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

Ainsi, probabilité de voir des valeurs plus élevées que 1simplement (selon la fonction de distribution cumulative complémentaire (ccdf) :

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

Veuillez noter que ccdf est étroitement lié à la fonction de survie (sf) , mais il est également défini avec des distributions discrètes, alors que sf n'est défini que pour les distributions contiguës.

Cela me semble être un problème d'estimation de densité de probabilité.

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

Regarde aussi http://jpktd.blogspot.com/2009/03/using-gaussian-kernel-density.html .

Essayez la distfitbibliothèque.

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

Notez que dans ce cas, tous les points seront significatifs en raison de la distribution uniforme. Vous pouvez filtrer avec dist.y_pred si nécessaire.

Avec OpenTURNS , j'utiliserais les critères BIC pour sélectionner la meilleure distribution qui correspond à ces données. En effet, ce critère ne donne pas trop d'avantages aux distributions qui ont plus de paramètres. En effet, si une distribution a plus de paramètres, il est plus facile pour la distribution ajustée d'être plus proche des données. De plus, le Kolmogorov-Smirnov peut ne pas avoir de sens dans ce cas, car une petite erreur dans les valeurs mesurées aura un impact énorme sur la valeur p.

Pour illustrer le processus, je charge les données El-Nino, qui contiennent 732 mesures de température mensuelles de 1950 à 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

Il est facile d'obtenir les 30 usines de distribution univariées intégrées avec la GetContinuousUniVariateFactoriesméthode statique. Une fois cela fait, la BestModelBICméthode statique renvoie le meilleur modèle et le score BIC correspondant.

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

qui imprime:

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

Afin de comparer graphiquement l'ajustement à l'histogramme, j'utilise les drawPDFméthodes de la meilleure distribution.

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)

Cela produit:

Plus de détails sur ce sujet sont présentés dans le document BestModelBIC . Il serait possible d'inclure la distribution Scipy dans SciPyDistribution ou même avec les distributions ChaosPy avec ChaosPyDistribution , mais je suppose que le script actuel remplit les objectifs les plus pratiques.

Pardonnez-moi si je ne comprends pas votre besoin, mais qu'en est-il de stocker vos données dans un dictionnaire où les clés seraient les nombres entre 0 et 47 et valoriseraient le nombre d'occurrences de leurs clés associées dans votre liste d'origine?
Ainsi, votre probabilité p (x) sera la somme de toutes les valeurs des clés supérieures à x divisée par 30000.