Xi（i=1～n）を正規分布（平均μ、分散σ2）に従う母集団からのi.i.d.サンプルとする。μ、σ、σ2等のベイズ推定を行う。なお、事前分布やハイパーパラメータの値には以下と同じものを用いる。
・R2WinBUGSによる正規母集団パラメータのベイズ推定
（参考文献）豊田秀樹編著「マルコフ連鎖モンテカルロ法」、朝倉書店、2008年。
１．パッケージBRugsのインストール
Rを起動して「パッケージ」－「パッケージのインストール」をクリックするとCRANのミラーサイトを選択する画面が表示される。適当なサイトを選択して「OK」を押すと、パッケージを選択する画面が表示される。
「BRugs」を選択して「OK」を押すとダウンロードが始まる。パッケージは約5.5MBのzipファイルで、一時ファイルフォルダにダウンロードされる。
ダウンロードされたパッケージをRに読み込むため、「パッケージ」－「パッケージの読み込み」をクリックする。
パッケージから「BRugs」を選択し「OK」を押すと読み込みが始まる。読み込み終了後、RのフォルダlibraryにフォルダBRugsがあることを確認すればインストールは完了する。
BUGSではモデルファイル、データファイル及び初期値ファイルを用いて事後分布からサンプリングを行う。事前分布や尤度関数等の情報はモデルファイルに、所与のデータ情報はデータファイルに、計算に用いる初期値の情報は初期値ファイルに収められている。使用する各ファイルを適当なエディタを用いて以下のように作成し、ワーキングディレクトリに保存する。
１）model.txt
model{
mu ~ dnorm(0, 100) # 事前分布
tau ~ dgamma(18, 0.68) # 事前分布
sig2 <- 1/tau
sig <- sqrt(1/tau)
for(i in 1:n) {
X[i] ~ dnorm(mu, tau) # 尤度関数（サンプルX[i]の母集団分布）
}
}
２）data.txt
list(n=6, X=c(4.59000E-02, 4.36000E-02, 2.07000E-02, 8.67000E-02, 1.67800E-01, 1.74800E-01))
３）inits.txt
list(mu=0, tau=1)
なお、初期値ファイルは連鎖の数だけ必要である。特に初期値の設定に条件がない場合は、連鎖の数に自動対応する初期値の生成関数modelGenInits()を用いると便利である。
３．BRugsの関数を用いたRコードの一例
以下のコードを適当なファイル名を付けてワーキングディレクトリに保存し、Rから実行する。
library(BRugs)
Ne <- 10000 # 有効サンプル数
Nb <- 1000 # burn-in
N <- Nb+Ne # 総繰返し数
modelCheck("model.txt")
modelData("data.txt")
modelCompile(numChains=1)
modelInits("inits.txt") # 初期値を自動生成する場合はmodelGenInits()を使用する。
samplesSet(c("mu","tau","sig","sig2")) # パラメータの指定。
modelUpdate(N)
samplesStats("*") # すべての変数の統計量を表示する。
samplesHistory("mu", beg=Nb)
samplesDensity("mu", beg=Nb)
samplesCoda("mu",stem="mu") # パラメータmuのCodaファイルをワーキングディレクトリに生成する。
CD <- read.openbugs(stem="mu") # パラメータmuのMCMCオブジェクトをCDとする。
codamenu()
４．結果の表示（「・・・・・」は一部割愛を表す。）
> modelCheck("model.txt")
model is syntactically correct
> modelData("data.txt")
data loaded
> modelCompile(numChains=1)
model compiled
> modelInits("inits1.txt")
Initializing chain 1: model is initialized
> samplesSet(c("mu","tau","sig","sig2"))
monitor set for variable 'mu'
monitor set for variable 'tau'
monitor set for variable 'sig'
monitor set for variable 'sig2'
> modelUpdate(N)
11000 updates took 0 s
> samplesStats("*")
mean sd MC_error val2.5pc median val97.5pc start sample
mu 0.05698 0.060660 5.836e-04 -0.06118 0.05783 0.17610 1 11000
sig 0.18650 0.020810 1.995e-04 0.15090 0.18470 0.23230 1 11000
sig2 0.03522 0.008037 7.804e-05 0.02277 0.03412 0.05396 1 11000
tau 29.80000 6.500000 5.965e-02 18.54000 29.31000 43.92000 1 11000
> samplesHistory("mu", beg=Nb)
CODA files written
>CD <- read.openbugs(stem="mu")
・・・・・
> codamenu()
CODA startup menu
1: Read BUGS output files
2: Use an mcmc object
3: Quit
Enter an item from the menu, or 0 to exit
選択： 2
Enter name of saved object (or type "exit" to quit)
1:CD
Checking effective sample size ...OK
CODA Main Menu
1: Output Analysis
2: Diagnostics
3: List/Change Options
4: Quit
選択： 1
CODA Output Analysis menu
1: Plots
2: Statistics
3: List/Change Options
4: Return to Main Menu
選択： 4
CODA Main Menu
1: Output Analysis
2: Diagnostics
3: List/Change Options
4: Quit
選択： 4
Are you sure you want to quit (y/N)?
:y
投稿
