webentwicklung-frage-antwort-db.com.de

Die empirische Verteilung mit Scipy (Python) an die theoretischen anpassen?

EINFÜHRUNG : Ich habe eine Liste von mehr als 30.000 Ganzzahlwerten 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 beispielsweise sehen können, nähert sich der p-Wert für 0 dem Wert 1, und der p-Wert für höhere Zahlen tendiert zu 0.

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

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

Vielen Dank!

110
s_sherly

Verteilungsanpassung mit Quadratsummenfehler (SSE)

Dies ist eine Aktualisierung und Modifikation von Saullos Antwort , die die vollständige Liste der aktuellen scipy.stats Distributionen und gibt die Verteilung mit dem geringsten SSE zwischen dem Histogramm der Verteilung und dem Histogramm der Daten zurück.

Beispiel Fitting

Mit dem El Niño-Datensatz von statsmodels werden die Verteilungen angepasst und der Fehler bestimmt. Die Distribution mit dem geringsten Fehler wird zurückgegeben.

Alle Distributionen

All Fitted Distributions

Best-Fit-Verteilung

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

Es gibt 82 implementierte Verteilungsfunktionen in SciPy 0.12. . Sie können testen, wie einige von ihnen zu Ihren Daten passen, indem Sie die fit() -Methode verwenden. Überprüfen Sie den Code unten für weitere Details:

enter image description here

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 dies mit Scipy (Python) möglich?

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

Die von @Saullo Castro erwähnte Methode fit() liefert Schätzungen der maximalen Wahrscheinlichkeit (Maximum Likelihood Estimation, MLE). Die beste Verteilung für Ihre Daten ist die, die Sie am besten erhalten

1, derjenige, der Ihnen die höchste Log-Wahrscheinlichkeit gibt.

2, derjenige, der Ihnen die kleinsten AIC-, BIC- oder BICc-Werte liefert (siehe wiki: http://en.wikipedia.org/wiki/Akaike_information_criterion ), kann im Grunde genommen als die Log-Wahrscheinlichkeit angepasst an die Anzahl angesehen werden von Parametern, da eine Verteilung mit mehr Parametern voraussichtlich besser passt)

3, derjenige, der 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 des Identifizierens der am besten passenden Verteilung.

scipy enthält keine Funktion zum Berechnen der Log-Wahrscheinlichkeit (obwohl die MLE-Methode bereitgestellt wird), aber der harte Code ist einfach: siehe Ist die Wahrscheinlichkeitsdichte in `scipy.stat eingebaut. Distributionen langsamer als von einem Benutzer angegeben?

10
CT Zhu

AFAICU, Ihre Distribution ist diskret (und nichts als diskret). Daher sollte es für Ihre Zwecke ausreichen, 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()

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

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

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

4
eat

Es klingt für mich wie ein Problem mit der Schätzung der Wahrscheinlichkeitsdichte.

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 .

2
emre

Verzeihen Sie mir, wenn ich Ihren Bedarf nicht verstehe, aber was ist mit dem Speichern Ihrer Daten in einem Wörterbuch, in dem Schlüssel die Zahlen zwischen 0 und 47 sind und die Anzahl der Vorkommen ihrer verwandten Schlüssel in Ihrer ursprünglichen Liste angeben?
Somit ist Ihre Wahrscheinlichkeit p(x) die Summe aller Werte für Schlüssel größer als x geteilt durch 30000.

0
PierrOz