二番目に簡単な rstan コード

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 <- stan(model_code=stancode, data=datastan, iter=1000, chain=4)
traceplot(fit, ask=T)
print(fit, digit=2)

          mean se_mean   sd     2.5%      25%      50%      75%    97.5% n_eff Rhat
alpha    11.60    0.06 1.11     9.29    10.90    11.60    12.33    13.73   396 1.01
beta      0.77    0.00 0.02     0.73     0.76     0.77     0.79     0.82   404 1.01
s         7.01    0.01 0.15     6.71     6.90     7.01     7.10     7.32   438 1.00
lp__  -2447.58    0.06 1.26 -2450.88 -2448.09 -2447.25 -2446.70 -2446.22   451 1.00