RとJAGS

  • このページでは、RとJAGSを使ったMCMC、ギブス・サンプラーの計算方法を紹介します
  • RのR2jagsパッケージを使います

JAGSのインストール

  • JAGSホームページにアクセス
  • Windowsの場合
    • .exeファイルをダウンロードして実行する
  • Mac OS X 10.4-10.6の場合
    • バイナリーファイルをダウンロードして実行する

R2jagsパッケージ

  • Rを起動し、R2jagsパッケージを呼び出す
    library(R2jags)

計算手順

BUGSファイルの準備

  • BUGSファイル(例えば、以下のサンプルファイル)をテキストファイルとして保存する(例:model1.txt)
    model{
    for(i in 1:n){
    y[i] ~ dnorm(mu[i], tau)
    mu[i] <- beta0 + beta1 * X[i,1] + beta2 * X[i,2] + beta3 * X[i,3] + beta4 * X[i,4] + beta5 * X[i,5]
    }
    beta0 ~ dnorm(0,1.0E-6)
    beta1 ~ dnorm(0,1.0E-6)
    beta2 ~ dnorm(0,1.0E-6)
    beta3 ~ dnorm(0,1.0E-6)
    beta4 ~ dnorm(0,1.0E-6)
    beta5 ~ dnorm(0,1.0E-6)
    tau ~ dgamma(0.01, 0.01)
    }
  • BUGSファイルを適当な場所に置く
    • R2jagsパッケージをインストールしている場合、Rインストール先フォルダの下にある、R2jagsフォルダ下のmodelフォルダに置く
    • Windowsの場合
      • 「C:\Program Files\R\R-2.x.x\library\R2jags\model」にファイルを置く
      • R画面の「ファイル」→「ディレクトリの変更」で上記フォルダを作業領域として指定する
    • Macの場合
      • 「Macintosh HD」→「ライブラリ」→「Frameworks」→「R.frameworks」→「Resources」→「library」→「R2jags」→「model」にファイルを置く

R2jagsを使った計算

  • データ、初期値、パラメータ、モデルファイルを指定する。
    data(swiss)
    y<-swiss[,1]
    k<-ncol(swiss)
    X<-as.matrix(swiss[,2:k])
    n<-nrow(X)
    data <- list("n", "y", "X")
    in1 <- list(beta0=0, beta1=0, beta2=0, beta3=0, beta4=0, beta5=0, tau=1)
    in2 <- list(beta0=0, beta1=1, beta2=1, beta3=1, beta4=1, beta5=1, tau=1)
    in3 <- list(beta0=0, beta1=2, beta2=2, beta3=2, beta4=2, beta5=2, tau=1)
    inits <- list(in1,in2,in3)
    parameters <- c("beta0", "beta1","beta2","beta3","beta4","beta5","tau")
    model.file <- system.file(package="R2jags", "model", "model1.txt")
  • jags()関数で計算する
    swiss.jags <- jags(data=data, inits=inits, parameters, n.iter=10000, 
    n.burnin=1000, n.chains=3, model.file=model.file)
  • 計算結果を表示する
    print(swiss.jags, digits=3)
    plot(swiss.jags)
    traceplot(swiss.jags)
  • MCMCが収束しなかった場合、update()関数で収束させる
    swiss.fit <- update(swiss.jags)
    print(swiss.fit, digits=3)

トップ   編集 凍結解除 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2014-07-31 (木) 00:20:45 (1596d)