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!

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

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

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:

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

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?

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.

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 .

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

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

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.

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:

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.

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.