Ajuster une distribution empirique à une distribution théorique avec Scipy (Python) ?

Question

INTRODUCTION : Je dispose d'une liste de plus de 30 000 valeurs entières comprises entre 0 et 47 inclus. [0,0,0,0,..,1,1,1,1,...,2,2,2,2,...,47,47,47,...] échantillonné à 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, j'aimerais calculer la valeur p (la probabilité de voir des valeurs plus élevées) pour toute valeur donnée. Par exemple, comme vous pouvez le voir, la valeur p pour 0 serait proche de 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 je dois 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'adéquation est nécessaire pour déterminer le meilleur modèle.

Existe-t-il un moyen de mettre en œuvre une telle analyse en Python ( Scipy o Numpy ) ? Pouvez-vous présenter des exemples ?

Answer 1

Le code suivant est la version de la réponse générale, mais avec des corrections et de la clarté.

import numpy as np
import pandas as pd
import scipy.stats as st
import statsmodels.api as sm
import matplotlib as mpl
import matplotlib.pyplot as plt
import math
import random

mpl.style.use("ggplot")

def danoes_formula(data):
    """
    DANOE'S FORMULA
    https://en.wikipedia.org/wiki/Histogram#Doane's_formula
    """
    N = len(data)
    skewness = st.skew(data)
    sigma_g1 = math.sqrt((6*(N-2))/((N+1)*(N+3)))
    num_bins = 1 + math.log(N,2) + math.log(1+abs(skewness)/sigma_g1,2)
    num_bins = round(num_bins)
    return num_bins

def plot_histogram(data, results, n):
    ## n first distribution of the ranking
    N_DISTRIBUTIONS = {k: results[k] for k in list(results)[:n]}

    ## Histogram of data
    plt.figure(figsize=(10, 5))
    plt.hist(data, density=True, ec='white', color=(63/235, 149/235, 170/235))
    plt.title('HISTOGRAM')
    plt.xlabel('Values')
    plt.ylabel('Frequencies')

    ## Plot n distributions
    for distribution, result in N_DISTRIBUTIONS.items():
        # print(i, distribution)
        sse = result[0]
        arg = result[1]
        loc = result[2]
        scale = result[3]
        x_plot = np.linspace(min(data), max(data), 1000)
        y_plot = distribution.pdf(x_plot, loc=loc, scale=scale, *arg)
        plt.plot(x_plot, y_plot, label=str(distribution)[32:-34] + ": " + str(sse)[0:6], color=(random.uniform(0, 1), random.uniform(0, 1), random.uniform(0, 1)))

    plt.legend(title='DISTRIBUTIONS', bbox_to_anchor=(1.05, 1), loc='upper left')
    plt.show()

def fit_data(data):
    ## st.frechet_r,st.frechet_l: are disbled in current SciPy version
    ## st.levy_stable: a lot of time of estimation parameters
    ALL_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.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.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
    ]

    MY_DISTRIBUTIONS = [st.beta, st.expon, st.norm, st.uniform, st.johnsonsb, st.gennorm, st.gausshyper]

    ## Calculae Histogram
    num_bins = danoes_formula(data)
    frequencies, bin_edges = np.histogram(data, num_bins, density=True)
    central_values = [(bin_edges[i] + bin_edges[i+1])/2 for i in range(len(bin_edges)-1)]

    results = {}
    for distribution in MY_DISTRIBUTIONS:
        ## Get parameters of distribution
        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_values = [distribution.pdf(c, loc=loc, scale=scale, *arg) for c in central_values]

        ## Calculate SSE (sum of squared estimate of errors)
        sse = np.sum(np.power(frequencies - pdf_values, 2.0))

        ## Build results and sort by sse
        results[distribution] = [sse, arg, loc, scale]

    results = {k: results[k] for k in sorted(results, key=results.get)}
    return results

def main():
    ## Import data
    data = pd.Series(sm.datasets.elnino.load_pandas().data.set_index('YEAR').values.ravel())
    results = fit_data(data)
    plot_histogram(data, results, 5)

if __name__ == "__main__":
    main()

Answer 2

Bien que de nombreuses réponses ci-dessus soient tout à fait valables, aucune ne semble répondre complètement à votre question, en particulier à la partie :

Je ne sais pas si j'ai raison, mais pour déterminer les probabilités, je pense que je dois ajuster 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'adéquation est nécessaire pour déterminer le meilleur modèle.

