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

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%的中心可信区间。

Answer 2

为了计算HPD，可以利用pymc3，下面是一个例子

import pymc3
from scipy.stats import norm
a = norm.rvs(size=10000)
pymc3.stats.hpd(a)

Answer 3

根据我的理解，“中心可信区”与计算置信区间没有任何不同；你所需要的只是cdf函数在alpha/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和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

here您可以看到scipy中包含的参数分布；我猜它们都有ppf函数

至于最高后密度区，则更为棘手，因为pdf函数不一定是可逆的；一般来说，这样的区域甚至可能是不相连的；例如，在β与a = b = .5的情况下（可以看到here）

但是，在高斯分布的情况下，很容易看到“最高后验密度区域”与“中心可信区域”一致；我认为这是所有对称的单峰分布的情况（即，如果pdf函数围绕分布模式对称）

一般情况下的一种可能的数值方法是使用pdf的numerical integration对p*的值进行二进制搜索；利用积分是p*的单调函数的事实

---

下面是混合高斯的一个例子：

[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]现在，您需要一个错误函数，该函数为p*提供一个测试值，集成上述p*的pdf函数，并从所需值返回平方误差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*。