このページでは、ベイズモデルのパラメータ推定に不可欠な（自分でプログラムを書ける人は別ですが）、MCMCサンプラーの「使い方」について説明します。
なお、Stan以外はほぼ共通のBUGS言語で動作します。BUGS言語の意味などを知りたい方は、「BUGSで学ぶ階層モデリング入門: 個体群のベイズ解析」をご覧下さい。
| WinBUGS | OpenBUGS | JAGS | Stan | |
| OS | Windows（頑張ればMacも） | Windows（頑張ればMacも） | マルチプラットフォーム | マルチプラットフォーム |
| 長所 | car.normal()が使える。 | WinBUGSと同様の操作。 | 多項モデルの計算が速い。エラー個所の特定が容易。 | 収束が速い。エラー個所の特定が容易。空間自己相関を考慮可能。 |
| 短所 | 開発は止まっている。エラーメッセージが理解不能。プログラムは非公開。 | 同左。 | 特になし（個人的には）。 | 離散パラメータを推定できない。1stepあたりの計算速度が遅い。 |
| 公式サイト | The BUGS project | OpenBUGS | JAGS | Stan |
| インストール方法 | 北大の久保先生のページ | 2012年度統計数理研究所共同研究集会「データ解析環境Rの整備と利用」の石倉さんのスライド、BUGSstan勉強会、R stan導入公開版 |
以下の例で共通して使うデータを用意します。
set.seed(1) N <- 100 x1 <- rnorm(N, 0, 1) x2 <- rnorm(N, 0, 1) y <- x1 + rnorm(N, 0, 1)
データの定義
データは、データ名のベクトルを与えます。
data <- c("N", "y", "x1", "x2")
初期値の設定
WinBUGSにまつわるエラーの大半はここが原因です。間違いのないように指定しましょう。
inits <- list(list(alpha = 0, bx1 = 0, bx2 = 0, sigma=70),
list(alpha = -5, bx1 = -5, bx2 = -5, sigma=40),
list(alpha = 5, bx1 = 5, bx2 = 5, sigma=50)
)
監視対象パラメータの定義
WinBUGSは、あらかじめ指定したパラメータの結果しかRに返してくれません。これを、監視対象パラメータといい、事前に指定します。
para <- c("alpha", "bx1", "bx2", "sigma")
BUGS言語でモデルを記述
解析したいモデルを、BUGS言語で記述します。ここでは非常に単純なモデルを扱いますので、R上に直接記述します（ただし、念のためtestmod.txtという名前のテキストファイルに同時に出力しています）。モデルのファイルを別に用意する場合は、model { }の部分だけ書いてファイルに保存して下さい。
modelFilename = "testmod.txt"
cat("
model {
for (i in 1:N) {
y[i] ~ dnorm(mu[i], tau)
mu[i] <- alpha + bx1*x1[i] + bx2*x2[i]
}
alpha ~ dnorm(0.0, 1.0E-3)
bx1 ~ dnorm(0.0, 1.0E-3)
bx2 ~ dnorm(0.0, 1.0E-3)
tau <- pow(sigma, -2)
sigma ~ dunif(0, 100)
}
", fill=TRUE, file=modelFilename)
計算の実行
必要なパッケージを読み込み、解析を実施します。
library(R2WinBUGS)
res <- bugs(data, inits, para,
model=modelFilename,
# ファイルから読み込む場合は、model="testmod.txt",
n.burnin = 500, n.iter=6500, n.chains=3, n.thin=6,
bugs.directory="C:/WinBUGS14/",
# WinBUGSをデフォルト以外の場所にインストール
# している場合のみ指定
working.directory=getwd(),
debug=TRUE
)
結果の確認
推定された係数の値や収束診断などを、要約量を見たり図を書いたりして確認します。
print(res$summary, digits=3) plot(res)
基本的にWinBUGSと同じです。計算の実行関数の一部の引数が異なるだけです。そのため、詳しくは説明しません。
# データの定義
data <- c("N", "y", "x1", "x2")
# 初期値の定義
inits <- list(list(alpha = 0, bx1 = 0, bx2 = 0, sigma=70),
list(alpha = -5, bx1 = -5, bx2 = -5, sigma=40),
list(alpha = 5, bx1 = 5, bx2 = 5, sigma=50)
)
# 監視対象パラメータの定義
para <- c("alpha", "bx1", "bx2", "sigma")
# モデルの記述
modelFilename = "testmod.txt"
cat("
model {
for (i in 1:N) {
y[i] ~ dnorm(mu[i], tau)
mu[i] <- alpha + bx1*x1[i] + bx2*x2[i]
}
alpha ~ dnorm(0.0, 1.0E-3)
bx1 ~ dnorm(0.0, 1.0E-3)
bx2 ~ dnorm(0.0, 1.0E-3)
tau <- pow(sigma, -2)
sigma ~ dunif(0, 100)
}
", fill=TRUE, file=modelFilename)
# 計算の実行
library(R2OpenBUGS)
res2 <- bugs(data, inits, para, model=modelFilename,
n.burnin = 500, n.iter=6500, n.chains=3, n.thin=6,
OpenBUGS.pgm="C:/OpenBUGS322/OpenBUGS.exe",
# プログラムまでフルパスを指定
working.directory=getwd(),
debug=TRUE
)
# 結果の確認
print(res$summary, digits=3)
plot(res)
前者2つと比べると、多少データの指定方法などが異なります。
データの定義
データは、リスト形式で準備して下さい。
list.data <- list(N=N, y=y, x1=x1, x2=x2)
初期値
乱数の種類や種を指定するので、若干WinBUGSと異なります。
inits <- list(alpha = 0, bx1 = 0, bx2 = 0, sigma = 50) inits <- list(inits, inits, inits) inits[[1]]$.RNG.name <- "base::Mersenne-Twister" inits[[1]]$.RNG.seed <- 1 inits[[2]]$.RNG.name <- "base::Mersenne-Twister" inits[[2]]$.RNG.seed <- 12 inits[[3]]$.RNG.name <- "base::Mersenne-Twister" inits[[3]]$.RNG.seed <- 123
監視対象パラメータの定義
ここはWinBUGSと一緒です。
para <- c("alpha", "bx1", "bx2", "sigma")
BUGS言語でモデルを記述
ここもWinBUGSと一緒です。
# モデルの記述
modelFilename = "testmod.txt"
cat("
model {
for (i in 1:N) {
y[i] ~ dnorm(mu[i], tau)
mu[i] <- alpha + bx1*x1[i] + bx2*x2[i]
}
alpha ~ dnorm(0.0, 1.0E-3)
bx1 ~ dnorm(0.0, 1.0E-3)
bx2 ~ dnorm(0.0, 1.0E-3)
tau <- pow(sigma, -2)
sigma ~ dunif(0, 100)
}
", fill=TRUE, file=modelFilename)
計算の実行（runjags）
いくつかパッケージがあります。まず、library(runjags)の場合です。
library(runjags)
res <- autorun.jags(model=modelFilename, data=list.data, inits=inits,
monitor=para,n.chains=3, method="parallel",
thin.sample=TRUE, psrf.target=1.1)
# methodで計算に使うパッケージを指定します。
# psrf.targetに指定したRhat値以下になるようにします。
計算の実行（rjags）
続いて、library(rjags)の場合です。
# 計算回数の設定
n.chains <- 3
n.iter <- 5000
n.update <- 2000
thin <- 5
# パッケージの読み込み
library(rjags)
# 計算
m <- jags.model(file = modelFilename, data = list.data,
inits = inits, n.chain = n.chains
)
update(m, n.update)
x <- coda.samples(m, para,
thin = thin, n.iter = n.iter
)
結果の確認（runjags）
このパッケージの場合はすでに収束していることは明らかですが、推定値などを見ておきます。
# 要約統計量を見る res$summary # MCMCの過程や収束状況などを図で見る par(ask=TRUE) plot(res)
結果の確認（rjags）
推定値や収束指標などを見ておきます。
# 計算結果の出力 res <- data.frame(summary(x)$statistics) ci <- data.frame(summary(x)$quantiles) # 95%信用区間が 0 をまたぐかどうかを計算 res$sig <- abs(sign(ci[, 1]) + sign(ci[, 5])) == 2 # Rhat 値の計算 rhat <- gelman.diag(x)[["psrf"]][, 1] res$Rhat <- rhat # 収束状況の図化 plot(x)
Stanは、モデルの書き方がBUGS言語と多少異なります。
# データの定義
list.data <- list(N=N, y=y, x1=x1, x2=x2)
# 初期値の定義
inits <- list(alpha = 0, bx1 = 0, bx2 = 0, sigma = 50)
inits <- list(inits, inits, inits)
inits[[1]]$.RNG.name <- "base::Mersenne-Twister"
inits[[1]]$.RNG.seed <- 1
inits[[2]]$.RNG.name <- "base::Mersenne-Twister"
inits[[2]]$.RNG.seed <- 12
inits[[3]]$.RNG.name <- "base::Mersenne-Twister"
inits[[3]]$.RNG.seed <- 123
# 監視対象パラメータの定義
para <- c("alpha", "bx1", "bx2", "sigma")
モデルの記述
ここが結構異なります。慣れるまで手間取るかも。
modelfile <- '
data {
int N; //整数はint
real y[N]; //実数はreal
real x1[N];
real x2[N];
}
parameters {
real alpha;
real bx1;
real bx2;
real sigma; //分散パラメータなので、負にならないように範囲を指定
}
model {
for (i in 1:N)
y[i] ~ normal(alpha + bx1*x1[i] + bx2*x2[i], sigma);
//事前分布を指定しない場合、自動的に一様分布が事前分布として与えられます
}
'
計算の実行
以下のようにして、計算を実行します。
library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
fit <- stan(model_code = modelfile, data = list.data, par = para,
iter = 1000, chains = 3, thin=1, init=inits)
結果の確認
結果の見方も、BUGS系列と多少異なります。
print(fit, digits=3)
plot(fit)
#パラメータが多くて部分表示したい場合は
#一度データフレーム型に変換すると便利
res <- as.data.frame(summary(fit))[1:10]
colnames(res) <- c("mean", "se", "sd", "ci2.5", "ci25",
"ci50", "ci75", "ci97.5", "neff", "Rhat")
head(res)
#同じように、パラメータが多くて一部のトレースだけを
#確認したい場合は以下のようにする
traceplot(fit, pars=c("alpha", "bx1"))
空間的自己相関
Hirokiさんによると、Stanでは空間的自己相関も考慮できるようです。