2014-03-28

openBUGS

BUGS ベイズ MCMC openBUGS Gibbs sampling

Bayesian inference using Gibbs sampling:BUGS は、ベイズ推定を計算機統計学的にやってしまいましょうという手法(Wiki記事)
計算機統計学的に…という部分は、「難しい関数記述(の省略できるところ)省略して、MCMC/Gibbs samplingで代用する」という風に読み替えてもよい
Rと連動する使い方として、DynacomさんのサイトにBUGSを使う場合(こちら)とStanを使う場合(こちら)との紹介があります
簡単に言うと、"30分環境設定！"で動かせたのはopenBUGSの方でrstanはエラーを吐きました
以下、openBUGS(rstanは後続の記事で)
基本的には、こちらをべたなぞりです(ただし、リンク先・バージョンアップなどがあり、少し、修正が必要でした)。べたなぞりで出来る、というのは本当にありがたいです。
openBUGSをこちらからとってきてインストール
記事に言われるがままにパッケージをインストール

install.packages(c("BRugs", "coda", "R2WinBUGS","Rcpp"))

library(R2WinBUGS) 
library(BRugs)

ごく簡単な例をRのlm()の結果と比較しながら実行してみます

# 適当にデータを作る
set.seed(314) 
x <- rnorm(100, mean=3, sd=1) 
y <- 3.00 + 0.25*x + rnorm(100, mean=0, sd=0.10)
# そのデータを見てみて、lm()で回帰して回帰直線を引きます
plot(y ~ x) 
lm1 <- lm(y ~ x) # Ordinal Least Square Method 
abline(lm1, col="blue", lwd=2) 
summary(lm1) 

# BUGSをするにあたって、モデルを扱うために変数・オブジェクトの整備をする
N <- length(x) 
data1 <- list(y=y, x=x, N=N) 
 
in1 <- list(beta=c(0, 0), S=1) 
inits <- list(in1) 
 
param <- c("beta", "S") 
# テキストファイルでBUGS用のモデルを記載(後で説明
model1 <- "lm1.txt" 
# 実行
# model1に指定したファイルのパスはもちろんRの実行ディレクトリに合わせましょう
# bugs.directory="C:/Program Files(x86)/OpenBUGS/OpenBUGS323"という指定はopenBUGSを動かすための指定なので、インストールした先のディレクトリの名前を確認しておきましょう(参考ＰＤＦとはインストールバージョンが違ったので、書き換えてあります
bugs1 <- bugs(data1, inits, param, model.file=model1, 
n.chains=1, n.iter=205000, n.burnin=5000, n.thin=2, debug=TRUE, 
program="OpenBUGS", 
bugs.directory="C:/Program Files(x86)/OpenBUGS/OpenBUGS323") 
# 終了(大きいデータ・複雑なモデルだと時間がかかりますが、この例は簡単)
print(bugs1)
post1 <- bugs1$sims.matrix
apply(post1,2,mean) 
apply(post1,2,median)
apply(post1,2,sd) 
plot(bugs1)

par(mfrow=c(1,2)) 
 
plot(post1[1:5000,1], type="l", col=4) 
plot(post1[1:5000,2], type="l", col=4) 


hist(post1[,1],freq=FALSE,xlab="beta[1]",main="Posterior 
distribution of beta[1]") 
lines(density(post1[,1], bw = 0.002),col=4) 
 
hist(post1[,2],freq=FALSE,xlab="beta[2]",main="Posterior 
distribution of beta[2]") 
lines(density(post1[,2], bw = 0.001),col=4)

BUGSのモデル用のファイル"lm1.txt"
- $y ~ N(\mu,S)$ で、yは色々な値 $\mu$ を真値としてsd=Sのばらつきのある値(ただし、Sはガンマ関数に基づいて決める)、と指定
- 真値 $\mu$ は固定項とx依存項の線形和で決まるものとし
- その係数は0周辺の正規乱数、と書いてある

model{ 
 
 for(i in 1:N) { 
 
 y[i] ~ dnorm(mu[i], S) 
 mu[i] <- beta[1] + beta[2] * x[i] 
 
 } 
 
 beta[1] ~ dnorm(0, 0.0001) 
 beta[2] ~ dnorm(0, 0.0001) 
 S ~ dgamma(0.0001, 0.0001) 
 
}

これをさらに大きく複雑にしたのが以下…なかなか終わらない…

install.packages(c("lme4","geepack"))
library(lme4) 
library(geepack) 
 
data(ohio)
fm1 <- glmer(ohio$resp ~ ohio$age + ohio$smoke + (1|ohio$id),  family=binomial) 
summary(fm1) 


N <- length(ohio$resp) 
ID <- as.numeric(factor(ohio$id)) 
NS <- length(unique(ID)) 
 
data1 <- list(resp=ohio$resp, age=ohio$age, smoke=ohio$smoke, ID=ID, N=N, NS=NS) 
 
in1 <- list(beta=c(0, 0, 0), gamma=numeric(NS), Tau=1) 
inits <- list(in1) 
 
param <- c("beta", "Tau") 
 
model1 <- "lm2.txt" 
 
bugs1 <- bugs(data1, inits, param, model.file=model1, 
 n.chains=1, n.iter=55000, n.burnin=5000, n.thin=1, 
 debug=TRUE, program="OpenBUGS", 
 bugs.directory="C:/Program Files (x86)/OpenBUGS/OpenBUGS323") 
post1 <- bugs1$sims.matrix 
 
print(bugs1)

モデルファイル"lm2.txt"

model{ 
 
 for(i in 1:N) { 
 
 resp[i] ~ dbin(p[i], 1) 
 logit(p[i]) <- gamma[ID[i]] + beta[1] + beta[2] * age[i] + beta 
[3] * smoke[i] 
 
 } 
 
 beta[1] ~ dnorm(0, 0.0001) 
 beta[2] ~ dnorm(0, 0.0001) 
 beta[3] ~ dnorm(0, 0.0001)
 for(j in 1:NS) { gamma[j] ~ dnorm(0, Tau) } 
 
 Tau ~ dgamma(0.0001, 0.0001) 
 
}

結果

> print(bugs1) 
Inference for Bugs model at "lm2.txt", fit using OpenBUGS,
 1 chains, each with 55000 iterations (first 5000 discarded)
 n.sims = 50000 iterations saved
           mean   sd   2.5%    25%    50%    75%  97.5%
beta[1]    -3.1  0.2   -3.6   -3.3   -3.1   -3.0   -2.7
beta[2]    -0.2  0.1   -0.3   -0.2   -0.2   -0.1    0.0
beta[3]     0.4  0.3   -0.1    0.2    0.4    0.6    0.9
Tau         0.2  0.0    0.1    0.2    0.2    0.2    0.3
deviance 1166.2 31.3 1107.2 1144.8 1165.4 1186.8 1229.4

DIC info (using the rule, pD = Dbar-Dhat)
pD = 248.0 and DIC = 1414.0
DIC is an estimate of expected predictive error (lower deviance is better).
>