python - Pystan 后验不确定区间

Question

我看到了on another forum PyStan 没有与 RStan 相同的功能，它们使用 posterior_interval()，但我们可以使用 numpy.percentile()。我目前正在使用 pystan.StanModel.optimizing() PyStan 中的函数以获得最大化后验似然的参数集。我现在还想获得后验结果的外部 95% 置信区间，所以我想知道 numpy.percentile() 函数是否会与 optimizing 一起使用功能？

我尝试找到参数分布的 95% 区间，但这并没有给出围绕结果的良好置信区间。特别是，我不认为它很好，因为当我期望后验呈现多峰分布时，我使用 numpy.percentile() 进行的置信区间位于后验 2-D 高斯补丁内。

我认为 95% 的间隔必须从后部获取。 我是否可以将百分位数函数与优化函数一起使用以获得 95% 置信度的后验结果？

Answer 1

要获得后验估计的界限，需要对后验进行采样，pystan.StanModel.optimizing 不会这样做。相反，使用 pystan.StanModel.sampling从后验生成 MCMC 绘图的方法。

如果您只需要标准置信区间的读数，那么 pystan.StanFit.stansummary() 方法可能就足够了，因为这将打印 2.5%、25%、50%、75% , 以及每个参数的 97.5% 分位数。例如，

fit = sm.sampling(...) # eight schools model
print(fit.stansummary())

Inference for Stan model: anon_model_19a09b474d1901f191444eaf8a6b8ce2.
4 chains, each with iter=10000; warmup=5000; thin=1;  post-warmup
draws per chain=5000, total post-warmup draws=20000.

           mean se_mean     sd   2.5%    25%    50%   75%  97.5%   n_eff   Rhat 
mu         7.98    0.05   5.04   -2.0   4.76   7.91  11.2   18.2   10614    1.0 
tau        6.54    0.08   5.65   0.24   2.49   5.25   8.98  20.65   4552    1.0 
eta[0]     0.39  6.7e-3   0.94  -1.53  -0.23   0.42   1.02   2.18  20000    1.0 
eta[1]   3.3e-4  6.2e-3   0.88  -1.74  -0.58-2.5e-3   0.57   1.75  20000    1.0 
eta[2]     -0.2  6.6e-3   0.93  -2.01  -0.84  -0.22   0.41   1.68  20000    1.0 
eta[3]    -0.03  6.3e-3   0.89   -1.8  -0.61  -0.03   0.56   1.75  20000    1.0 
eta[4]    -0.35  6.7e-3   0.88  -2.04  -0.94  -0.36   0.22   1.44  17344    1.0 
eta[5]    -0.22  6.6e-3    0.9  -1.96  -0.81  -0.24   0.35   1.59  18298    1.0 
eta[6]     0.34  6.8e-3   0.88  -1.43  -0.23   0.36   0.93   2.04  16644    1.0 
eta[7]     0.05  6.6e-3   0.93  -1.77  -0.58   0.04   0.66   1.88  20000    1.0 
theta[0]   11.4    0.07   8.23  -1.83   6.04  10.22  15.42  31.52  13699    1.0 
theta[1]   7.93    0.04   6.21  -4.58   4.09   7.89  11.79  20.48  20000    1.0 
theta[2]   6.17    0.06   7.72 -11.43   2.06   6.65  10.85  20.53  16367    1.0 
theta[3]   7.72    0.05   6.53  -5.29   3.68    7.7  11.75  21.04  20000    1.0 
theta[4]   5.14    0.04   6.35  -9.35   1.44   5.64   9.38  16.49  20000    1.0 
theta[5]   6.11    0.05   6.66  -8.44   2.22   6.44  10.41  18.52  20000    1.0 
theta[6]  10.63    0.05   6.76  -1.25   6.04  10.08  14.51  25.66  20000    1.0 
theta[7]    8.4    0.06   7.85  -7.56   3.89   8.17  12.76   25.3  16584    1.0 
lp__      -4.89    0.04   2.63 -10.79  -6.47  -4.66  -3.01  -0.43   5355    1.0

或者，如果您需要特定的分位数，您可以像您提到的那样使用 numpy.percentile。

但是，正如您正确观察到的那样，这不适用于多峰分布。此案例在 a different answer 中得到解决。 ，但请注意，如果人们先验期望多个模式，混合模型通常用于将模式分离为不同的单峰随机变量。