PENDAHULUAN : Saya memiliki daftar lebih dari 30.000 nilai integer mulai dari 0 hingga 47, inklusif, misalnya [0,0,0,0,..,1,1,1,1,...,2,2,2,2,...,47,47,47,...]sampel dari beberapa distribusi berkelanjutan. Nilai dalam daftar belum tentu berurutan, tetapi urutan tidak menjadi masalah untuk masalah ini.

MASALAH : Berdasarkan distribusi saya, saya ingin menghitung nilai-p (kemungkinan melihat nilai yang lebih besar) untuk nilai apa pun. Misalnya, seperti yang Anda lihat, nilai p untuk 0 akan mendekati 1 dan nilai p untuk angka yang lebih tinggi akan cenderung ke 0.

Saya tidak tahu apakah saya benar, tetapi untuk menentukan probabilitas, saya rasa saya perlu menyesuaikan data saya dengan distribusi teoretis yang paling sesuai untuk menggambarkan data saya. Saya berasumsi bahwa diperlukan semacam uji kesesuaian yang baik untuk menentukan model terbaik.

Apakah ada cara untuk mengimplementasikan analisis seperti itu dengan Python ( Scipyatau Numpy)? Bisakah Anda memberikan contoh?

Terima kasih!

Pemasangan Distribusi dengan Sum of Square Error (SSE)

Ini adalah pembaruan dan modifikasi jawaban Saullo , yang menggunakan daftar lengkap dari scipy.statsdistribusi saat ini dan mengembalikan distribusi dengan SSE terkecil antara histogram distribusi dan histogram data.

Contoh Pemasangan

Menggunakan dataset El Niño daristatsmodels , distribusinya sesuai dan kesalahan ditentukan. Distribusi dengan kesalahan paling sedikit dikembalikan.

Semua Distribusi

Distribusi Paling Sesuai

Kode Contoh

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

Ada 82 fungsi distribusi yang diimplementasikan di SciPy 0.12.0 . Anda dapat menguji bagaimana beberapa di antaranya sesuai dengan data Anda menggunakan fit()metode mereka . Cek kode di bawah ini untuk lebih jelasnya:

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

Referensi:

- Distribusi pas, kesesuaian, nilai-p. Apakah mungkin melakukan ini dengan Scipy (Python)?

- Distribusi pas dengan Scipy

Dan di sini daftar dengan nama semua fungsi distribusi yang tersedia di 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()Metode yang disebutkan oleh @Saullo Castro memberikan perkiraan kemungkinan maksimum (MLE). Distribusi terbaik untuk data Anda adalah yang memberikan Anda tertinggi dapat ditentukan dengan beberapa cara berbeda: seperti

1, salah satu yang memberi Anda kemungkinan log tertinggi.

