「研究ソフト/MCMCサンプラー」の編集履歴（バックアップ）一覧に戻る
研究ソフト/MCMCサンプラー - (2016/07/25 (月) 15:27:52) のソース
* MCMCサンプラー
このページでは、ベイズモデルのパラメータ推定に不可欠な（自分でプログラムを書ける人は別ですが）、MCMCサンプラーの「使い方」について説明します。
| 名称 | WinBUGS | OpenBUGS | JAGS | Stan |
| 利用可能OS | Windows（頑張ればMacも）| Windows（頑張ればMacも）| マルチプラットフォーム | マルチプラットフォーム |
| 長所 | car.normal()が使える。 | WinBUGSとほぼ同じ操作。 | 多項モデルの計算が速い。マルチプラットフォーム。エラー個所の特定が容易。 | 収束が速い。マルチプラットフォーム。エラー個所の特定が容易。 |
| 短所 | 開発は止まっている。エラーメッセージが理解不能。プログラムは非公開。 | WinBUGSと同じ。 | 推定の実行が複数段階にわたり面倒くさい。結果のオブジェクトがそのままでは使いにくい。| 離散パラメータを推定できない。1stepあたりの計算速度が遅い。 |
| 公式サイト | [[The BUGS project>http://www.mrc-bsu.cam.ac.uk/software/bugs/the-bugs-project-winbugs/]] | [[OpenBUGS>http://www.openbugs.net/w/FrontPage]] | [[JAGS>http://mcmc-jags.sourceforge.net/]] | [[Stan>http://mc-stan.org/]] |
| インストール方法 | [[北大の久保先生のページ>http://hosho.ees.hokudai.ac.jp/~kubo/ce/HowtoInstallWinBUGS.html]] | | | [[2012年度統計数理研究所共同研究集会「データ解析環境Rの整備と利用」>http://prcs.ism.ac.jp/useRjp/hiki.cgi?2012%C7%AF%C5%D9+%A5%C7%A1%BC%A5%BF%B2%F2%C0%CF%B4%C4%B6%ADR%A4%CE%C0%B0%C8%F7%A4%C8%CD%F8%CD%D1]]の石倉さんのスライド、[[BUGSstan勉強会>http://atnd.org/events/43260]]、[[R stan導入公開版>http://www.slideshare.net/KojiKosugi/r-stan]] |
#contents
** 使うデータ
#データの生成
set.seed(1)
N <- 100
Nplot <- 4
Nsigma <- 2
y <- rnorm(N, 0, 5)
x1 <- y + rnorm(N, 0, 5)
x2 <- rnorm(N, 0, 1)
x3 <- as.numeric(cut(y + rnorm(N, 0, 5), Nplot, labels=1:Nplot))
** WinBUGS
*** データの定義
データは、データ名のベクトルを与えます。
data <- c("N", "Nplot", "Nsigma", "y", "x1", "x2", "x3")
*** 初期値の設定
WinBUGSにまつわるエラーの大半はここが原因です。間違いのないように指定しましょう。
inits <- list(list(alpha = 0, bx1 = 0, bx2 = 0, bx3 = rnorm(Nplot, 0, 1), sigma = rnorm(Nsigma, 50, 1)),
list(alpha = -5, bx1 = -5, bx2 = -5, bx3 = rnorm(Nplot, 0, 1), sigma = rnorm(Nsigma, 50, 1)),
list(alpha = 5, bx1 = 5, bx2 = 5, bx3 = rnorm(Nplot, 0, 1), sigma = rnorm(Nsigma, 50, 1))
)
*** 結果として得たいパラメータの定義
para <- c("alpha", "bx1", "bx2", "bx3", "sigma")
*** モデルの記述
#モデルのファイルを別に用意する場合は、model { }の部分だけ書いてファイルに保存して下さい。
modelFilename = "testmod.txt"
cat("
model {
for (i in 1:N) {
y[i] ~ dnorm(mu[i], tau[1])
mu[i] <- alpha + bx1*x1[i] + bx2*x2[i] + bx3[x3[i]]
}
alpha ~ dnorm(0.0, 1.0E-6)
bx1 ~ dnorm(0.0, 1.0E-6)
bx2 ~ dnorm(0.0, 1.0E-6)
for (i in 1:Nplot) {
bx3[i] ~ dnorm(0.0, tau[2])
}
for (i in 1:Nsigma) {
tau[i] <- pow(sigma[i], -2)
sigma[i] ~ dunif(0, 100)
}
}
", fill=TRUE, file=modelFilename)
*** 計算の実行
library(R2WinBUGS)
res <- bugs(data, inits, para, model=modelFilename,
n.burnin = 500, n.iter=6500, n.chains=3, n.thin=6,
bugs.directory="C:/WinBUGS14/", #デフォルトのフォルダから変更している場合のみ指定
working.directory=getwd(),
debug=TRUE
)
*** 結果の確認
print(res$summary, digits=3)
plot(res)
- 昔の[[R2WinBUGSに関するページ>研究ソフト/R2WinBUGS]]
** R2OpenBUGS
基本的に、WinBUGSの場合とほとんど一緒です。データの定義、初期値の定義、パラメータの指定は完全に一緒。
*** データの定義
データは、データ名のベクトルを与えます。
data <- c("N", "Nplot", "Nsigma", "y", "x1", "x2", "x3")
*** 初期値の設定
inits <- list(list(alpha = 0, bx1 = 0, bx2 = 0, bx3 = rnorm(Nplot, 0, 1), sigma = rnorm(Nsigma, 50, 1)),
list(alpha = -5, bx1 = -5, bx2 = -5, bx3 = rnorm(Nplot, 0, 1), sigma = rnorm(Nsigma, 50, 1)),
list(alpha = 5, bx1 = 5, bx2 = 5, bx3 = rnorm(Nplot, 0, 1), sigma = rnorm(Nsigma, 50, 1))
)
*** 結果として得たいパラメータの定義
para <- c("alpha", "bx1", "bx2", "bx3", "sigma")
*** モデルの記述
#モデルのファイルを別に用意する場合は、model { }の部分だけ書いてファイルに保存して下さい。
modelFilename = "testmod.txt"
cat("
model {
for (i in 1:N) {
y[i] ~ dnorm(mu[i], tau[1])
mu[i] <- alpha + bx1*x1[i] + bx2*x2[i] + bx3[x3[i]]
}
alpha ~ dnorm(0.0, 1.0E-6)
bx1 ~ dnorm(0.0, 1.0E-6)
bx2 ~ dnorm(0.0, 1.0E-6)
for (i in 1:Nplot) {
bx3[i] ~ dnorm(0.0, tau[2])
}
for (i in 1:Nsigma) {
tau[i] <- pow(sigma[i], -2)
sigma[i] ~ 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
)
** JAGS
前者2つと比べると、多少データの指定方法などが異なります。
*** データの定義
データは、リスト形式で準備して下さい。
list.data <- list(N=N, Nplot=Nplot, Nsigma=Nsigma, y=y, x1=x1, x2=x2, x3=x3)
*** 初期値
ここはWinBUGSなどと一緒です。
inits <- list(list(alpha = 0, bx1 = 0, bx2 = 0, bx3 = rnorm(Nplot, 0, 1), sigma = rnorm(Nsigma, 50, 1)),
list(alpha = -5, bx1 = -5, bx2 = -5, bx3 = rnorm(Nplot, 0, 1), sigma = rnorm(Nsigma, 50, 1)),
list(alpha = 5, bx1 = 5, bx2 = 5, bx3 = rnorm(Nplot, 0, 1), sigma = rnorm(Nsigma, 50, 1))
)
*** 結果として得たいパラメータの指定
ここもWinBUGSなどと一緒です。
para <- c("alpha", "bx1", "bx2", "bx3", "sigma")
*** モデルの記述
#モデルのファイルを別に用意する場合は、model { }の部分だけ書いてファイルに保存して下さい。
modelFilename = "testmod.txt"
cat("
model {
for (i in 1:N) {
y[i] ~ dnorm(mu[i], tau[1])
mu[i] <- alpha + bx1*x1[i] + bx2*x2[i] + bx3[x3[i]]
}
alpha ~ dnorm(0.0, 1.0E-6)
bx1 ~ dnorm(0.0, 1.0E-6)
bx2 ~ dnorm(0.0, 1.0E-6)
for (i in 1:Nplot) {
bx3[i] ~ dnorm(0.0, tau[2])
}
for (i in 1:Nsigma) {
tau[i] <- pow(sigma[i], -2)
sigma[i] ~ dunif(0, 100)
}
}
", fill=TRUE, file=modelFilename)
*** 計算の実行
#パッケージの読み込み
library(rjags)
#計算開始時間を記録
start.time <- Sys.time()
#初期化
m <- jags.model(file = modelFilename,
data = list.data,
# inits = inits,
#JAGSは初期値を与えなくても上手く計算してくれます。こける場合は指定してください。
n.chain = 3
)
update(m, 500)
x <- coda.samples(m, para,
thin = 6, n.iter = 6000
)
#終了時間を記録、表示
end.time <- Sys.time()
elapsed.time <- difftime(end.time, start.time, units='hours')
cat(paste(paste('Posterior computed in ', elapsed.time, sep=''), ' hours\n', sep=''))
*** 結果の確認
res <- data.frame(summary(x)$statistics)
ci <- data.frame(summary(x)$quantiles)
res$sig <- abs(sign(ci[, 1]) + sign(ci[, 5])) == 2
# Calculation of "R-hat" for convergence checking
rhat <- gelman.diag(x)[["psrf"]][, 1]
res$Rhat <- rhat
** Stan
StanはBUGS系と比べると、色々な点で違いがあります。
*** データの定義
list.data <- list(N=N, Nplot=Nplot, Nsigma=Nsigma, y=y, x1=x1, x2=x2, x3=x3)
*** 初期値の設定
inits <- list(list(alpha = 0, bx1 = 0, bx2 = 0, bx3 = rnorm(Nplot, 0, 1), sigma = rnorm(Nsigma, 50, 1)),
list(alpha = -5, bx1 = -5, bx2 = -5, bx3 = rnorm(Nplot, 0, 1), sigma = rnorm(Nsigma, 50, 1)),
list(alpha = 5, bx1 = 5, bx2 = 5, bx3 = rnorm(Nplot, 0, 1), sigma = rnorm(Nsigma, 50, 1))
)
*** 結果として得たいパラメータの定義
para <- c("alpha", "bx1", "bx2", "bx3", "sigma")
*** モデルの記述
modelfile <- '
data {
int<lower=0> N; //整数はint
int<lower=0> Nplot;
int<lower=0> Nsigma;
real y[N]; //実数はreal
real x1[N];
real x2[N];
int<lower=0> x3[N];
}
parameters {
real alpha;
real bx1;
real bx2;
real bx3[Nplot];
real<lower=0> sigma[Nsigma]; //分散パラメータなので、負にならないように範囲を指定
}
model {
for (i in 1:Nplot)
bx3[i] ~ normal(0.0, sigma[2]);
for (i in 1:N)
y[i] ~ normal(alpha + bx1*x1[i] + bx2*x2[i] + bx3[x3[i]], sigma[1]);
//事前分布を指定しない場合、自動的に一様分布が事前分布として与えられます
}
'
*** 計算の実行
library(rstan)
#BUGS系と引数の名称が微妙に異なります
fit <- stan(model_code = modelfile, data = list.data, par = para, inits=inits,
iter = 1000, chains = 3, thin=1, init=inits)
[[Hirokiさん>http://ito-hi.blog.so-net.ne.jp/2015-07-26]]によると、Stanのバージョン2.7.0から並列化が簡単にできるようになったとのことです。便利！
#並列化する場合
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
fit <- stan(model_code = modelfile, data = list.data, par = para, inits=inits,
iter = 1000, chains = 3, thin=1, init=inits)
*** 結果の確認
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"))
