RStan

stanファイルの準備

次のようにファイルを作る。
ここではstan_file.stanとしておく。

data{
	int N;
	vector[N] vals;
}

parameters{
	real mu;
	real<lower=0> sigma;
}

model{
	vals ~ normal(mu, sigma);
}

データとして、Nとvalsを読み込む造りである。

RStanのインポート

RでRStanをインポートする。

library(rstan)
rstan_options(auto_write=T)
options(mc.cores=parallel::detectCores())

ここで、高速化のためにRDSファイルの再コンパイル節約と並列化を併せて設定した。

データの読み込み

stanファイルではデータとしてNとvalsを読み込むので、それらのリストを次の形式で作る。

sample_data <- list(N=xxxxx, vals=xxxxx)

MCMCの実行

result <- stan(
	file="stan_file.stan"
	data=sample_data
)

必要に応じて、以下も引数に加える。

• iter：乱数生成の繰り返し数（収束へ向かわせる）
• warmup：バーンイン期間（初期値への依存性を緩和する）
• thin：生成する乱数の間引き（自己相関を緩和する）
• chains：何セット初期値から収束へ向かわせてみるか（チェーンごとに分布が異ならないかチェックするため）

結果の出力

ci <- c(0.025, 0.975)
print(result, probs=ci)
traceplot(result, inc_warmup=T)  #バーンイン期間を含めたトレースプロット

サンプル抽出

sample <- rstan::extract(result, permuted=F)  # permuted=Fはサンプルの順番が変更されるのを防ぐため

でサンプルが抽出できる。
このとき、sampleは3次元になっており、iter・chains・parametersに対応する値が格納されている。

（バーンイン期間を経過した）全ての繰り返しについて、第2チェーンのmuを取り出すなら

sample[,"chain:2","mu"]

となる。
これを利用して代表値を計算することができる。