WPROWADZENIE : Mam listę ponad 30 000 wartości całkowitych z przedziału od 0 do 47 włącznie, np. [0,0,0,0,..,1,1,1,1,...,2,2,2,2,...,47,47,47,...]Pobranych z jakiegoś ciągłego rozkładu. Wartości na liście niekoniecznie są w kolejności, ale kolejność nie ma znaczenia dla tego problemu.

PROBLEM : Na podstawie mojego rozkładu chciałbym obliczyć wartość p (prawdopodobieństwo zobaczenia większych wartości) dla dowolnej podanej wartości. Na przykład, jak widać, wartość p dla 0 zbliżałaby się do 1, a wartość p dla wyższych liczb dążyłaby do 0.

Nie wiem, czy mam rację, ale myślę, że aby określić prawdopodobieństwa, muszę dopasować moje dane do teoretycznego rozkładu, który jest najbardziej odpowiedni do opisania moich danych. Zakładam, że do określenia najlepszego modelu potrzebny jest jakiś test dopasowania.

Czy istnieje sposób na zaimplementowanie takiej analizy w Pythonie ( Scipylub Numpy)? Czy mógłbyś przedstawić jakieś przykłady?

Dziękuję Ci!

Dopasowanie rozkładu z sumą błędu kwadratowego (SSE)

Jest to aktualizacja i modyfikacja odpowiedzi Saullo , która wykorzystuje pełną listę bieżących scipy.statsrozkładów i zwraca rozkład z najmniejszą liczbą SSE między histogramem dystrybucji a histogramem danych.

Przykładowe dopasowanie

Korzystając ze zbioru danych El Niño zstatsmodels , rozkłady są zgodne, a błąd jest określany. Zwracana jest dystrybucja z najmniejszą liczbą błędów.

Wszystkie dystrybucje

Dystrybucja najlepszego dopasowania

Przykładowy kod

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

W SciPy 0.12.0 zaimplementowano 82 funkcje dystrybucji . Możesz sprawdzić, jak niektóre z nich pasują do Twoich danych, używając ich fit()metody . Sprawdź poniższy kod, aby uzyskać więcej informacji:

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

Bibliografia:

- Rozkłady dopasowania, dobroć dopasowania, wartość p. Czy można to zrobić w Scipy (Python)?

- Oprawa rozprowadzająca ze Scipy

A tutaj lista z nazwami wszystkich funkcji dystrybucyjnych dostępnych w 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'] 

fit()metoda wspomniana przez @Saullo Castro zapewnia oszacowanie maksymalnego prawdopodobieństwa (MLE). Najlepszy rozkład danych to ten, który daje najwyższy, można określić na kilka różnych sposobów: na przykład

1, ten, który daje największe prawdopodobieństwo logowania.

2, ten, który daje najmniejsze wartości AIC, BIC lub BICc (patrz wiki: http://en.wikipedia.org/wiki/Akaike_information_criterion , zasadniczo można go postrzegać jako prawdopodobieństwo dziennika dostosowane do liczby parametrów, jako dystrybucję z większą parametry powinny być lepiej dopasowane)

3, ten, który maksymalizuje późniejsze prawdopodobieństwo bayesowskie. (patrz wiki: http://en.wikipedia.org/wiki/Posterior_probability )

Oczywiście, jeśli masz już rozkład, który powinien opisywać twoje dane (w oparciu o teorie z twojej konkretnej dziedziny) i chcesz się tego trzymać, pominiesz krok identyfikacji najlepiej dopasowanego rozkładu.

scipynie zawiera funkcji do obliczania prawdopodobieństwa logów (chociaż zapewniona jest metoda MLE), ale twardy kod jest łatwy: zobacz Czy wbudowane funkcje gęstości prawdopodobieństwa w `scipy.stat.distributions` są wolniejsze niż te podane przez użytkownika?

AFAICU, twoja dystrybucja jest dyskretna (i tylko dyskretna). Dlatego samo policzenie częstotliwości różnych wartości i ich normalizacja powinno wystarczyć do twoich celów. A więc przykład, aby to zademonstrować:

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

Zatem prawdopodobieństwo zobaczenia wartości wyższych niż 1jest po prostu (zgodnie z komplementarną funkcją dystrybucji skumulowanej (ccdf) :

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

Proszę to zanotować ccdf jest ściśle powiązany z funkcją przetrwania (sf) , ale jest również definiowany za pomocą dystrybucji dyskretnych, podczas gdy sf jest definiowany tylko dla ciągłych dystrybucji.

Brzmi to jak problem z oszacowaniem gęstości prawdopodobieństwa.

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

Zobacz także http://jpktd.blogspot.com/2009/03/using-gaussian-kernel-density.html .

Spróbuj distfit bibliotekę.

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

Zauważ, że w tym przypadku wszystkie punkty będą znaczące ze względu na równomierny rozkład. W razie potrzeby możesz filtrować za pomocą dist.y_pred.

W OpenTURNS użyłbym kryteriów BIC, aby wybrać najlepszą dystrybucję, która pasuje do takich danych. Dzieje się tak, ponieważ kryteria te nie dają zbyt dużej przewagi rozkładom, które mają więcej parametrów. Rzeczywiście, jeśli rozkład ma więcej parametrów, łatwiej jest, aby dopasowany rozkład był bliżej danych. Co więcej, Kołmogorov-Smirnov może nie mieć sensu w tym przypadku, ponieważ mały błąd w zmierzonych wartościach będzie miał ogromny wpływ na wartość p.

Aby zilustrować proces ładuję dane El-Nino, które zawierają 732 miesięczne pomiary temperatury od 1950 do 2010 roku:

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

Łatwo jest uzyskać 30 wbudowanych jednozmiennych fabryk dystrybucji GetContinuousUniVariateFactoriesmetodą statyczną. Po zakończeniu BestModelBICmetoda statyczna zwraca najlepszy model i odpowiadający mu wynik BIC.

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

który drukuje:

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

Aby graficznie porównać dopasowanie do histogramu, korzystam z drawPDFmetod najlepszego rozkładu.

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)

To daje:

Więcej szczegółów na ten temat znajduje się w dokumencie BestModelBIC . Byłoby możliwe włączenie dystrybucji Scipy do SciPyDistribution lub nawet z dystrybucjami ChaosPy do ChaosPyDistribution , ale myślę, że obecny skrypt spełnia najbardziej praktyczne cele.

Wybacz mi, jeśli nie rozumiem Twojej potrzeby, ale co z przechowywaniem danych w słowniku, w którym klucze byłyby liczbami od 0 do 47 i wartościami liczby wystąpień powiązanych z nimi kluczy na Twojej oryginalnej liście?
Zatem Twoje prawdopodobieństwo p (x) będzie sumą wszystkich wartości kluczy większych niż x podzieloną przez 30000.