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

Ajustement de la distribution avec la somme des erreurs quadratiques (SSE)

Il s'agit d'une mise à jour et d'une modification de Réponse de Saullo qui utilise la liste complète des scipy.stats distributions et renvoie la distribution ayant le moins de ESS entre l'histogramme de la distribution et l'histogramme des données.

Exemple d'ajustement

L'utilisation de la Jeu de données El Niño de statsmodels les distributions sont ajustées et l'erreur est déterminée. La distribution présentant l'erreur la plus faible est renvoyée.

Toutes les distributions

Distribution la mieux adaptée

Exemple de code

%matplotlib inline

import warnings
import numpy as np
import pandas as pd
import scipy.stats as st
import statsmodels.api as sm
from scipy.stats._continuous_distns import _distn_names
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

    # Best holders
    best_distributions = []

    # Estimate distribution parameters from data
    for ii, distribution in enumerate([d for d in _distn_names if not d in ['levy_stable', 'studentized_range']]):

        print("{:>3} / {:<3}: {}".format( ii+1, len(_distn_names), distribution ))

        distribution = getattr(st, distribution)

        # 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
                best_distributions.append((distribution, params, sse))

        except Exception:
            pass

    return sorted(best_distributions, key=lambda x:x[2])

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, density=True, alpha=0.5, color=list(matplotlib.rcParams['axes.prop_cycle'])[1]['color'])

# Save plot limits
dataYLim = ax.get_ylim()

# Find best fit distribution
best_distibutions = best_fit_distribution(data, 200, ax)
best_dist = best_distibutions[0]

# 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[0], best_dist[1])

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

param_names = (best_dist[0].shapes + ', loc, scale').split(', ') if best_dist[0].shapes else ['loc', 'scale']
param_str = ', '.join(['{}={:0.2f}'.format(k,v) for k,v in zip(param_names, best_dist[1])])
dist_str = '{}({})'.format(best_dist[0].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

Il y a plus de 90 fonctions de distribution implémentées dans SciPy v1.6.0 . Vous pouvez tester l'adéquation de certains d'entre eux à vos données en utilisant leur fit() méthode . Consultez le code ci-dessous pour plus de détails :

import matplotlib.pyplot as plt
import numpy as np
import scipy
import scipy.stats
size = 30000
x = np.arange(size)
y = scipy.int_(np.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)
    params = dist.fit(y)
    arg = params[:-2]
    loc = params[-2]
    scale = params[-1]
    if arg:
        pdf_fitted = dist.pdf(x, *arg, loc=loc, scale=scale) * size
    else:
        pdf_fitted = dist.pdf(x, loc=loc, scale=scale) * size
    plt.plot(pdf_fitted, label=dist_name)
    plt.xlim(0,47)
plt.legend(loc='upper right')
plt.show()

Références :

• Ajustement des distributions, qualité de l'ajustement, valeur p. Est-il possible de faire cela avec Scipy (Python) ?

• Ajustement de la 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']

Answer 3

Vous pouvez essayer le bibliothèque distfit . Si vous avez d'autres questions, faites-le moi savoir, je suis également le développeur de cette bibliothèque open-source.

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

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.

Des informations plus détaillées et des exemples peuvent être trouvés sur le site Web de la Commission européenne. pages de documentation .

Answer 4

fit() 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 le plus grand nombre d'estimations, ce qui peut être déterminé de différentes manières : par exemple

1, celui qui donne la plus grande probabilité logarithmique.

2, celui qui vous donne les valeurs AIC, BIC ou BICc les plus faibles (voir wiki : http://en.wikipedia.org/wiki/Akaike_information_criterion , peut être considéré comme la log-vraisemblance ajustée en fonction du nombre de paramètres, étant donné que la distribution avec plus de paramètres est censée mieux s'ajuster)

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

Bien entendu, si vous disposez déjà d'une distribution qui devrait décrire vos données (sur la base des théories en vigueur dans votre domaine) et que vous souhaitez vous en tenir à cette distribution, vous passerez l'étape de l'identification de la distribution la mieux adaptée.

scipy n'est pas accompagné d'une fonction permettant de calculer la log-vraisemblance (bien que la méthode MLE soit fournie), mais il est facile d'en coder une : voir Les fonctions de densité de probabilité intégrées dans `scipy.stat.distributions` sont-elles plus lentes que celles fournies par l'utilisateur ?

Answer 5

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

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, la probabilité d'observer des valeurs supérieures à 1 est simplement (selon le fonction de distribution cumulative complémentaire (ccdf) :

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

Veuillez noter que ccdf est étroitement lié à 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.