必要なパッケージのロード
library(rstan) library(bridgesampling)
パッケージがなければ，
install.packages("bridgesampling")
でインストールします．
stanでのコンパイルの効率化のために，以下を指定します（いつものやつ）．
rstan_options(auto_write=TRUE)#コンパイルを保存する設定 options(mc.cores=parallel::detectCores())#マルチコア計算用設定
データの生成
平均0分散1の正規分布にしたがう乱数を1000個生成し、stanコードでパラメタリカバリします． データの生成
set.seed(123) y <- rnorm(1000, mean = 0, sd = 1)
Stanコード
stanでモデルを定義します．コードは，BF比較用に事前分布も含めて全てtarget記法で書きます． stanファイルを使わずに，すべてR内で完結するように，Rコード上でstanモデルを定義します．
# Stan code for parameter recovery m0 <- " data { int<lower=0> N; vector[N] y; } parameters { real mu; real<lower=0> sigma; } model { target += normal_lpdf(mu | 0, 1); // prior for mu target += normal_lpdf(sigma | 0, 1); // prior for sigma target += normal_lpdf(y | mu, sigma); // likelihood } " # 平均位置がズレている m1 <- " data { int<lower=0> N; vector[N] y; } parameters { real<lower=0> sigma; } model { target += normal_lpdf(sigma | 0, 1); // prior for sigma target += normal_lpdf(y | 0.5, sigma); // likelihood } " # m0に比べて事前分布のサポートが広い場合 m2 <- " data { int<lower=0> N; vector[N] y; } parameters { real mu; real<lower=0> sigma; } model { target += normal_lpdf(mu | 0, 10); // prior for mu target += normal_lpdf(sigma | 0, 10); // prior for sigma target += normal_lpdf(y | mu, sigma); // likelihood } "
Prepare data for Stan
stan用にデータを定義します．
data <- list(N = length(y), y =y)
Compile the Stan model
stanモデルをコンパイルします
model0 <- stan_model(model_code = m0) model1 <- stan_model(model_code = m1) model2 <- stan_model(model_code = m2)
Fit the model
サンプリング結果を格納します．
fit0 <- sampling(model0, data = data, iter = 5000, chains = 3) fit1 <- sampling(model1, data = data, iter = 5000, chains = 3) fit2 <- sampling(model2, data = data, iter = 5000, chains = 3)
Bridge sampling
ブリッジサンプリングします
bs0 <- bridge_sampler(fit0) bs1 <- bridge_sampler(fit1) bs2 <- bridge_sampler(fit2)
自由エネルギー
ブリッジサンプリングした結果を使って自由エネルギを計算します
-logml(bs0) -logml(bs1) -logml(bs2)
Compare free energies
ベイズファクタを計算します．直感的に言うと，bf(bs0, bs1)はbs1（後ろ）に対してbs0（前）が，相対的にどれくらいよいのかを表します．
慣習的な目安として10より大きければ，強い支持を与えると考えられている．
bf(bs0, bs1) bf(bs0, bs2)