2, yang memberi Anda nilai AIC, BIC, atau BICc terkecil (lihat wiki: http://en.wikipedia.org/wiki/Akaike_information_criterion , pada dasarnya dapat dilihat sebagai kemungkinan log yang disesuaikan dengan jumlah parameter, sebagai distribusi dengan lebih banyak parameter diharapkan lebih cocok)

3, salah satu yang memaksimalkan probabilitas posterior Bayesian. (lihat wiki: http://en.wikipedia.org/wiki/Posterior_probability )

Tentu saja, jika Anda sudah memiliki distribusi yang harus mendeskripsikan data Anda (berdasarkan teori di bidang khusus Anda) dan ingin tetap berpegang pada itu, Anda akan melewatkan langkah mengidentifikasi distribusi yang paling sesuai.

scipytidak dilengkapi dengan fungsi untuk menghitung kemungkinan log (meskipun metode MLE disediakan), tetapi kode sulitnya mudah: lihat Apakah fungsi kepadatan probabilitas build-in dari `scipy.stat.distributions` lebih lambat daripada yang disediakan pengguna?

AFAICU, distribusi Anda terpisah (dan tidak lain adalah diskrit). Oleh karena itu hanya menghitung frekuensi dari nilai yang berbeda dan menormalkannya sudah cukup untuk keperluan Anda. Jadi, contoh untuk menunjukkan ini:

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

Dengan demikian, probabilitas melihat nilai yang lebih tinggi daripada 1hanya (menurut fungsi distribusi kumulatif pelengkap (CCDF) :

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

Harap dicatat bahwa ccdf terkait erat dengan fungsi survival (sf) , tetapi juga didefinisikan dengan distribusi diskrit, sedangkan sf hanya didefinisikan untuk distribusi yang berdekatan.

Coba distfitperpustakaan.

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

Perhatikan bahwa dalam kasus ini, semua poin akan menjadi signifikan karena distribusi seragam. Anda dapat memfilter dengan dist.y_pred jika diperlukan.

Dengan OpenTURNS , saya akan menggunakan kriteria BIC untuk memilih distribusi terbaik yang sesuai dengan data tersebut. Ini karena kriteria ini tidak memberikan banyak keuntungan bagi distribusi yang memiliki parameter lebih banyak. Memang, jika suatu distribusi memiliki lebih banyak parameter, maka distribusi yang dipasang akan lebih mudah untuk lebih dekat dengan data. Selain itu, Kolmogorov-Smirnov mungkin tidak masuk akal dalam kasus ini, karena kesalahan kecil dalam nilai yang diukur akan berdampak besar pada nilai p.

Untuk mengilustrasikan prosesnya, saya memuat data El-Nino, yang berisi 732 pengukuran suhu bulanan dari tahun 1950 hingga 2010:

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

Mudah untuk mendapatkan 30 pabrik distribusi univariat built-in dengan GetContinuousUniVariateFactoriesmetode statis. Setelah selesai, BestModelBICmetode statis mengembalikan model terbaik dan skor BIC yang sesuai.

sample = ot.Sample([[p] for p in data]) # data reshaping
tested_factories = ot.DistributionFactory.GetContinuousUniVariateFactories()
best_model, best_bic = ot.FittingTest.BestModelBIC(sample,
                                                   tested_factories)
print("Best=",best_model)

yang mencetak:

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

Untuk membandingkan secara grafis kesesuaian dengan histogram, saya menggunakan drawPDFmetode distribusi terbaik.

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)

Ini menghasilkan:

Detail lebih lanjut tentang topik ini disajikan di dokumen BestModelBIC . Dimungkinkan untuk menyertakan distribusi Scipy di SciPyDistribution atau bahkan dengan distribusi ChaosPy dengan ChaosPyDistribution , tetapi saya rasa skrip saat ini memenuhi sebagian besar tujuan praktis.

Kedengarannya seperti masalah estimasi kepadatan probabilitas bagi saya.

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

Lihat juga http://jpktd.blogspot.com/2009/03/using-gaussian-kernel-density.html .

Maafkan saya jika saya tidak mengerti kebutuhan Anda, tetapi bagaimana dengan menyimpan data Anda dalam kamus di mana kuncinya adalah angka antara 0 dan 47 dan menilai jumlah kemunculan kunci terkait di daftar asli Anda?
Jadi kemungkinan Anda p (x) adalah jumlah dari semua nilai untuk kunci yang lebih besar dari x dibagi dengan 30000.

Meskipun banyak dari jawaban di atas benar-benar valid, sepertinya tidak ada yang menjawab pertanyaan Anda secara lengkap, khususnya bagian:

Saya tidak tahu apakah saya benar, tetapi untuk menentukan probabilitas, saya rasa saya perlu menyesuaikan data saya dengan distribusi teoretis yang paling sesuai untuk menggambarkan data saya. Saya berasumsi bahwa diperlukan semacam uji kesesuaian yang baik untuk menentukan model terbaik.

Pendekatan parametrik

Ini adalah proses yang Anda gambarkan dengan menggunakan beberapa distribusi teoretis dan menyesuaikan parameter dengan data Anda dan ada beberapa jawaban yang sangat bagus bagaimana melakukan ini.

Pendekatan non-parametrik

Namun, mungkin juga untuk menggunakan pendekatan non-parametrik untuk masalah Anda, yang berarti Anda sama sekali tidak mengasumsikan distribusi yang mendasari.

Dengan menggunakan apa yang disebut fungsi distribusi Empiris yang sama dengan: Fn (x) = SUM (I [X <= x]) / n . Jadi proporsi nilainya di bawah x.

Seperti yang ditunjukkan dalam salah satu jawaban di atas adalah bahwa yang Anda minati adalah CDF terbalik (fungsi distribusi kumulatif), yang sama dengan 1-F (x)

Dapat ditunjukkan bahwa fungsi distribusi empiris akan menyatu dengan CDF 'sebenarnya' apa pun yang menghasilkan data Anda.

Selain itu, sangat mudah untuk membuat interval kepercayaan 1-alfa dengan:

L(X) = max{Fn(x)-en, 0}
U(X) = min{Fn(x)+en, 0}
en = sqrt( (1/2n)*log(2/alpha)

Kemudian P (L (X) <= F (X) <= U (X))> = 1-alpha untuk semua x.

Saya cukup terkejut bahwa jawaban PierrOz memiliki 0 suara, sedangkan itu adalah jawaban yang sepenuhnya valid untuk pertanyaan menggunakan pendekatan non-parametrik untuk memperkirakan F (x).

Perhatikan bahwa masalah yang Anda sebutkan tentang P (X> = x) = 0 untuk setiap x> 47 hanyalah preferensi pribadi yang mungkin mengarahkan Anda untuk memilih pendekatan parametrik di atas pendekatan non-parametrik. Namun kedua pendekatan tersebut merupakan solusi yang sepenuhnya valid untuk masalah Anda.

Untuk detail lebih lanjut dan bukti dari pernyataan di atas, saya akan merekomendasikan untuk melihat 'Semua Statistik: Kursus Singkat dalam Inferensi Statistik oleh Larry A. Wasserman'. Buku bagus tentang inferensi parametrik dan non-parametrik.

EDIT: Karena Anda secara khusus meminta beberapa contoh python, itu dapat dilakukan menggunakan numpy:

import numpy as np

def empirical_cdf(data, x):
    return np.sum(x<=data)/len(data)

def p_value(data, x):
    return 1-empirical_cdf(data, x)

# Generate some data for demonstration purposes
data = np.floor(np.random.uniform(low=0, high=48, size=30000))

print(empirical_cdf(data, 20))
print(p_value(data, 20)) # This is the value you're interested in