uH35Nw  <a href="http://syglcuiopvus.com/">syglcuiopvus</a>, [url=http://czhixexwoxkp.com/]czhixexwoxkp[/url], [link=http://rwwbyfoznmva.com/]rwwbyfoznmva[/link], http://wkrbplcelrph.com/
* RとJAGS [#u10fa3b4]
- このページでは、RとJAGSを使ったMCMC、ギブス・サンプラーの計算方法を紹介します
- RのR2jagsパッケージを使います
** JAGSのインストール [#qf7f22e0]
- [[JAGSホームページ>http://mcmc-jags.sourceforge.net/]]にアクセス
- Windowsの場合
-- .exeファイルをダウンロードして実行する
- Mac OS X 10.4-10.6の場合
-- バイナリーファイルをダウンロードして実行する

** R2jagsパッケージ [#p78e3eb2]
- Rを起動し、R2jagsパッケージを呼び出す
 library(R2jags)

** 計算手順 [#f800ad35]
*** BUGSファイルの準備 [#nc648877]
- 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を使った計算[#u6f5698e]
- データ、初期値、パラメータ、モデルファイルを指定する。
 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