소개 : 0에서 47까지의 정수 값 30,000 개가 넘는 목록이 있습니다 (예 : [0,0,0,0,..,1,1,1,1,...,2,2,2,2,...,47,47,47,...]
일부 연속 분포에서 샘플링). 목록의 값이 반드시 순서대로있는 것은 아니지만 순서는이 문제와 관련이 없습니다.
문제 : 분포에 따라 주어진 값에 대해 p- 값 (더 큰 값을 볼 확률)을 계산하고 싶습니다. 예를 들어, 0에 대한 p- 값이 1에 가까워지고 높은 숫자에 대한 p- 값이 0에 가까워지는 것을 볼 수 있습니다.
나는 옳은지 모르겠지만 확률을 결정하기 위해 내 데이터를 설명하기에 가장 적합한 이론적 분포에 내 데이터를 맞출 필요가 있다고 생각합니다. 최고의 모델을 결정하기 위해서는 어떤 종류의 적합도 테스트가 필요하다고 가정합니다.
파이썬 ( Scipy
또는 Numpy
)으로 그러한 분석을 구현하는 방법이 있습니까? 예를 제시해 주시겠습니까?
감사합니다!
답변
SSE (Sum of Square Error)를 갖는 분포 피팅
이것은 현재 분포 의 전체 목록을 사용 하고 분포 히스토그램과 데이터 히스토그램 사이에서 SSE 가 가장 작은 분포를 반환하는 Saullo의 답변에 대한 업데이트 및 수정 입니다.scipy.stats
피팅 예
의 El Niño 데이터 세트를statsmodels
사용하여 분포가 적합하고 오차가 결정됩니다. 오류가 가장 적은 분포가 반환됩니다.
모든 배포
최적의 분포
예제 코드
%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')
답변
SciPy 0.12.0 에는 82 개의 구현 된 분포 함수가 있습니다. fit()
방법을 사용하여 일부 데이터가 데이터에 얼마나 적합한 지 테스트 할 수 있습니다 . 자세한 내용은 아래 코드를 확인하십시오.
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()
참고 문헌 :
-적합 분포, 적합도, p- 값. Scipy (Python) 로이 작업을 수행 할 수 있습니까?
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()
@Saullo Castro가 언급 한 방법은 최대 가능성 추정치 (MLE)를 제공합니다. 귀하의 데이터에 가장 적합한 분포는 다음과 같은 여러 가지 방법으로 결정될 수 있습니다.
1, 가장 높은 로그 가능성을 제공하는 것.
2, 가장 작은 AIC, BIC 또는 BICc 값을 제공하는 값 (wiki : http://en.wikipedia.org/wiki/Akaike_information_criterion 참조)은 기본적으로 더 많은 분포를 가진 매개 변수 수에 대해 조정 된 로그 가능성으로 볼 수 있습니다 매개 변수가 더 잘 맞을 것으로 예상됩니다)
3, 베이지안 후 확률을 최대화하는 것. (wiki : http://en.wikipedia.org/wiki/Posterior_probability 참조 )
물론 특정 분야의 이론을 기반으로 데이터를 설명하는 분포가 이미 있고이를 고수하려는 경우 가장 적합한 분포를 식별하는 단계를 건너 뜁니다.
scipy
로그 우도를 계산하는 기능이 없지만 (MLE 방법이 제공되지만) 하드 코드 하나는 쉽다 : ‘scipy.stat.distributions’의 내장 확률 밀도 함수가 사용자가 제공 한 것보다 느리다?를 참조하십시오.
답변
AFAICU, 당신의 배포판은 이산 적입니다. 따라서 다른 값의 빈도를 세고 정규화하는 것만으로도 충분합니다. 따라서 이것을 보여주는 예제 :
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()
따라서 상보 누적 누적 분포 함수 (ccdf)1
에 따르면 단순히 보다 높은 값을 볼 확률은 다음과 같습니다.
In []: 1- cdf[1]
Out[]: 0.40000000000000002
있습니다 CCDF가 밀접하게 관련되어 생존 기능 (SF) 하지만, 반면에 또한 불연속 분포로 정의되어 SF는 연속하고 분배를 위해 정의된다.
답변
그것은 확률 밀도 추정 문제처럼 들립니다.
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:])
http://jpktd.blogspot.com/2009/03/using-gaussian-kernel-density.html 도 참조하십시오 .
답변
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()
이 경우 균일 분포로 인해 모든 점이 중요합니다. 필요한 경우 dist.y_pred로 필터링 할 수 있습니다.
답변
함께 OpenTURNS , 나는 맞는 이러한 데이터가 최선의 분포를 선택 BIC 기준을 사용합니다. 이는이 기준이 더 많은 모수를 갖는 분포에 큰 이점을 제공하지 않기 때문입니다. 실제로 분포에 더 많은 매개 변수가있는 경우 적합 된 분포가 데이터에 더 가깝습니다. 또한 Kolmogorov-Smirnov는이 경우 의미가 없을 수 있습니다. 측정 된 값의 작은 오류가 p- 값에 큰 영향을 미치기 때문입니다.
프로세스를 설명하기 위해 1950 년에서 2010 년까지 732 개의 월별 온도 측정 값이 포함 된 El-Nino 데이터를로드합니다.
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
GetContinuousUniVariateFactories
정적 방법으로 30 개의 빌트인 일 변량 분포 팩토리를 쉽게 얻을 수 있습니다. 완료되면 BestModelBIC
정적 메소드는 최상의 모델과 해당 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)
어떤 인쇄 :
Best= Beta(alpha = 1.64258, beta = 2.4348, a = 18.936, b = 29.254)
히스토그램에 대한 적합치를 그래픽으로 비교하기 위해 drawPDF
최상의 분포 방법을 사용합니다 .
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)
이것은 다음을 생성합니다.
이 주제에 대한 자세한 내용은 BestModelBIC 문서에 나와 있습니다. Scipy 배포를 SciPyDistribution 또는 ChaosPyDistribution과 함께 ChaosPy 배포에 포함 시킬 수 는 있지만 현재 스크립트가 가장 실용적인 목적을 달성한다고 생각합니다.