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

Question

给定一些参数Θ的后p(Θ| D),可以定义以下内容:

最高后部密度区域:

的最高后验密度区域是集合Θ的最可能值,在总构成后部100质量(1-α)%的.

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

然后获得最高后部密度区域作为集合:

中央可信区域:

使用与上述相同的表示法,可信区域(或区间)定义为:

根据分布,可能有许多这样的间隔.中心可信区间定义为每个尾部有(1-α)/ 2质量的可信区间.

计算:

• 对于一般分布,从分布中给出样本,是否有任何内置函数可以在Python或PyMC中获得上述两个量？

• 对于常见的参数分布(例如Beta,Gaussian等),是否有任何内置函数或库可以使用SciPy或statsmodels进行计算？

Answer 1

根据我的理解,"中心可信区域"与置信区间的计算方式没有任何不同; 你需要的是逆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

类似于Beta后验的说法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

在这里你可以看到scipy中包含的参数分布; 我想他们都有ppf功能;

至于最高的后密度区域,它更棘手,因为pdf功能不一定是可逆的; 一般来说,这样的地区甚至可能没有连接; 例如,在Beta的情况下a = b = .5(如这里可以看到);

但是,在高斯分布的情况下,很容易看出"最高后部密度区域"与"中央可信区域"一致; 我认为所有对称单模态分布都是这种情况(即pdf函数是否围绕分布模式对称)

对于一般情况的可能数值方法将超过的值的二进制搜索p*使用数值积分的pdf; 利用积分是单调函数的事实p*;

以下是混合Gaussian的示例:

[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*积分上面的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*.

Answer 2

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

Answer 3

另一个选项(改编自R到Python)取自John K. Kruschke的"做贝叶斯数据分析"一书,内容如下:

from scipy.optimize import fmin
from scipy.stats import *

def HDIofICDF(dist_name, credMass=0.95, **args):
    # freeze distribution with given arguments
    distri = dist_name(**args)
    # initial guess for HDIlowTailPr
    incredMass =  1.0 - credMass

    def intervalWidth(lowTailPr):
        return distri.ppf(credMass + lowTailPr) - distri.ppf(lowTailPr)

    # find lowTailPr that minimizes intervalWidth
    HDIlowTailPr = fmin(intervalWidth, incredMass, ftol=1e-8, disp=False)[0]
    # return interval as array([low, high])
    return distri.ppf([HDIlowTailPr, credMass + HDIlowTailPr])

我们的想法是创建一个函数intervalWidth,它返回从lowTailPr开始并具有credMass质量的区间宽度.intervalWidth函数的最小值是通过使用scipy中的fmin minimizer建立的.

例如:结果:

print HDIofICDF(norm, credMass=0.95, loc=0, scale=1)

是

    [-1.95996398  1.95996398]

传递给HDIofICDF的分发参数的名称必须与scipy中使用的完全相同.

Answer 4

要计算HPD,你可以利用pymc3,这是一个例子

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

Answer 5

我偶然发现了这篇文章试图找到一种从MCMC样本中估算HDI的方法,但没有一个答案对我有用.像aloctavodia一样,我在"做贝叶斯数据分析"一书中改编了一个R例子.我需要从MCMC样本计算95%的HDI.这是我的解决方案:

import numpy as np
def HDI_from_MCMC(posterior_samples, credible_mass):
    # Computes highest density interval from a sample of representative values,
    # estimated as the shortest credible interval
    # Takes Arguments posterior_samples (samples from posterior) and credible mass (normally .95)
    sorted_points = sorted(posterior_samples)
    ciIdxInc = np.ceil(credible_mass * len(sorted_points)).astype('int')
    nCIs = len(sorted_points) - ciIdxInc
    ciWidth = [0]*nCIs
    for i in range(0, nCIs):
    ciWidth[i] = sorted_points[i + ciIdxInc] - sorted_points[i]
    HDImin = sorted_points[ciWidth.index(min(ciWidth))]
    HDImax = sorted_points[ciWidth.index(min(ciWidth))+ciIdxInc]
    return(HDImin, HDImax)

上面的方法根据我的数据给出了逻辑答案!