Passen Sie mit Scipy (Python) die empirische Verteilung an die theoretische an?

139

EINLEITUNG : Ich habe eine Liste von mehr als 30.000 ganzzahligen Werten im Bereich von 0 bis einschließlich 47, z. B. [0,0,0,0,..,1,1,1,1,...,2,2,2,2,...,47,47,47,...]aus einer kontinuierlichen Verteilung. Die Werte in der Liste sind nicht unbedingt in der richtigen Reihenfolge, aber die Reihenfolge spielt für dieses Problem keine Rolle.

PROBLEM : Basierend auf meiner Verteilung möchte ich den p-Wert (die Wahrscheinlichkeit, größere Werte zu sehen) für einen bestimmten Wert berechnen. Wie Sie sehen können, nähert sich der p-Wert für 0 beispielsweise 1 und der p-Wert für höhere Zahlen tendiert zu 0.

Ich weiß nicht, ob ich Recht habe, aber um Wahrscheinlichkeiten zu bestimmen, muss ich meine Daten an eine theoretische Verteilung anpassen, die für die Beschreibung meiner Daten am besten geeignet ist. Ich gehe davon aus, dass eine Art Fit-Test erforderlich ist, um das beste Modell zu ermitteln.

Gibt es eine Möglichkeit, eine solche Analyse in Python ( Scipyoder Numpy) zu implementieren ? Könnten Sie Beispiele vorstellen?

Danke dir!

s_sherly
quelle
2
Sie haben nur diskrete empirische Werte, möchten aber eine kontinuierliche Verteilung? Verstehe ich das richtig
Michael J. Barber
1
Es scheint unsinnig. Was bedeuten die Zahlen? Messungen mit begrenzter Genauigkeit?
Michael J. Barber
1
Michael, ich erklärte, was die Zahlen in meiner vorherigen Frage darstellen: stackoverflow.com/questions/6615489/…
s_sherly
6
Das sind Zähldaten. Es ist keine kontinuierliche Verteilung.
Michael J. Barber
1
Überprüfen Sie die akzeptierte Antwort auf diese Frage stackoverflow.com/questions/48455018/…
Ahmad Suliman

Antworten:

209

Verteilungsanpassung mit Summe der quadratischen Fehler (SSE)

Dies ist eine Aktualisierung und Änderung der Antwort von Saullo , die die vollständige Liste der aktuellen scipy.statsVerteilungen verwendet und die Verteilung mit der geringsten SSE zwischen dem Histogramm der Verteilung und dem Histogramm der Daten zurückgibt.

Beispielanpassung

Unter Verwendung des El Niño-Datensatzes vonstatsmodels werden die Verteilungen angepasst und Fehler ermittelt. Die Verteilung mit dem geringsten Fehler wird zurückgegeben.

Alle Distributionen

Alle angepassten Distributionen

Best Fit Distribution

Best Fit Distribution

Beispielcode

