背景

我有一个例子试图在正态测量模型的背景下证明后验预测分布。使用的数据如下：

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模型如下：

```{stan output.var="NMM_PPD"}
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);
  }
}
```

然后我们调用以下采样命令：

```{r}
data.in <- list(y=speed, n=length(speed))
model.fit <- sampling(NMM_PPD, data=data.in)
```

这个例子表明，正常的测量模型似乎不适用于这些数据。为什么？因为虽然原始数据的平均值和中位数y几乎位于这些统计数据的中心，这些统计数据是根据从后验预测分布采样的复制数据集计算得出的，但对于最小值、最大值或四分位数范围而言，情况并非如此。此外，与来自后验预测分布的复制数据集的直方图相比，原始数据集的直方图看起来明显不同。如下图所示。

我们首先使用该extract()函数从model.fit对象中提取我们复制的数据集：

```{r}
yrep <- extract(model.fit, pars = "y_rep")[[1]]
```

直方图：

```{r}
ppc_hist(speed, yrep[sample(NROW(yrep), 11), ])
```

意思是：

```{r}
ppc_stat(speed, yrep)
```

最大：

```{r}
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)

问题

我现在想拟合以下模型：

$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$代表一个$t$随机变量和符号$\nu$表示自由度。

我想尝试不同的值$\nu$查看哪个值适合于对上述统计数据（最大值、平均值、中位数、最小值、分位数）进行建模。

我目前的斯坦代码如下：

```{stan output.var="NMM_PPD"}
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);
  }
}
```

我使用以下代码从模型中抽取样本：

```{r}
data.in <- list(y=speed, n=length(speed))
model.fit <- sampling(NMM_PPD, data=data.in)
```

结果如下：

```{r}
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/

如果人们能花时间审查我的工作，我将不胜感激。