后向分布Pystan的最高密度间隔(HDI)

Question

我看到，在Pystan，HDI功能可以用来提供95%可信的间隔围绕后分布。然而，他们说，这将只适用于单峰分布。如果我的模型可能有一个多模分布(最多4个峰值)，我有没有办法在Pystan找到HDI？谢谢!

Answer 1

我不认为这是斯坦/派斯坦特有的问题。根据定义，最高密度区间是一个单一区间，因此不适合描述多峰分布。Rob，https://www.jstor.org/stable/2684423做了一项很好的工作，它将这个概念扩展到多模式分布，这已经在hdrcde包下的R中实现了。

至于Python，有PyMC语篇网站对此的探讨，建议使用Osvaldo为他的“”一书编写的函数(hpd_grid)。该函数的源在hpd.py文件中，95%的区域将用于以下

from hpd import hpd_grid

intervals, x, y, modes = hpd_grid(samples, alpha=0.05)

其中，samples是一个参数的后验样本，而intervals是表示密度最高的区域的元组列表。

图例

下面是一个使用假多模数据的示例图。

import numpy as np
from matplotlib import pyplot as plt
from hpd import hpd_grid

# include two modes
samples = np.random.normal(loc=[-4,4], size=(1000, 2)).flatten()

# compute high density regions
hpd_mu, x_mu, y_mu, modes_mu = hpd_grid(samples)

plt.figure(figsize=(8,6))

# raw data
plt.hist(samples), density=True, bins=29, alpha=0.5)

# estimated distribution
plt.plot(x_mu, y_mu)

# high density intervals
for (x0, x1) in hpd_mu:
    plt.hlines(y=0, xmin=x0, xmax=x1, linewidth=5)
    plt.axvline(x=x0, color='grey', linestyle='--', linewidth=1)
    plt.axvline(x=x1, color='grey', linestyle='--', linewidth=1)

# modes
for xm in modes_mu:
    plt.axvline(x=xm, color='r')

plt.show()

​

​

警告注意事项

应该注意的是，在适当建模的参数上的多模后验分布通常是罕见的，但在非收敛的MCMC抽样中却经常出现，特别是当使用多个链时(这是最好的做法)。如果一个人先验地期望多模，那么通常这会导致某种形式的混合模型，从而消除多模。如果一个人不期待多模，但无论如何，后座展示它，那么这是一个对结果表示怀疑的危险信号。