最高后密度区和中心可信区

Question

给定一些参数上的后验p（D），可以define如下：

最高后密度区：

最大后密度区是Θ的一组最可能值，总共构成后肿块的100（1-α）%。

换句话说，对于给定的α，我们寻找满足以下条件的ap*

然后得到最高后密度区域作为集合：

中部可信区域：

使用与上述相同的符号，可信区域（或区间）定义为：

根据分布情况，可能有许多这样的间隔。中心可信区间定义为每个尾部上都有质量（1-α）/2的可信区间。

计算：

• 对于一般的发行版，给定发行版的样本，是否有内置的Python或PyMC来获取上面的两个数量？

• 对于常见的参数分布（如Beta、Gaussian等），是否有内置的或库可以使用SciPy或statsmodels来计算？

Answer 1

PyMC有一个用于计算hpd的内置函数。在2.3版本中是在utils中。请参阅源here。作为线性模型的一个例子，它是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%的中心可信区间。