r - 使用 Stan 实现预测后验分布
问题描述
背景
我有一个例子试图在正态测量模型的背景下证明后验预测分布。使用的数据如下:
speed <- c(28, 26, 33, 24, 34, -44, 27, 16, 40, -2, 29, 22, 24, 21, 25, 30, 23, 29, 31, 19, 24, 20, 36, 32, 36, 28, 25, 21, 28, 29, 37, 25, 28, 26, 30, 32, 36, 26, 30, 22, 36, 23, 27, 27, 28, 27, 31, 27, 26, 33, 26, 32, 32, 24, 39, 28, 24, 25, 32, 25, 29, 27, 28, 29, 16, 23)
提供的Stan模型如下:
data{
int<lower=1> n;
vector[n] y;
}
parameters{
real y_mu;
real y_lsd;
}
transformed parameters{
real<lower=0> y_sd;
y_sd = exp(y_lsd);
}
model{
y ~ normal(y_mu, y_sd);
}
generated quantities{
vector[n] y_rep;
for(i in 1:n){
y_rep[i] = normal_rng(y_mu, y_sd);
}
}
然后我们调用以下采样命令:
data.in <- list(y=speed, n=length(speed))
model.fit <- sampling(NMM_PPD, data=data.in)
这个例子表明,正常的测量模型似乎不适用于这些数据。为什么?因为虽然原始数据的平均值和中位数y
几乎位于这些统计数据的中心,这些统计数据是根据从后验预测分布采样的复制数据集计算得出的,但对于最小值、最大值或四分位数范围而言,情况并非如此。此外,与来自后验预测分布的复制数据集的直方图相比,原始数据集的直方图看起来明显不同。如下图所示。
我们首先使用该extract()
函数从model.fit
对象中提取我们复制的数据集:
yrep <- extract(model.fit, pars = "y_rep")[[1]]
直方图:
ppc_hist(speed, yrep[sample(NROW(yrep), 11), ])
意思是:
ppc_stat(speed, yrep)
最大:
ppc_stat(speed, yrep, stat = "max")
其他计算如下:
ppc_stat(speed, yrep, stat = "median")
ppc_stat(speed, yrep, stat = "min")
stat <- function(x) diff(quantile(x, probs = c(0.25, 0.75)))
ppc_stat(speed, yrep, stat = stat)
问题
我现在想拟合以下模型:
(TeX 表示)
$Y_i | \mu, \sigma \sim t_{\nu}(\mu, \sigma)$, $i = 1, \dots, n$ 独立
$\mu \sim N(0, 1000^2)$
$\sigma \sim \text{Half_Cauchy}(0, 5)$
(图像表示)
其中 t 代表随机变量,符号 $\nu$ 代表自由度。
我想尝试不同的 $\nu$ 值,看看哪个值适合对上述统计数据(最大值、平均值、中位数、最小值、分位数)进行建模。
我目前的斯坦代码如下:
data{
int<lower=1> n;
vector[n] y;
}
parameters{
real y_mu;
real y_sd;
real nu;
}
model{
y ~ student_t(nu, y_mu, y_sd);
y_mu ~ normal(0, 1000);
y_sd ~ cauchy(0, 5);
}
generated quantities{
vector[n] y_rep;
for(i in 1:n){
y_rep[i] = student_t_rng(nu, y_mu, y_sd);
}
}
我使用以下代码从模型中抽取样本:
data.in <- list(y=speed, n=length(speed))
model.fit <- sampling(NMM_PPD, data=data.in)
结果如下:
print(model.fit, pars = c("y_mu", "y_sd", "nu"), digits = 5)
所以我们有 nu = 2.56。
但是,我不确定我是否正确地解决了这个问题。这就是我们获得nu
最适合模型的值的方式吗?
我一直在研究预测后验分布的其他 Stan 实现,但我仍然不能 100% 确定我是否正确实现了这一点。
https://magesblog.com/post/2015-05-19-postterior-predictive-output-with-stan/
https://pdfs.semanticscholar.org/4e97/66371e7572609594a4f68fc74b7c6fe70767.pdf
https://magesblog.com/post/2015-05-19-postterior-predictive-output-with-stan/
解决方案
推荐阅读
- ios - iOS在块内崩溃
- hive - 小数点后的 Hive 数字
- azure-devops - Azure DevOps 分支策略:如何解决从 id 到用户的必需审查?
- python - python3.x importing gi... gi包的路径
- python - 如何在 Keras 中为多个序列 LSTM 准备输入数据?
- python - Tesseract OCR - 指定模式
- scala - 隐式添加到集合
- recaptcha - 如何强制 Google reCaptcha 显示挑战图像?
- python - Python在Windows中更改目录
- batch-file - 通过批处理文件比较两个文件夹的内容而不是两个文件的内容