为什么 stan 抽样与理论值不符？

Question

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

顺便说一句，我试过 "optimizing" 得到 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 
...

Answer 1

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

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

但是只有十次观察，后验抽取的平均值不会那么接近生成它们的参数。

Answer 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

为什么 stan 抽样与理论值不符？

why stan sampling do not match theoretical values?

stan

rstan