[R program]Rでマルコフ連鎖モンテカルロ法

RでMarkov chain Monte Carlo（MCMC、マルコフ連鎖モンテカルロ）を行うには大きく2つの方法があります。

1つ目はパッケージR2WinBUGSを使って、Rの後ろでWinBugsを走らせる方法。

もう1つはパッケージMCMCpackを使って直接R上で計算する方法。

ちなみにMCMCの利点は、混合効果モデルのように複雑なモデリングになった場合の複雑なパラメータ推定を安定にさせるというもの。

ロジスティック回帰で交互作用がいっぱいのときも安定するらしい。

それではパッケージの説明です。

• MCMCpack

#MCMCpack
library(MCMCpack)

y <- c(1,3,3,3,5)
x1 <- c(1,2,3,5,5)
x2 <- c(2,3,1,5,2)
x <- matrix(c(x1,x2),ncol=2,nrow=5)

result3.sim <- MCMCregress(y~x, data = parent.frame(), burnin = 1000, mcmc = 10000) 
plot(result3.sim)
raftery.diag(result3.sim)
summary(result3.sim)


#普通の線形回帰
summary(glm(y ~ x))


#MCMCpackのロジスティック回帰
data(birthwt)
posterior <- MCMClogit(low~age+as.factor(race)+smoke, data=birthwt)
plot(posterior)
summary(posterior)

#普通のロジスティック回帰
summary(glm(low~age+as.factor(race)+smoke, family=binomial, data=birthwt))

MCMCでのパラメータと、普通の最尤法の結果の比較です。

MCMCの中央値と最尤法の点推定値がほぼ同じですね。

plotで、パラメータの推移状況を見れます。

• R2WinBUGS

まずはここからWinBugsをインストールします。

私の場合はインストール以外の操作は必要ありませんでした。

Rでのコードは次のようになります（一部Rjpwiki参考）。

library(R2WinBUGS)

#↓"R2WinBUGS/model/schools.txt"という意味
model.file <- system.file(package="R2WinBUGS", "model", "schools.txt")	

#"R/R-210~1.0/library/R2WinBUGS/model/schools.txt"にあるschools.txtの中身
#↓WinBugsのコンソールに書かなきゃいけないモデルをテキストで保存している
#   model {
#       for (j in 1:J){
#         y[j] ~ dnorm (theta[j], tau.y[j])
#         theta[j] ~ dnorm (mu.theta, tau.theta)
#         tau.y[j] <- pow(sigma.y[j], -2)
#       }
#       mu.theta ~ dnorm (0.0, 1.0E-6)
#       tau.theta <- pow(sigma.theta, -2)
#       sigma.theta ~ dunif (0, 1000)
#     }  
file.show(model.file)


#平均値と標準偏差を指定しているけどどんなモデルなのだろ？
data(schools)
schools
J <- nrow(schools)
y <- schools$estimate
sigma.y <- schools$sd
data <- list ("J", "y", "sigma.y")

inits <- function(){
	list(theta=rnorm(J, 0, 100), mu.theta=rnorm(1, 0, 100),
	sigma.theta=runif(1, 0, 100))
}
parameters <- c("theta", "mu.theta", "sigma.theta")
schools.sim <- bugs(data, inits, parameters, model.file,
	n.chains=3, n.iter=5000, bugs.directory="c:/Program Files/WinBUGS14/")
print(schools.sim)
plot(schools.sim)


#Rjpwikiの例
#以前のワーキングスペースの保存
OldWd <- getwd()
#Bugsのモデルが保存してあるディレクトリを指定
setwd("C:/Documents and Settings/Issei/デスクトップ/")

#↓をmodel.txtという名前でワーキングディレクトリに保存（#は消す）
#model{
#for(i in 1:n){y[i]~dnorm(mu[i],tau)
#mu[i] <-alpha + beta1*x[i,1]+beta2*x[i,2]
#}
#alpha ~dnorm(0,0.0001)
#beta1 ~dnorm(0,0.0001)
#beta2 ~dnorm(0,0.0001)
#tau ~dgamma(0.001, 0.001)
#}

#データの作成
y <- c(1,3,3,3,5)
n <- NROW(y)
x1 <- c(1,2,3,5,5)
x2 <- c(2,3,1,5,2)
x <- matrix(c(x1,x2),ncol=2,nrow=5)
data <- list ("n", "y", "x") 

in1 <- list(alpha=0, beta1=0, beta2=0, tau=1)
in2 <- list(alpha=1, beta1=1, beta2=1, tau=1)
in3 <- list(alpha=2, beta1=2, beta2=2, tau=1)
inits <- list(in1,in2,in3)
parameters <- c("alpha", "beta1", "beta2", "tau")

result.sim <- bugs(data, inits, parameters, model.file="model.txt",
	n.chains = 3, n.iter = 1000, debug=FALSE, bugs.directory = "c:/Program Files/WinBUGS14/", 
	working.directory = NULL)

result.sim$sims.list$beta1
print(result.sim, digits=3)
plot(result.sim)
setwd(OldWd)

前半がRのexampleコードで、後半がRjpwikiのコードです。

後半のデータはMCMCpackのコードと同じものですので、同じような結果が出ています。

こちらのplotは中央値と80%区間が出ます。

でも結局WinBugsのコードを書かなきゃいけないので、これをやるくらいならMCMCpackを使うか直接Bugsを動かした方が楽ですかね。

ただ、混合効果ロジスティックモデルとか混合効果ポアソンモデルとかをMCMCでやりたい場合はMCMCpackでできるのかどうか調査中。

WinBugsでは出来るけど、コードの書き方を忘れたので今思い出している最中です。。。。

また、BPHOパッケージでは説明変数が全て（？）カテゴリの場合で、交互作用も多いときのMCMCロジスティックができたと思います。