R2WinBUGSによる正規母集団パラメータのベイズ推定

BUGS（Bayesian Inference Using Gibbs Sampling）は以下のような方法で利用できる。
１）OpenBUGS やWinBUGS14をインストールし、単体で使用する。
２）RからライブラリR2WinBUGSを介してOpenBUGSやWinBUGS14を使用する。
３）RからライブラリBRugsを介してOpenBUGSを使用する。
ここでは２）の方法による一例を示す。
Xi（i=1～n）を正規分布（平均μ、分散σ²）に従う母集団からのi.i.d.サンプルとする。μ、σ、σ²等のベイズ推定を行う。
・参考文献
[1]大森裕浩「MCMCの基礎と統計科学への応用」、
http://www.e.u-tokyo.ac.jp/~omori/MCMC/mcmc-ism04.pdf、2004。
[2]中妻照雄「入門ベイズ統計学」、朝倉書店、p.145。

１．事前分布

以下では、ハイパーパラメータの値をμ₀=0、τ₀²=0.01、ν₀=36、λ₀＝1.36（文献[2]）と設定したときの例を示す。

２．modelファイル
次のmodelコードをファイルmodel.txtとして作業ディレクトリに保存する。
model{
mu ~ dnorm(0, 100) # 事前分布。BUGSではdnorm(mean, tau)=N(mean, σ²)、
# ここに、 tau = 1/σ² であることに注意。
tau ~ dgamma(18, 0.68) # 事前分布
var <- 1/tau # vaｒ ~ 逆gamma分布
sig <- sqrt(1/tau)
for(i in 1:n) {
X[i] ~ dnorm(mu, tau) # 尤度関数
}
}

３．Rコードの一例
# bugs関数ではmodel.file="model.txt"として２．で作成したファイルmodel.txtを指定する。
# データは参考文献[2]を使用。
# 結果の表示のためcodamenu()（ライブラリはcoda）を用いる。
# このRコードは適当なファイル名を付けて２．と同じ作業ディレクトリに保存する。
library(R2WinBUGS)
library(coda)
X <- c(4.590E-02, 4.360E-02, 2.070E-02, 8.670E-02, 1.678E-01, 1.748E-01)
n <- NROW(X)
data <- list("n","X")
ini1 <- list(mu=10, tau=0.4) # 初期値
param <- c("mu","tau","sig","var")
N <- 10000 # サンプル数
Nb <- 1000 # burn-in
out <- bugs(data,inits=ini1,param,model.file="model.txt", debug=T,
n.chains=2, n.thin=1, n.iter=N, n.burnin=Nb,
bugs.directory="C:/Program Files/WinBUGS14/", working.directory=NULL)
print(out,digits=4)
codamenu()

４．結果の表示（「・・・・・」は一部割愛を表す。）
> print(out,digits=4)
Inference for Bugs model at "model.txt", fit using WinBUGS,
2 chains, each with 10000 iterations (first 1000 discarded)
n.sims = 18000 iterations saved
mean sd 2.5%
mu 0.0567 0.0602 -0.0622
tau 29.7734 6.6214 18.2000
sig 0.1867 0.0213 0.1505
var 0.0353 0.0082 0.0227 ・・・・・
97.5% Rhat n.eff
mu 0.1760 1.0010 18000
tau 44.1400 1.0010 18000
sig 0.2344 1.0010 18000
var 0.0549 1.0010 18000
For each parameter, n.eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
・・・・・
> codamenu()
CODA startup menu
1: Read BUGS output files
2: Use an mcmc object
3: Quit
選択： 1
Enter CODA index file name 注） WinBUGS14はCODAファイルを自動生成する。
(or a blank line to exit)
1: codaIndex.txt
Enter CODA output file names, separated by return key
(leave a blank line when you have finished)
1: coda1.txt
2: coda2.txt
3:
Abstracting deviance ... 9000 valid values
Abstracting mu ... 9000 valid values
・・・・・
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
選択： 1

Save plots as a postscript file (y/N) ?
N
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