☆ WBUGS(Mac,Win向け、R2WinBUGSサポート関数)

R2WinBUGSの使用が面倒で、利用を簡易化するため自分好みのサポート関数を作りました。
おもにMacでWine経由での使用を想定していますが、Windows7でも動作を確認できています
(cf. MacへのWinBUGSインストール方法

(2015.12.11 追記)OSX10.11 El Capitanに対応
(2014.05.23 追記)burninやthinの自動生成法を変更、dcloneに一部対応
(2013.01.08 追記)推定する変数がベクトルや行列、配列のとき、変数名と各次元の番地の列を追加するように変更
(2012.12.26 追記)mcmc.plotで、パラメータを1つずつ一枚の図にまとめるように変更
(2012.12.14 追記)mcmc.plot, mcmc.est, mcmc.zeroをより柔軟に改良、パラメータが、a[1], a[2]...のようなベクトルの場合にも対応するようにした。
(2012.12.06 追記)初期値の返し方にエラーを発見、修正…
ifelseが何故か、パラメータベクトルの最初の値しか返さなかったことが判明
→ if () else に書き直したら解決(ifelseは条件式と同一長分の結果しか返さないようだ)。
(2012.10.23 追記)modelnameがある場合だけ推定結果のファイルを生成するように変更
また、推定結果は 3 桁で出力するように変更
今更ながら、旧名は既存のものがあるのが分かって、当方のはWBUGSと改称
(2012.9.24 追記)パラメータのFやNAが全然無い場合にエラーになることに気づき修正


まず、このページの末尾の関数定義をコピペして、WBUGS.R、とでもファイル名をつけて作業ディレクトリに保存して使います。すると、このコードでwbugsを呼び出すことができます。
source("WBUGS.R") 


############ MCMC推定をするメインの関数 #############
wbugs(data, parameter.set, model)   # 括弧内は最低限必要な変数


## 通常のMCMC計算に加え、計算結果の出力と計算時間の表示をします。
#### データの形式
## data: 

Data <- list(X=X, Y=Y) # 通常のbugs()で用いるものと同様のデータリストを用意する。

#### パラメータの初期値、推定すべきパラメータの指定
## parameter.set: 

Parms <- list(a=0, b=0, tau=1, F, sigma=NA)  # パラメータリストはこういう形式で

# 通常のinitsとparametersに相当するものを統合している。
# 初期値を与えるパラメータは、=0 や =1のように数値を与える。
# 初期値を指定しないパラメータ(tauの逆数の平方根であるsigmaなど、推定値がほしいだけのもの)は =NA としておく。
# 初期値は与えるけれど、推定値は不要のパラメータ(分散の逆数 tauなど)は、次の項に F を入れておくと推定値を返さない。
#(この F の項は入れようが入れまいが、計算実行に害はありません。お好みで。)
# これで、同じパラメータを複数回記述したり、パラメータ名をいちいち "" で囲う手間から解放されます。

#### モデル式の記述
## model: 

Model <- function() { モデル式 }
のようにRコードとして記述してください。

注: WinBUGSの I() 関数を用いる場合、R上でエラーになってしまうため、ダミーコードを使用し、 %_% I() のようにすると動作する(cf. ?write.model)。

#### 他の項(指定しなくとも適当な値がデフォルトで割り当てられています):
# modelname: 

# "モデル名"の項。無くてもOKです、入れると計算結果のcsvファイルを出力します。
# 出力されるファイル名は「"モデル名"+MCMC回数など+DIC」となります。

#### その他の項はbugs()同様にユーザが指定可能。

n.chain, n.iter, n.burnin, n.thin, debug
などは適当なデフォルト値を設定しているが、自前で記述することも可能
デフォルトではMCMCサンプル数が3000になるように設定





############ 実行例 ###########################
#### 計算に必要なデータ、パラメータ、モデ
X <- c(1:100)
Y <- rnorm(100, mean=(3 + 2*X), sd=1)
Data <- list(X=X, Y=Y)
Parms <- list(a=0, b=0, tau=1, F, sigma=NA)
# 推定値が不要なパラメータ → 次の項に F を入れる(入れずとも推定には何も影響はしない)
# 初期値を設定する必要がなく、推定値のみが必要なパラメータ → =NAとする

Model <- function() {
   a ~ dnorm(0, 1.0E-6)
   b ~ dnorm(0, 1.0E-6)
   tau ~ dgamma(1.0E-2, 1.0E-2)
      for (i in 1:100) {
         Y[i] ~ dnorm(mean[i], tau)
         mean[i] <- a + b*X[i] }
   sigma <- 1/sqrt(tau) }

#### 下記のようにして推定計算を実行
## モデル名を"testmodel"のように付けると、MCMC回数やDICをファイル名に付けて.csvに出力します
source("WBUGS.R") 
mcmc <- wbugs(Data, Parms, Model, "testmodel", n.iter=11000)

#(モデル名を付けなければcsvファイルを出力しない)




#### 推定結果を見る(下記のように打つだけ)
mcmc

#### MCMC プロセス図示の実行(a,b,sigmaについての例)
mcmc.plot(mcmc, c("a","b","sigma"))

#### 事後分布をまとめて抽出(各MCMCの値、a,b,sigmaについての例
# 単に値を抽出。実行するとコンソールがあふれるので注意
mcmc.est(mcmc, c("a","b","sigma"))

# 二次元以上のパラメータには対応せず。その場合は地道に1つずつ取り出し↓
mcmc$sims.list$a    # aを取り出す場合(applyを使って料理しますか…)
mcmc$mean$a    # aの平均値を取り出す場合
mcmc$median$a    # aの中央値を取り出す場合(summaryの"50%"の値がmedian)


#### パラメータ事後分布間の散布図を図示
pairs(mcmc.est(mcmc, c("a","b","sigma")), panel=panel.smooth)

#### 事後分布の確率密度曲線("a"についての例:自分でグラフの体裁をカスタマイズしたい場合、mcmc.plot(2つ上の)よりもこちらが便利でしょう)
plot(density(mcmc.est(mcmc, "a")), main="", xlab="a", lwd=3, xlim=c(2,4), ylim=c(0,2))

#### 事後分布がある値以上または以下になる確率
# WinBUGSのstep()関数の後出しバージョン
#(デフォルトは 0、a,b,sigmaについての例
mcmc.zero(mcmc, c("a","b","sigma"))




############# WBUGS.R(ここから)##############
require(R2WinBUGS)
require(coda)
suppressWarnings(require(dclone))

wbugs <- function(data, # 通常のdata形式(list)
   parameter.set, # list:上記実行例のParmsを参照
   model, # model <- function() {} :上記実行例参照
   modelname=NULL, # 明記すると保存ファイル名の冒頭にモデル名を追加できる
   n.chains=3, n.iter=11000, # デフォルトの値
   n.burnin=1000, n.thin=10,
   debug=F, ...) {

  #print(paste("Start", Sys.time()))
   modelpath <- file.path(tempdir(), paste("model", ".bug", sep=""))
   write.model(model, modelpath) # Rで記述したモデルからファイルを生成

# デフォルトのMCMCサンプル数が3000になるように設定
   n.burnin0 <- as.numeric(substr(n.iter,2,nchar(n.iter)))
   if (is.null(n.burnin) & n.burnin0 == 0) n.burnin <- n.iter/10 else n.burnin <- n.burnin0
   if (is.null(n.thin)) n.thin <- round(as.numeric(substr(n.iter,1,1))*10^(nchar(n.iter)-1) / 1000) 
   if (is.null(n.thin) & n.burnin0 == 0) n.thin <- round((n.iter - n.burnin) / 1000)
   if (is.null(n.thin) & n.thin == 0) n.thin <- 1

 n.clones <- dclone::nclones.list(data) # Data clonig (Lele et al. 2007) に対応するコード
    data <- lapply(data, function(z) {
        attr(z, "n.clones") <- NULL
        z })

   false <- which(names(parameter.set)=="") # F の位置
   if (length(false)==0) { # F がない場合
      na <- which(is.na(parameter.set)) # NAの位置
      na <- if (length(na)==0) { length(parameter.set) + 1 } else { na } # NAの有無を判断
      inits <- lapply(1:n.chains, function(i) parameter.set[-na]) # 初期値のリストをn.chain数へ増やす
      parms.est <- names(parameter.set) # 推定したいパラメータ名
            } else { # F がある場合
      parameter.set <- parameter.set[-false] # F を除く(Fの位置はわかったのでもはや不要)
      except <- false - c(1:length(false)) # F:推定値不要のパラメータの位置(Fの一個前、累積和)
      na <- which(is.na(parameter.set)) # NAの位置
      na <- if (length(na)==0) { length(parameter.set) + 1 } else { na } # NAの有無を判断
      inits <- lapply(1:n.chains, function(i) parameter.set[-na]) # 初期値のリストをn.chain数へ増やす
      parms.est <- names(parameter.set)[-except] # 推定したいパラメータ名
      }

   if (Sys.info()["sysname"]=="Darwin") { # プラットフォームがMac OSXの場合
      if(as.numeric(substr(Sys.info()["release"],1,4)) >= 15.2) path.name <- "/usr/local/bin/" else path.name <- "/opt/local/bin/" # OSXのversionによってパスを変える
      time <- system.time(
      mcmc <- bugs(
            data=data, inits=inits, parameters=parms.est, model.file=modelpath,
            n.chains=n.chains, n.iter=n.iter, n.burnin=n.burnin, n.thin=n.thin,
            working.directory=NULL, clearWD=T, debug=debug, useWINE=T, newWINE=T,
            WINE=paste(path.name, "wine", sep=""), WINEPATH=paste(path.name, "winepath", sep=""),
            ...) )
      print(time)
            } else { # プラットフォームがMac以外の場合(Windowsを想定)
   time <- system.time(
      mcmc <- bugs(
            data=data, inits=inits, parameters=parms.est, model.file=modelpath,
            n.chains=n.chains, n.iter=n.iter, n.thin=n.thin, n.burnin=n.burnin,
            bugs.directory="c:/Program Files/WinBUGS14/",
            working.directory=NULL, clearWD=T, debug=debug, ...) )
            }
   ### 変数名の整理(推定する変数がベクトルや行列、配列になっている時の対処)
   results <- mcmc$summary
   est0 <- unlist(strsplit(rownames(results), "\\]")) # 最後の ] を削除
   est1 <- sapply(1:nrow(results), function(i) unlist(strsplit(est0[i], "\\["))[1]) # 変数名を抜き出し
        if (max(table(est1)) > 1) { # そういう変数が無ければ以降はやらない
     est2 <- sapply(1:nrow(results), function(i) unlist(strsplit(est0[i], "\\["))[2]) # 番地の抜き出し
     est.max <- est2[sapply(1:length(unique(est1)), function(i) min(which(est1==unique(est1)[i])))]
     est.max <- max(sapply(1:length(est.max), function(i) length(strsplit(est.max, ",")[[i]]))) # 最大次元数
     est2 <- sapply(1:nrow(results), function(i) strsplit(est2[i], ",")) # 番地を分解
     est2 <- data.frame(t(sapply(1:length(est2), # 番地が少ない項にNAを追加した上で結合
           function(i) c(as.numeric(est2[[i]]), rep(NA, (est.max - length(est2[[i]])))))))
     colnames(est2) <- 1:ncol(est2) # 次元数をナンバリング
     if (dim(est2)[1]==1) est2 <- unlist(est2) # ベクトルの場合、一次元に戻す
     results <- data.frame(estimates=est1, dim=est2, results) } # 変数名と次元数の列を付加
                  if (!is.null(modelname))
                  write.csv(results, paste(modelname,"Iter",n.iter,"Burn",n.burnin,"Thin",n.thin,
                  "pD",mcmc$pD,"DIC",mcmc$DIC,".csv",sep=""))
  #   print(paste("End", Sys.time()))
  print(mcmc, digit=3)
   }



## MCMCプロセスと事後分布の図示 (パラメータが、a[1], a[2]...のようなベクトルの場合にも対応)

mcmc.plot <- function(mcmc, parameters) {
 parm.name <- numeric(0)
 for (i in 1:length(parameters)) {
  estimates <- mcmc$sims.list[[parameters[i]]]
  parm.name <- c(parm.name, if (length(estimates) == mcmc$n.sims) { parameters[i]
  } else { paste(parameters[i], "[", c(1:ncol(estimates)), "]", sep="") } )
    }
   for (j in 1:length(parm.name)) {
   dev.set(1) 
   plot(mcmc.list(lapply(1:mcmc$n.chains,
   function(k) as.mcmc(mcmc$sims.array[, k, parm.name[j]]))))
   mtext(parm.name[j], side=3, cex=2)
   }}



## 事後分布の推定値抽出(全chainの値をプール)

mcmc.est <- function(mcmc, parameters) {
est <- numeric(0)
for (i in 1:length(parameters)) {
estimates <- mcmc$sims.list[[parameters[i]]]
if (length(estimates) == mcmc$n.sims) {
estimates <- matrix(estimates)
colnames(estimates) <- parameters[i]
} else {
colnames(estimates) <- paste(parameters[i], "[", c(1:ncol(estimates)), "]", sep="") }
est <- cbind(est, estimates) }
return(est) }



## 事後分布がある値以上または以下になる確率(デフォルトは 0)

mcmc.zero <- function(mcmc, parameters, value=0) {
est <- numeric(0)
for (i in 1:length(parameters)) {
estimates <- mcmc$sims.list[[parameters[i]]]
if (length(estimates) == mcmc$n.sims) {
estimates <- matrix(estimates)
colnames(estimates) <- parameters[i]
} else {
colnames(estimates) <- paste(parameters[i], "[", c(1:ncol(estimates)), "]", sep="") }
est <- cbind(est, estimates) }
  Ps <- apply(est, 2, function(p) length(p[p > value])/nrow(est))
  ifelse(Ps > 0.5, 1 - Ps, Ps) }


## 2つのパラメータ推定値(est1、est2)が異なる確率
# 差のパラメータ(d.est <- est1 - est2)を推定するBUGSコードを書いておき、d.estが0以上である確率を求めればよい




## 変数名と次元数を別々に取り出して付加する関数(wbugs内で定義しているが、二次利用するために、外部関数化)
est.index <- function(mcmc) {
   results <- mcmc$summary
   est0 <- unlist(strsplit(rownames(results), "\\]")) # 最後の ] を削除
   est1 <- sapply(1:nrow(results), function(i) unlist(strsplit(est0[i], "\\["))[1]) # 変数名を抜き出し
   est2 <- sapply(1:nrow(results), function(i) unlist(strsplit(est0[i], "\\["))[2]) # 番地を抜き出し
   est.max <- est2[sapply(1:length(unique(est1)), function(i) min(which(est1==unique(est1)[i])))] # 各変数の一番目
   est.max <- max(sapply(1:length(est.max), function(i) length(strsplit(est.max, ",")[[i]]))) # 最大次元数
   est2 <- sapply(1:nrow(results), function(i) strsplit(est2[i], ",")) # 番地を分解
   est2 <- data.frame(t(sapply(1:length(est2), # 番地が少ない項にNAを追加した上で結合
           function(i) c(as.numeric(est2[[i]]), rep(NA, (est.max - length(est2[[i]])))))))
   colnames(est2) <- 1:ncol(est2) # 次元数をナンバリング
         if (dim(est2)[1]==1) est2 <- unlist(est2) # ベクトルの場合、一次元に戻す
     data.frame(estimates=est1, dim=est2, results)  } # 変数名と次元数の列を付加


############# WBUGS.R(ここまで)##############

0 件のコメント:

コメントを投稿