这篇来自 Rbloggers 的博客文章描述了如何编写一个简单的三部分正态混合模型，该模型具有已知的混合系数、均值和标准差。虽然它详细描述了 R 中许多可用的贝叶斯采样工具的过程，但作者报告说她/他无法使用rstan. 是否可以在 中构建此模型rstan？如果是这样，怎么办？
该模型是具有已知均值和标准差的法线的混合：
用于法线简单混合的 Rstan Stan 模型
机器算法验证
贝叶斯
马尔可夫链蒙特卡罗
斯坦
2022-03-16 06:14:54
1个回答
因为这个分布的所有参数都是已知的，我们只想从这个分布中抽取样本，所以对模型进行编码rstan很简单。请注意，到目前为止，就我花在编码上的时间（15 分钟）而言，这是从这个特定模型中采样效率最低的方法之一。原帖的作者是正确的，他指出从该模型中采样的最简单方法是sample创造性地使用该函数。
library(rstan)
mix_model <- "
data{
int J;
vector<lower=0>[J] weights;
vector[J] means;
vector<lower=0>[J] sdevs;
}
transformed data{
vector[J] ln_weights;
ln_weights <- log(weights);
}
parameters{
real y;
}
model{
vector[J] probs;
for(j in 1:J){
probs[j] <- exp(ln_weights[j]+normal_log(y,means[j],sdevs[j]));
}
increment_log_prob(log(sum(probs)));
}
"
mixdata <- list(J=3, weights=c(0.3,0.4,0.3),means=c(-3,2,10),sdevs=c(2,1,4))
testfit <- stan(model_code=mix_model, data=mixdata, iter=10)
fit <- stan(fit=testfit, data=mixdata, iter=25000, chains=5)
我采取了将混合物的每个参数作为数据读取的步骤，以便“几个已知法线的总和”模型很容易扩展到任意数量的混合物成分的情况。
将混合权重转换为对数比例是transformed data因为在这个模型中它是已知的。在那里而不是在model块中转换它意味着我们只是读取存储的值，而不是在每次迭代时重新计算日志。
这个模型中我唯一不满意的部分是每个混合成分的对数概率循环。通常，人们更喜欢使用 的本机组合函数，rstan因为它们已经计算出衍生函数，因此您不必使用较慢的 autodiff 例程。另一方面，本例中的组合函数只接受两个参数，而不是 3 个或更多......
y <- extract(fit, "y")[[1]]
plot(density(y))
x <- seq(-10,25,by=0.01)
y1 <- 0.3*dnorm(x, mean=-3,sd=2)
y2 <- 0.4*dnorm(x, mean=2,sd=1)
y3 <- 0.3*dnorm(x, mean=10,sd=4)
lines(x,y1, col="red", lty="dashed")
lines(x,y2, col="red", lty="dashed")
lines(x,y3, col="red", lty="dashed")
从视觉上看，结果似乎是混合物密度的合理近似值。
