后验分布 Pystan 的最高密度区间 (HDI)
Highest Density Interval (HDI) for Posterior Distribution Pystan
我看到在 Pystan 中，HDI 函数可用于提供围绕后验分布的 95% 可信区间。然而，他们说它只适用于单峰分布。如果我的模型可能具有多峰分布（最多 4 个峰），有没有办法在 Pystan 中找到 HDI？谢谢！
我认为这不是 Stan/PyStan-specific 问题。根据定义，最高密度区间是单个区间，因此不适合表征多峰分布。 Rob Hyndman 的作品很好，Computing and Graphing Highest Density Regions, that extends the concept to multimodal distributions, and this has been implemented in R under the hdrcde package。
至于 Python，有 a discussion of this on the PyMC Discourse site, where it is recommended to use a function (hpd_grid) that Osvaldo Martin wrote for his "Bayesian Analysis with Python" book. The source for that function is in 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 采样中却非常频繁地出现，尤其是在使用多个链时（这是最佳实践）。如果人们期望多模态 先验 ，那么通常会导致某种形式的混合模型消除多模态。如果一个人不期望多模态，但后验者无论如何都会展示它，那么这是一个对结果持怀疑态度的危险信号。
我看到在 Pystan 中，HDI 函数可用于提供围绕后验分布的 95% 可信区间。然而，他们说它只适用于单峰分布。如果我的模型可能具有多峰分布（最多 4 个峰），有没有办法在 Pystan 中找到 HDI？谢谢！
我认为这不是 Stan/PyStan-specific 问题。根据定义，最高密度区间是单个区间，因此不适合表征多峰分布。 Rob Hyndman 的作品很好，Computing and Graphing Highest Density Regions, that extends the concept to multimodal distributions, and this has been implemented in R under the hdrcde package。
至于 Python，有 a discussion of this on the PyMC Discourse site, where it is recommended to use a function (hpd_grid) that Osvaldo Martin wrote for his "Bayesian Analysis with Python" book. The source for that function is in 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 采样中却非常频繁地出现，尤其是在使用多个链时（这是最佳实践）。如果人们期望多模态 先验 ，那么通常会导致某种形式的混合模型消除多模态。如果一个人不期望多模态，但后验者无论如何都会展示它，那么这是一个对结果持怀疑态度的危险信号。