Как использовать stan в rmarkdown

Я хотел бы получить оценочные коэффициенты модели, используя rstan в rnotebook

У меня есть следующий блок stan:

```{stan output.var="rats"}
data {
  int<lower=0> N;
  int<lower=0> T;
  real x[T];
  real y[N,T];
  real xbar;
}
parameters {
  real alpha[N];
  real beta[N];

  real mu_alpha;
  real mu_beta;          // beta.c in original bugs model

  real<lower=0> sigmasq_y;
  real<lower=0> sigmasq_alpha;
  real<lower=0> sigmasq_beta;
}
transformed parameters {
  real<lower=0> sigma_y;       // sigma in original bugs model
  real<lower=0> sigma_alpha;
  real<lower=0> sigma_beta;

  sigma_y = sqrt(sigmasq_y);
  sigma_alpha = sqrt(sigmasq_alpha);
  sigma_beta =  sqrt(sigmasq_beta);
}
model {
  mu_alpha ~ normal(0, 100);
  mu_beta ~ normal(0, 100);
  sigmasq_y ~ inv_gamma(0.001, 0.001);
  sigmasq_alpha ~ inv_gamma(0.001, 0.001);
  sigmasq_beta ~ inv_gamma(0.001, 0.001);
  alpha ~ normal(mu_alpha, sigma_alpha); // vectorized
  beta ~ normal(mu_beta, sigma_beta);  // vectorized
  for (n in 1:N)
    for (t in 1:T) 
      y[n,t] ~ normal(alpha[n] + beta[n] * (x[t] - xbar), sigma_y);

}
generated quantities {
  real alpha0;
  alpha0 = mu_alpha - xbar * mu_beta;
}
```

У меня также есть следующие данные

```{r}
df <- read_delim("https://raw.githubusercontent.com/wiki/stan-dev/rstan/rats.txt",delim = " ")

y <- as.matrix(df)
x <- c(8,15,22,29,36)
xbar <- mean(x)
N <- nrow(y)
T <- ncol(y)
```

документация на github показывает rats_fit <- stan(file = 'https://raw.githubusercontent.com/stan-dev/example-models/master/bugs_examples/vol1/rats/rats.stan'), но поскольку я использую кусок У меня нет файла, на который можно сослаться.

Я пробовал stan(rats), summary(rats), print(rats), но ни один из них не работает.


r rstan r-markdown stan
person Alex    schedule 03.02.2018    source источник

Ответы (2)


arrow_upward
7
arrow_downward

Первый фрагмент RMarkdown вызывает rats <- rstan::stan_model(model_code=the_text) за кулисами, поэтому для выборки из этого апостериорного распределения вам нужно в конечном итоге сделать rats_fit <- sampling(rats, data = list()), остальные аргументы которого в значительной степени такие же, как для stan. Но перед всем этим вам нужно позвонить library(rstan).

person Ben Goodrich    schedule 03.02.2018

arrow_upward
1
arrow_downward

Спасибо! с вашей помощью мне удалось придумать следующее

library(tidyverse)
df <- read_delim("https://raw.githubusercontent.com/wiki/stan-dev/rstan/rats.txt",delim = " ")
y <- as.matrix(df)
x <- c(8,15,22,29,36)
xbar <- mean(x)
N <- nrow(y)
T <- ncol(y)
library(rstan)
options(mc.cores = parallel::detectCores())
rats_fit <- rstan::sampling(rats, 
                     data = list(y,x,xbar,N,T))
rstan::summary(rats_fit)
person Alex    schedule 04.02.2018