L'approche paramétrique

Il s'agit du processus que vous décrivez, à savoir l'utilisation d'une distribution théorique et l'ajustement des paramètres à vos données.

L'approche non paramétrique

Cependant, il est également possible d'utiliser une approche non paramétrique pour résoudre votre problème, ce qui signifie que vous ne supposez aucune distribution sous-jacente.

En utilisant ce que l'on appelle la fonction de distribution empirique, qui est égale : Fn(x)= SOMME( I[X<=x] ) / n . Donc la proportion de valeurs inférieures à x.

Comme indiqué dans l'une des réponses ci-dessus, ce qui vous intéresse est l'inverse de la FDC (fonction de distribution cumulative), qui est égale à 1-F(x)

Il est possible de démontrer que la fonction de distribution empirique convergera vers la "vraie" FCD qui a généré vos données.

En outre, il est facile de construire un intervalle de confiance à 1 alpha par :

L(X) = max{Fn(x)-en, 0}
U(X) = min{Fn(x)+en, 0}
en = sqrt( (1/2n)*log(2/alpha)

Dans ce cas P( L(X) <= F(X) <= U(X) ) >= 1-alpha pour tout x.

Je suis assez surpris que PierrOz a obtenu 0 voix, alors qu'il s'agit d'une réponse tout à fait valable à la question utilisant une approche non paramétrique pour estimer F(x).

Notez que la question que vous mentionnez de P(X>=x)=0 pour tout x>47 est simplement une préférence personnelle qui pourrait vous amener à choisir l'approche paramétrique plutôt que l'approche non paramétrique. Les deux approches sont cependant des solutions tout à fait valables à votre problème.

Pour plus de détails et de preuves des affirmations ci-dessus, je vous recommande de consulter le site suivant Tout sur les statistiques : A Concise Course in Statistical Inference de Larry A. Wasserman". Il s'agit d'un excellent ouvrage sur l'inférence paramétrique et non paramétrique.

EDIT : Puisque vous avez spécifiquement demandé des exemples en python, il est possible d'utiliser numpy :

import numpy as np

def empirical_cdf(data, x):
    return np.sum(x<=data)/len(data)

def p_value(data, x):
    return 1-empirical_cdf(data, x)

# Generate some data for demonstration purposes
data = np.floor(np.random.uniform(low=0, high=48, size=30000))

print(empirical_cdf(data, 20))
print(p_value(data, 20)) # This is the value you're interested in

Answer 3

La méthode la plus simple que j'ai trouvée consiste à utiliser le module "fitter". pip install fitter . Tout ce que vous avez à faire est d'importer le jeu de données par pandas. Il dispose d'une fonction intégrée pour rechercher les 80 distributions de Scipy et obtenir le meilleur ajustement aux données par différentes méthodes. Exemple :

f = Fitter(height, distributions=['gamma','lognorm', "beta","burr","norm"])
f.fit()
f.summary()

L'auteur a fourni ici une liste de distributions, car l'analyse des 80 distributions peut prendre beaucoup de temps.

f.get_best(method = 'sumsquare_error')

Vous obtiendrez ainsi les 5 meilleures distributions avec leurs critères d'adéquation :

            sumsquare_error    aic          bic       kl_div
chi2             0.000010  1716.234916 -1945.821606     inf
gamma            0.000010  1716.234909 -1945.821606     inf
rayleigh         0.000010  1711.807360 -1945.526026     inf
norm             0.000011  1758.797036 -1934.865211     inf
cauchy           0.000011  1762.735606 -1934.803414     inf

Vous avez également distributions=get_common_distributions() qui contient environ 10 distributions les plus couramment utilisées, et qui les ajuste et les vérifie pour vous.

Il dispose également d'un grand nombre d'autres fonctions, telles que le traçage d'histogrammes, et une documentation complète peut être consultée à l'adresse suivante aquí .

Il s'agit d'un module très sous-estimé pour les scientifiques, les ingénieurs et en général.

Answer 4

Pourquoi ne pas stocker vos données dans un dictionnaire dont les clés seraient les nombres compris entre 0 et 47 et les valeurs le nombre d'occurrences des clés correspondantes dans votre liste initiale ?

Ainsi, votre probabilité p(x) sera la somme de toutes les valeurs des clés supérieures à x divisée par 30000.

Answer 5

Il me semble qu'il s'agit d'un problème d'estimation de la 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:])

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