Реализация прогнозирующего последующего распределения с использованием Stan

Question

Фон

У меня был пример, который пытался продемонстрировать апостериорное предиктивное распределение в контексте нормальной модели измерения. Данные, которые были использованы следующим образом:

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 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)

Проблема

Теперь я хочу соответствовать следующей модели:

(Представление 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 представляет t случайную величину, а символ $ \ nu $ представляет степени свободы.

Я хочу поэкспериментировать с различными значениями $ \ nu $, чтобы увидеть, какое значение подходит для моделирования вышеупомянутой статистики (максимум, среднее, медиана, минимум, квантиль).

Мой текущий код Stan выглядит следующим образом:

```{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-posterior-predictive-output-with-stan/

https://pdfs.semanticscholar.org/4e97/66371e7572609594a4f68fc74b7c6fe70767.pdf

https://magesblog.com/post/2015-05-19-posterior-predictive-output-with-stan/

Я был бы очень признателен, если бы люди могли уделить время рассмотрению моей работы.