%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')
tmthydvnprt
quelle
2
Genial. Erwägen Sie die Verwendung density=Trueanstelle von normed=Truein np.histogram(). ^^
Peque
1
@tmthydvnprt Vielleicht könnten Sie die Änderungen in den .plot()Methoden rückgängig machen , um zukünftige Verwirrung zu vermeiden. ^^
Peque
10
So erhalten Sie Distributionsnamen : from scipy.stats._continuous_distns import _distn_names. Sie können dann getattr(scipy.stats, distname)für jeden distnamein _distn_names` so etwas wie verwenden. Nützlich, da die Distributionen mit verschiedenen SciPy-Versionen aktualisiert werden.
Brad Solomon
1
Können Sie bitte erklären, warum dieser Code nur nach der besten Anpassung kontinuierlicher Verteilungen sucht und nicht nach diskreten oder multivariaten Verteilungen. Danke dir.
Adam Schroeder
6
Sehr cool. Ich musste den ax = data.plot(kind='hist', bins=50, normed=True, alpha=0.5, color=list(matplotlib.rcParams['axes.prop_cycle'])[1]['color'])
Farbparameter
147

In SciPy 0.12.0 sind 82 Verteilungsfunktionen implementiert . Sie können testen, wie einige von ihnen mit ihrer fit()Methode zu Ihren Daten passen . Überprüfen Sie den Code unten für weitere Details:

Geben Sie hier die Bildbeschreibung ein

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

Verweise:

- Anpassungsverteilungen, Anpassungsgüte, p-Wert. Ist das mit Scipy (Python) möglich?

- Verteiler mit Scipy

Und hier eine Liste mit den Namen aller in Scipy 0.12.0 (VI) verfügbaren Verteilungsfunktionen:

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'] 
Saullo GP Castro
quelle
7
Was ist, wenn normed = Truebeim Zeichnen des Histogramms? Sie würden mehrfach nicht pdf_fitteddurch das size, nicht wahr?
Aloha
3
Sehen Sie sich diese Antwort an, wenn Sie sehen möchten, wie alle Distributionen aussehen, oder wenn Sie eine Vorstellung davon haben möchten, wie Sie auf alle zugreifen können.
tmthydvnprt
@ SaulloCastro Was bedeuten die 3 Werte in param in der Ausgabe von dist.fit
shaifali Gupta
2
So erhalten Sie Distributionsnamen : from scipy.stats._continuous_distns import _distn_names. Sie können dann getattr(scipy.stats, distname)für jeden distnamein _distn_names` so etwas wie verwenden. Nützlich, da die Distributionen mit verschiedenen SciPy-Versionen aktualisiert werden.
Brad Solomon
1
Ich würde color = 'w' aus dem Code entfernen, sonst wird das Histogramm nicht angezeigt.
Eran
12

fit()Die von @Saullo Castro erwähnte Methode liefert Maximum-Likelihood-Schätzungen (MLE). Die beste Verteilung für Ihre Daten ist diejenige, die Ihnen die höchste gibt. Sie kann auf verschiedene Arten ermittelt werden: z

1, diejenige, die Ihnen die höchste Protokollwahrscheinlichkeit gibt.

