RStan の使い方

こちらのサンプルを実行して見ました。
Example
the Eight Schools examples というものです。

フォルダー構造

$ tree
.
├── ex_school.R
└── schools.stan

stan のプログラム

schools.stan

data {
  int<lower=0> J;          // number of schools 
  real y[J];               // estimated treatment effects
  real<lower=0> sigma[J];  // s.e. of effect estimates 
}
parameters {
  real mu; 
  real<lower=0> tau;
  vector[J] eta;
}
transformed parameters {
  vector[J] theta;
  theta = mu + tau * eta;
}
model {
  target += normal_lpdf(eta | 0, 1);
  target += normal_lpdf(y | theta, sigma);
}

R のプログラム

ライブラリーのインストール

install.packages('rstan',dependencies=TRUE)

ex_school.R

schools_data <- list(
  J = 8,
  y = c(28,  8, -3,  7, -1,  1, 18, 12),
  sigma = c(15, 10, 16, 11,  9, 11, 10, 18)
)

library(rstan)
fit1 <- stan(
  file = "schools.stan",  # Stan program
  data = schools_data,    # named list of data
  chains = 4,             # number of Markov chains
  warmup = 1000,          # number of warmup iterations per chain
  iter = 2000,            # total number of iterations per chain
  cores = 1,              # number of cores (could use one per chain)
  refresh = 0             # no progress shown
  )

print(fit1, pars=c("theta", "mu", "tau", "lp__"), probs=c(.1,.5,.9))

plot(fit1)

実行結果

$ Rscript ex_school.R 
 要求されたパッケージ StanHeaders をロード中です 
 要求されたパッケージ ggplot2 をロード中です 
rstan (Version 2.21.7, GitRev: 2e1f913d3ca3)
For execution on a local, multicore CPU with excess RAM we recommend calling
options(mc.cores = parallel::detectCores()).
To avoid recompilation of unchanged Stan programs, we recommend calling
rstan_options(auto_write = TRUE)
Inference for Stan model: schools.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

           mean se_mean   sd    10%    50%    90% n_eff Rhat
theta[1]  11.32    0.14 8.14   2.17  10.31  21.76  3625    1
theta[2]   7.99    0.09 6.37   0.11   8.09  15.53  5526    1
theta[3]   6.26    0.11 7.87  -3.59   6.67  15.35  4702    1
theta[4]   7.63    0.09 6.44  -0.31   7.73  15.45  5250    1
theta[5]   5.15    0.09 6.37  -3.34   5.66  12.81  4756    1
theta[6]   6.19    0.09 6.46  -1.95   6.60  13.86  5437    1
theta[7]  10.67    0.10 6.71   2.67  10.05  19.64  4352    1
theta[8]   8.54    0.12 7.72  -0.28   8.27  17.64  4300    1
mu         8.00    0.09 5.06   1.83   8.09  14.32  3295    1
tau        6.32    0.12 5.21   0.94   5.16  13.06  2048    1
lp__     -39.61    0.07 2.63 -43.10 -39.40 -36.46  1270    1

Samples were drawn using NUTS(diag_e) at Thu Oct 27 15:57:53 2022.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).
'pars' not specified. Showing first 10 parameters by default.
ci_level: 0.8 (80% intervals)
outer_level: 0.95 (95% intervals)

作成された PDF ファイルを表示

evince Rplots.pdf

次のバージョンで確認しました。

$ Rscript --version
Rscript (R) version 4.2.1 (2022-06-23)