我在另一个论坛上看到PyStan的功能与RStan不同,他们使用posterior_interval()
,但我们可以使用numpy.percentile()
。我目前正在使用PyStan中的pystan.StanModel.optimizing()
函数来获得最大后验似然的参数集。我现在也想获得后验结果的外部95%置信区间,所以我想知道numpy.percentile()
函数是否会与optimizing
函数一起使用?
我试图找到参数分布的95%区间,但这并没有给出一个很好的置信区间。特别是,我认为这不好,因为当我期望后验呈现多峰分布时,我使用numpy.percentile()
进行的置信区间在后验2-D高斯贴片内。
我认为95%的间隔必须从后部开始我会使用百分比函数和优化函数来获得95%置信的后验结果吗
要获得后验估计的边界,需要对后验进行采样,而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
。
然而,正如您正确地观察到的,这不适用于多模式分布。这种情况在不同的答案中得到了解决,但请注意,如果预期有多个模式先验,则通常使用混合模型将模式分离为不同的单峰随机变量。
您可以直接从pystan.stansummary
:中检索您想要的百分位数
percentiles = (0.025, 0.25, 0.5, 0.75, 0.975) # edit these at will
pystan.stansummary(fit=your_fit, probs=percentiles, digits_summary=2)
这应该很好用。