2, derjenige, der Ihnen die kleinsten AIC-, BIC- oder BICc-Werte liefert (siehe Wiki: http://en.wikipedia.org/wiki/Akaike_information_criterion) , kann grundsätzlich als Protokollwahrscheinlichkeit angesehen werden, die an die Anzahl der Parameter angepasst ist, als Verteilung mit mehr Es wird erwartet, dass die Parameter besser passen.

3, diejenige, die die Bayes'sche hintere Wahrscheinlichkeit maximiert. (siehe Wiki: http://en.wikipedia.org/wiki/Posterior_probability )

Wenn Sie bereits eine Verteilung haben, die Ihre Daten beschreiben sollte (basierend auf den Theorien in Ihrem speziellen Bereich) und sich daran halten möchten, überspringen Sie natürlich den Schritt der Ermittlung der am besten passenden Verteilung.

scipyEs wird keine Funktion zum Berechnen der Protokollwahrscheinlichkeit mitgeliefert (obwohl die MLE-Methode bereitgestellt wird), aber der Hardcode 1 ist einfach: Siehe Sind die integrierten Wahrscheinlichkeitsdichtefunktionen von "scipy.stat.distributions" langsamer als die von einem Benutzer bereitgestellten?

CT Zhu
quelle
1
Wie würde ich diese Methode auf eine Situation anwenden, in der die Daten bereits zusammengefasst wurden - das heißt, es handelt sich bereits um ein Histogramm, anstatt aus den Daten ein Histogramm zu erstellen?
Pete
@pete, das wäre eine Situation von intervallzensierten Daten, es gibt Maximum-Likelihood-Methode dafür, aber es ist derzeit nicht implementiert inscipy
CT Zhu
Vergessen Sie nicht die Beweise
jtlz2
5

AFAICU, Ihre Verteilung ist diskret (und nichts als diskret). Daher sollte es für Ihre Zwecke ausreichen, nur die Frequenzen verschiedener Werte zu zählen und zu normalisieren. Ein Beispiel, um dies zu demonstrieren:

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

Somit ist die Wahrscheinlichkeit, Werte zu sehen, die höher sind als 1einfach (gemäß der komplementären kumulativen Verteilungsfunktion (ccdf) :

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

Bitte beachten Sie, dass ccdf eng mit der Überlebensfunktion (sf) verwandt ist , aber auch mit diskreten Verteilungen definiert ist , während sf nur für zusammenhängende Verteilungen definiert ist.

Essen
quelle
2

Es klingt für mich nach einem Problem der Wahrscheinlichkeitsdichteschätzung.

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

Siehe auch http://jpktd.blogspot.com/2009/03/using-gaussian-kernel-density.html .

emre
quelle
1
Für zukünftige Leser: Diese Lösung (oder zumindest die Idee) bietet die einfachste Antwort auf die Fragen des OP („Was ist der p-Wert?“). Es wäre interessant zu wissen, wie sich dies im Vergleich zu einigen der komplizierteren Methoden verhält, die passen eine bekannte Verbreitung.
Greg
Funktionieren Gaußsche Kernel-Regressionen für alle Distributionen?
@mikey In der Regel funktionieren keine Regressionen für alle Distributionen. Sie sind aber nicht schlecht
TheEnvironmentalist
2

Probieren Sie die distfitBibliothek aus.

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

Geben Sie hier die Bildbeschreibung ein

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

Beste Passform

Beachten Sie, dass in diesem Fall alle Punkte aufgrund der gleichmäßigen Verteilung von Bedeutung sind. Sie können bei Bedarf mit dist.y_pred filtern.

erdogant
quelle
1

Mit OpenTURNS würde ich die BIC-Kriterien verwenden, um die beste Verteilung auszuwählen, die zu solchen Daten passt. Dies liegt daran, dass dieses Kriterium den Verteilungen mit mehr Parametern keinen allzu großen Vorteil verschafft. Wenn eine Verteilung mehr Parameter enthält, ist es für die angepasste Verteilung einfacher, näher an den Daten zu sein. Darüber hinaus ist Kolmogorov-Smirnov in diesem Fall möglicherweise nicht sinnvoll, da ein kleiner Fehler in den gemessenen Werten einen großen Einfluss auf den p-Wert hat.

Um den Prozess zu veranschaulichen, lade ich die El-Nino-Daten, die 732 monatliche Temperaturmessungen von 1950 bis 2010 enthalten:

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

Mit der GetContinuousUniVariateFactoriesstatischen Methode ist es einfach, die 30 integrierten univariaten Verteilungsfabriken abzurufen. Anschließend gibt die BestModelBICstatische Methode das beste Modell und den entsprechenden BIC-Score zurück.

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

welche druckt:

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

Um die Anpassung grafisch mit dem Histogramm zu vergleichen, verwende ich die drawPDFMethoden der besten Verteilung.

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)

Dies erzeugt:

Beta passt zu den El-Nino-Temperaturen

Weitere Details zu diesem Thema finden Sie im BestModelBIC- Dokument. Es wäre möglich, die Scipy-Distribution in die SciPyDistribution oder sogar in ChaosPy-Distributionen mit ChaosPyDistribution aufzunehmen , aber ich denke, dass das aktuelle Skript die meisten praktischen Zwecke erfüllt.

Michael Baudin
quelle
2
Sie sollten wahrscheinlich ein Interesse anmelden?
jtlz2
0

Verzeihen Sie mir, wenn ich Ihre Bedürfnisse nicht verstehe, aber wie wäre es, wenn Sie Ihre Daten in einem Wörterbuch speichern, in dem Schlüssel die Zahlen zwischen 0 und 47 sind und die Anzahl der Vorkommen der zugehörigen Schlüssel in Ihrer ursprünglichen Liste bewertet werden?
Somit ist Ihre Wahrscheinlichkeit p (x) die Summe aller Werte für Schlüssel größer als x geteilt durch 30000.

PierrOz
quelle
In diesem Fall ist p (x) für jeden Wert größer als 47 gleich (gleich 0). Ich benötige eine kontinuierliche Wahrscheinlichkeitsverteilung.
s_sherly
2
@s_sherly - Es wäre wahrscheinlich eine gute Sache, wenn Sie bearbeiten können und Ihre Frage besser verdeutlichen, wie in der Tat die „die Wahrscheinlichkeit für größere Werte zu sehen“ - wie Sie es nennen - IS Null für Werte , die über dem höchsten Wert im Pool .
Mac