library(rstan) N <- 1000 x <- rnorm(N, mean = 50, sd = 10) y <- 10 + 0.8 * x + rnorm(N, mean =0, sd = 7) stancode <- ' data{ int<lower=0> N; real x[N]; real y[N]; } parameters { real alpha; real beta; real<lower=0> s; } model{ for(i in 1:N) y[i] ~ normal(alpha + beta * x[i], s); alpha ~ normal(0, 100); beta ~ normal(0, 100); s ~ inv_gamma(0.001, 0.001); } ' datastan <- list(N=N, x=x, y=y) fit <- s