Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

How to use stan in rmarkdown

I would like to get the estimated coefficients of a model using rstan in an rnotebook

I have the following stan chunk:

```{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;
}
```

I also have the following data

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

The documentation on github shows rats_fit <- stan(file = 'https://raw.githubusercontent.com/stan-dev/example-models/master/bugs_examples/vol1/rats/rats.stan'), but since I am using a chunk I don't have a file to refer to.

I have tried stan(rats), summary(rats), print(rats), but none of these seem to work.

like image 252
Alex Avatar asked Feb 03 '18 01:02

Alex


People also ask

How do I embed code into rmarkdown?

You can insert an R code chunk either using the RStudio toolbar (the Insert button) or the keyboard shortcut Ctrl + Alt + I ( Cmd + Option + I on macOS).

How do I exclude output in rmarkdown?

You use results="hide" to hide the results/output (but here the code would still be displayed). You use include=FALSE to have the chunk evaluated, but neither the code nor its output displayed.

How do I Preview in rmarkdown?

You can preview your file by using the shortcut shift + ctrl + k on rmarkdown::render("file.


1 Answers

The first RMarkdown chunk calls rats <- rstan::stan_model(model_code=the_text) behind the scenes, so in order to sample from that posterior distribution you need to ultimately do rats_fit <- sampling(rats, data = list()), whose remaining arguments are pretty much the same as for stan. But you do have to call library(rstan) before all that.

like image 159
Ben Goodrich Avatar answered Sep 23 '22 23:09

Ben Goodrich