Zerohours: 2月 2009

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

MHアルゴリズムによるN(0,1)からのサンプリングを行う。提案分布は一様分布とする。乱数の自己相関の状況について、特に採択率やthinningを変えた場合の影響を例示する。

１．Rコードの一例
thin <- 0 #サンプリング間隔（0のとき間引きなし）
thin1 <- thin+1
Ne <- 1000 # 有効サンプル数
N <- thin1*Ne # 総サンプル数
samp <- numeric(Ne)
count <- 0 # 採択数のカウント
x0 <- 4.5 # 初期値
a <- 4 # 採択率の調整
samp[1] <- x0
for (i in 1:N) {
cand <- runif(1, -a, a) # 提案分布（一様分布）から乱数生成
x1 <- x0 + cand # random walk
alpha <- min(1, dnorm(x1)/dnorm(x0)) # 採択確率の計算、提案分布は対称
u <- runif(1)
if (u < alpha){
x0 <- x1
count <- count + 1
}
if(i %% thin1 == 0){
j <- (i+thin1)/thin1
samp[j] <- x0
}
}
count/N #採択率
mean(samp)
sd(samp)
par(mfrow=c(2,1))
plot(samp,type="l")
hist(samp) # 密度関数を描くときは plot(density(samp))
par(mfrow=c(1,1))
par(ask=TRUE)
NN <- N/2
for(i in 1:NN){
j <- 2*i-1
plot(samp[j],samp[j+1],xlim=c(-5,5),ylim=c(-5,5),col="blue")
par(new=T)
}
par(new=F)
acf(samp) # 自己相関係数
# http://aoki2.si.gunma-u.ac.jp/R/acf.html には自己相関係数を計算する
# 関数acf2(x,k)のRコードが掲載されている。

２．採択率と調整パラメータの関係など（有効サンプル数=10000）