2

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

4

1 回答 1

3

我不会认为这是特定于 Stan/PyStan 的问题。根据定义，最高密度区间是单个区间，因此不适用于表征多峰分布。Rob Hyndman 的一项出色的工作，计算和绘制最高密度区域，将概念扩展到多模态分布，这已在hdrcde 包下的R 中实现。

至于 Python，在 PyMC Discourse 网站上对此进行了讨论，建议使用hpd_gridOsvaldo Martin 为他的“使用 Python 进行贝叶斯分析”一书所写的函数 ( )。该函数的源代码在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 采样中却非常频繁地出现，尤其是在使用多个链时（这是最佳实践）。如果人们先验地期望多模态，那么通常这会导致某种形式的混合模型，从而消除多模态。如果一个人不期望多模态，但后验仍然表现出它，那么这是一个对结果持怀疑态度的危险信号。

于 2018-12-09T03:23:25.797 回答