为什么 stan 采样与理论值不匹配？答案

【问题标题】：why stan sampling do not match theoretical values?为什么 stan 采样与理论值不匹配？
【发布时间】：2019-09-14 21:38:00
【问题描述】：

我正在学习 stan，只是尝试了一个非常简单的模型（bernoulli），如下所示，我希望后验采样的平均值为 0.3，因为先验只是均匀分布，但 stan 实际上给出了平均值为 0.33。这是怎么回事？

顺便说一句，我尝试了“优化”，它给出了 0.3，这是我所期望的。

感谢您的帮助！

model_code = "
data {
  int N;
  int y[N];
}

parameters {
  real theta;
}

model {
  theta ~ uniform(0, 1);
  y ~ bernoulli(theta);
}

"

data <- list(
  N = 10, 
  y = c(0, 1, 1, 0, 0, 1, 0, 0, 0, 0)
)

fit = stan(model_code=model_code, data=data, iter=5000)
print(fit)

model = stan_model(model_code=model_code)
mle = optimizing(model, data=data)
print(mle, digits=3)

> print(fit)
...
       mean se_mean   sd  2.5%   25%   50%   75% 97.5% n_eff Rhat
theta  0.33    0.00 0.13  0.11  0.24  0.32  0.42  0.61  6920    1
lp__  -6.56    0.01 0.63 -8.32 -6.71 -6.31 -6.15 -6.11  6813    1


> print(mle, digits=3)
$par
theta 
  0.3 
...

【问题讨论】：

我想通了。基本上我最初的期望是错误的。它应该是 0.33333。这就是贝叶斯的工作原理。

标签： stan rstan

【解决方案1】：

一个问题是参数缺少上下界，应该像这样声明

real<lower = 0, upper = 1> theta;

但只有 10 次观察，后验图的平均值不会接近生成它们的参数。

【讨论】：

【解决方案2】：

这是我在介绍 Stan 课程中使用的示例，以准确解释为什么 Stan 确实符合理论。简短的回答是 beta 分布是倾斜的，所以平均值与众数不匹配。

使用统一的先验，这与theta ~ beta(theta | 1, 1) 相同。以 3 次成功和 7 次失败作为观察结果，分析后验为beta(theta | 4, 8)。众数（最佳）是 3/10，而平均值是 4/12。优化为您提供正确的后验模式估计值 0.30，采样提供正确的后验均值估计值 0.33。

随着成功 (a) 和失败 (b) 观察次数的增加（即 a, b -> 无穷大），后验模式 a / (a + b) 和均值 (a + 1) / (a + b + 2) 接近相同的极限，即经验比 a / (a + b)。在该限制之前，取平均值提供的估计值比取众数具有更低的预期平方误差。

见：https://en.wikipedia.org/wiki/Beta_distribution

【讨论】：