35

いくつかのパラメータ Θ に対する事後 p(Θ|D) が与えられると、次のように定義できます。

後部密度が最も高い領域:

最高後部密度領域は、合計で後部質量の 100(1-α) % を構成する Θ の最も可能性の高い値のセットです。

言い換えれば、与えられた α に対して、以下を満たすp *を探します。

ここに画像の説明を入力

次に、最高事後密度領域をセットとして取得します。

ここに画像の説明を入力

中央信用地域:

上記と同じ表記法を使用すると、信頼できる領域(または間隔) は次のように定義されます。

ここに画像の説明を入力

分布によっては、そのような間隔が多数存在する可能性があります。中心信頼区間は、各尾部に(1-α)/2質量がある信頼区間として定義されます。

計算:

  • 一般的なディストリビューションの場合、ディストリビューションからのサンプルが与えられた場合、Python またはPyMCで上記の 2 つの量を取得するための組み込み機能はありますか?

  • 一般的なパラメトリック分布 (ベータ、ガウスなど) の場合、SciPyまたはstatsmodelsを使用してこれを計算する組み込みまたはライブラリはありますか?

4

7 に答える 7

23

私の理解では、「中央信頼領域」は、信頼区間の計算方法と何ら変わりはありません。必要なのは、cdf関数 atalpha/2と1-alpha/2;の逆関数だけです。scipyこれで呼び出されますppf( パーセンテージポイント関数 ); ガウス事後分布の場合は次のようになります。

>>> from scipy.stats import norm
>>> alpha = .05
>>> l, u = norm.ppf(alpha / 2), norm.ppf(1 - alpha / 2)

事後密度が[l, u]カバーされていることを確認するには:(1-alpha)

>>> norm.cdf(u) - norm.cdf(l)
0.94999999999999996

a=1同様に、 sayとを使用した Beta posterior の場合b=3:

>>> from scipy.stats import beta
>>> l, u = beta.ppf(alpha / 2, a=1, b=3), beta.ppf(1 - alpha / 2, a=1, b=3)

そしてまた：

>>> beta.cdf(u, a=1, b=3) - beta.cdf(l, a=1, b=3)
0.94999999999999996

ここでは、scipy に含まれているパラメトリック分布を確認できます。ppfそして、それらすべてに機能があると思います。

事後密度が最も高い領域に関しては、pdf機能が必ずしも可逆的ではないため、より注意が必要です。一般に、そのような領域は接続されていない場合もあります。たとえば、ベータ版の場合(ここa = b = .5で見られるように);

しかし、ガウス分布の場合、「最高事後密度領域」が「中央信頼領域」と一致することは容易にわかります。そして、それはすべての対称単峰分布の場合だと思います（つまり、pdf関数が分布のモードを中心に対称である場合）

一般的なケースで考えられる数値的アプローチは、 の数値積分をp*使用した の値に対する二分探索です。積分が の単調関数であるという事実を利用する;pdfp*

混合ガウスの例を次に示します。

[ 1 ]最初に必要なのは、分析用の pdf 関数です。簡単な混合ガウスの場合：

def mix_norm_pdf(x, loc, scale, weight):
    from scipy.stats import norm
    return np.dot(weight, norm.pdf(x, loc, scale))

たとえば、場所、スケール、および重量の値については、次のように

loc    = np.array([-1, 3])   # mean values
scale  = np.array([.5, .8])  # standard deviations
weight = np.array([.4, .6])  # mixture probabilities

手をつないで 2 つの素敵なガウス分布が得られます。

ここに画像の説明を入力

[2]p*ここで、上記のpdf関数を統合するためのテスト値を指定p*し、目的の値から二乗誤差を返すエラー関数が必要です1 - alpha。

def errfn( p, alpha, *args):
    from scipy import integrate
    def fn( x ):
        pdf = mix_norm_pdf(x, *args)
        return pdf if pdf > p else 0

    # ideally integration limits should not
    # be hard coded but inferred
    lb, ub = -3, 6 
    prob = integrate.quad(fn, lb, ub)[0]
    return (prob + alpha - 1.0)**2

[ 3 ]ここで、 の特定の値に対してalpha、誤差関数を最小化して を取得できp*ます。

alpha = .05

from scipy.optimize import fmin
p = fmin(errfn, x0=0, args=(alpha, loc, scale, weight))[0]

その結果、p* = 0.0450、および HPD は次のようになります。赤い領域1 - alphaは分布を表し、水平の破線はp*です。

ここに画像の説明を入力

于 2014-03-10T00:19:46.547 に答える
9

PyMC には、hpd を計算する組み込み関数があります。v2.3 では utils にあります。ここのソースを参照してください。線形モデルと HPD の例として

import pymc as pc  
import numpy as np
import matplotlib.pyplot as plt 
## data
np.random.seed(1)
x = np.array(range(0,50))
y = np.random.uniform(low=0.0, high=40.0, size=50)
y = 2*x+y
## plt.scatter(x,y)

## priors
emm = pc.Uniform('m', -100.0, 100.0, value=0)
cee = pc.Uniform('c', -100.0, 100.0, value=0) 

#linear-model
@pc.deterministic(plot=False)
def lin_mod(x=x, cee=cee, emm=emm):
    return emm*x + cee 

#likelihood
llhy = pc.Normal('y', mu=lin_mod, tau=1.0/(10.0**2), value=y, observed=True)

linearModel = pc.Model( [llhy, lin_mod, emm, cee] )
MCMClinear = pc.MCMC( linearModel)
MCMClinear.sample(10000,burn=5000,thin=5)
linear_output=MCMClinear.stats()

## pc.Matplot.plot(MCMClinear)
## print HPD using the trace of each parameter 
print(pc.utils.hpd(MCMClinear.trace('m')[:] , 1.- 0.95))
print(pc.utils.hpd(MCMClinear.trace('c')[:] , 1.- 0.95))

分位点の計算を検討することもできます

print(linear_output['m']['quantiles'])
print(linear_output['c']['quantiles'])

ここで、2.5% から 97.5% の値を取るだけで、95% の中央信頼区間が得られると思います。

于 2014-03-10T17:20:33.727 に答える