「研究ソフト/MCMCサンプラー」の編集履歴（バックアップ）一覧はこちら

「研究ソフト/MCMCサンプラー」(2017/06/22 (木) 20:55:24) の最新版変更点

追加された行は緑色になります。

削除された行は赤色になります。

* 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"))

このページは、[[こちら>https://hayatoiijima.jimdo.com/%E7%B5%B1%E8%A8%88%E8%A7%A3%E6%9E%90/mcmc%E3%82%B5%E3%83%B3%E3%83%97%E3%83%A9%E3%83%BC/]]に移転しました。