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