R为什么斯坦采样与理论值不匹配?



我正在学习斯坦,只是尝试了一个非常简单的模型(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 
...

一个问题是参数缺乏下限和上限,应该像这样声明

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

但是只有十个观测值,后验绘制的平均值不会那么接近生成它们的参数。

这是我在介绍斯坦课程中使用的示例,用于解释为什么斯坦确实符合该理论。 简短的回答是 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

最新更新