python - 后分布 Pystan 的最高密度区间 (HDI)

Question

我发现在 Pystan 中，HDI 函数可用于提供后验分布周围 95% 的可信区间。然而，他们说它只适用于单峰分布。如果我的模型可能具有多峰分布(最多 4 个峰值)，有没有办法在 Pystan 中找到 HDI？谢谢!

Answer 1

我不认为这是 Stan/PyStan 特定的问题。根据定义，最高密度区间是单个区间，因此不适合表征多峰分布。 Rob Hyndman 的作品很好， Computing and Graphing Highest Density Regions ，将概念扩展到多模态分布，这已在 R 中实现 the hdrcde package .

对于Python，有a discussion of this on the PyMC Discourse site ，其中建议使用 Osvaldo Martin 为其“Python 贝叶斯分析”一书编写的函数 (hpd_grid)。该函数的来源位于 the hpd.py file ，对于 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 采样中确实经常出现，特别是在使用多个链时(这是最佳实践)。如果人们先验地期望多模态，那么通常会导致某种形式的混合模型，从而消除多模态。如果人们不期望多模态，但后验无论如何都表现出它，那么这是一个对结果持怀疑态度的危险信号。