MCMCの収束判定
MCMCによる正規分布のベイズ推定
- マルコフ連鎖モンテカルロ法（メトロポリス法）による正規母集団のサンプリング（第３回演習の内容と同じ）
# 正規母集団の平均と分散⇒事前分布超パラメータの初期値
mu <- 0; sig <- 1
# 正規母集団から抽出された標本数
n <- 30
# 正規母集団から抽出した標本データ
y <- rnorm(n, mu, sig)
# 標本平均
yvar <- mean(y)
# 事前分布の超パラメータ
mu0 <- 0 # 平均
sig0 <- 4 # 標準偏差
nu0 <- 1 # 精度の事前ガンマ分布超パラメータ
S0 <- 4 # 同上
# 繰り返し計算回数
iter <- 10000
# MCMCのチェーン数
n.chains <- 3
# 事後分布
posterior <- array(dim=c(n.chains,2,iter))
# 受容される分布
accepted <- array(dim=c(n.chains, iter-1))
# メトロポリス法
for(c in 1:n.chains){
# 事後分布
theta.post <- array(dim=c(2, iter))
# 初期値
mu.init <- rnorm(1,mu0,sqrt(sig0))
#tau.init <- rgamma(1,nu0/2,S0/2)
#sig.init <- sqrt(1/tau.init)
sig.init <- runif(1, 0, 2)
theta.post[1,1] <- mu.init
theta.post[2,1] <- sig.init
for(t in 2:iter){
# 需要棄却法
# 提案分布
h.theta <- c(theta.post[1,t-1]+rnorm(1,0,0.1),
theta.post[2,t-1]+rnorm(1,0,0.1))
# 事後分布(対数)=対数尤度+事前平均(対数)+事前分散(対数)
p.post <- sum(dnorm(y, h.theta[1], h.theta[2], log=T))
+ dnorm(h.theta[1], mu0, sig0, log=T)
+ dunif(h.theta[2], 0, 2, log=T)
p.prev <- sum(dnorm(y, theta.post[1,t-1], theta.post[2,t-1], log=T))
+ dnorm(theta.post[1,t-1], mu0, sig0, log=T)
+ dunif(theta.post[2,t-1], 0, 2, log=T)
lr <- p.post - p.prev
r <- exp(lr)
accept <- rbinom(1,1,prob=min(r,1))
accepted[c,t-1] <- accept
if(accept==1){
theta.post[,t] <- h.theta
}
else {
theta.post[,t] <- c(theta.post[1,t-1], theta.post[2,t-1])
}
}
posterior[c,,] <- theta.post
}
事後分布の可視化
library(coda)
# MCMCへの変換
data <- as.mcmc(posterior[1,1,])
summary(data)
##
## Iterations = 1:10000
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 10000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## 0.221045 0.278778 0.002788 0.018655
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## -0.1648 0.0737 0.2062 0.3345 0.6202
plot(data)
autocorr.plot(data)
rejectionRate(data)
## var1
## 0.3110311
Gelman&Rubinの診断方法
- Within-chain variance (W)とBetween-chain variance (B)を計算したのち、Potential Scale Reduction Factor (PSRF)を計算する。chain数は２つ以上あることが前提。
thita <- rbind(posterior[1,1,], posterior[2,1,], posterior[3,1,])
thita <- as.matrix(thita)
n <- ncol(thita)
k <- nrow(thita)
B <- n*var(rowMeans(thita))
W <- mean(apply(thita,1,"var"))
Vhat <- (n-1)/n*W+B/n
PSRF <- sqrt(Vhat/W)
print(PSRF)
## [1] 1.00257
Gewekeの診断方法
- Geweke統計量Z値を計算し、MCMCのburn-in期間（前半10%）と後半50%を比較し、|Z|値が十分に小さければMCMCが収束したと考えます。以下のコードはパッケージcodaのgeweke.diag()関数を加筆修正しています。
c <- 1 # n.chain
n <- length(posterior[c,1,])
# MCMCの最初の10%をg1、後半50%をg2とする
p1 <- 0.1
p2 <- 0.5
n1 <- ceiling(1+p1*(n-1))
n2 <- floor(n-p2*(n-1))
g1 <- posterior[c,1,1:(n1)]
g2 <- posterior[c,1,n2:n]
# g1とg2の平均
g1.mean <- mean(g1)
g2.mean <- mean(g2)
# g1とg2のスペクトラム密度
x1 <- as.matrix(g1)
g1.spec <- numeric(ncol(x1))
z1 <- 1:nrow(x1)
for (i in 1:ncol(x1)) {
lm.out <- lm(x1[, i] ~ z1)
ar.out <- ar(x1[, i], aic = TRUE)
g1.spec[i] <- ar.out$var.pred/(1 - sum(ar.out$ar))^2
}
x2 <- as.matrix(g2)
g2.spec <- numeric(ncol(x2))
z2 <- 1:nrow(x2)
for (i in 1:ncol(x2)) {
lm.out <- lm(x2[, i] ~ z2)
ar.out <- ar(x2[, i], aic = TRUE)
g2.spec[i] <- ar.out$var.pred/(1 - sum(ar.out$ar))^2
order[i] <- ar.out$order
}
## Error in order[i] <- ar.out$order: object of type 'closure' is not subsettable
# Geweke統計量Z値の計算
g1.var <- g1.spec/length(g1)
g2.var <- g2.spec/length(g2)
z <- (g1.mean-g2.mean)/sqrt(g1.var+g2.var)
print(list(z=z, p=c(p1, p2)))
## $z
## [1] 0.9028606
##
## $p
## [1] 0.1 0.5
Raftery&Lewisの診断方法
- MCMCの標本期間が適切かどうかを判断します。以下のコードはパッケージcodaのraftery.diag()関数を加筆修正しています。
# Raftery-Lewis
data <- as.mcmc(posterior[1,1,])
# 初期値設定
q = 0.025 # 四分位偏差
r = 0.005 # 推定精度
s = 0.95 # 区間[q-r,q+r]
converge.eps = 0.001 # 収束判定条件
# Minimum sample size (Nmin)
phi <- qnorm(0.5*(1+s))
nmin <- as.integer(ceiling((q*(1-q)*phi^2)/r^2))
# 結果出力用行列
out <- matrix(nrow = 1, ncol = 4,
dimnames = list(varnames(data, allow.null = TRUE),
c("Burnin(M)", "Required sample size(N)",
"Minimum sample size(Nmin)", "Dependence factor(I=(M+N)/Nmin)")))
# 四分位偏差
quant <- quantile(data, probs = q)
# 四分位偏差を下回るMCMC 下回るとTRUE
dichot <- mcmc(data <= quant, start = start(data),
end = end(data), thin = thin(data))
# Thinning intervalの設定
kthin <- 0
# Transision matrixの計算
bic <- 1
while (bic >= 0) {
kthin <- kthin + thin(data)
testres <- as.vector(window(dichot, thin = kthin))
testres <- factor(testres, levels = c(FALSE, TRUE))
newdim <- length(testres)
testtran <- table(testres[1:(newdim - 2)],
testres[2:(newdim - 1)], testres[3:newdim])
testtran <- array(as.double(testtran), dim = dim(testtran))
g2 <- 0
for (i1 in 1:2) {
for (i2 in 1:2) {
for (i3 in 1:2) {
if (testtran[i1, i2, i3] != 0) {
fitted <- (sum(testtran[i1, i2, 1:2]) *
sum(testtran[1:2, i2, i3]))/(sum(testtran[1:2, i2, 1:2]))
g2 <- g2 + testtran[i1, i2, i3] * log(testtran[i1, i2, i3]/fitted) * 2
}
}
}
}
bic <- g2 - log(newdim - 2) * 2
}
# Transision Matrix
finaltran <- table(testres[1:(newdim - 1)], testres[2:newdim])
alpha <- finaltran[1, 2]/(finaltran[1, 1]+finaltran[1, 2])
beta <- finaltran[2, 1]/(finaltran[2, 1]+finaltran[2, 2])
tempburn <- log((converge.eps * (alpha + beta))/max(alpha,beta))/(log(abs(1 - alpha - beta)))
nburn <- as.integer(ceiling(tempburn) * kthin)
tempprec <- ((2 - alpha - beta) * alpha * beta * phi^2)/(((alpha + beta)^3) * r^2)
nkeep <- as.integer(ceiling(tempprec) * kthin)
iratio <- (nburn + nkeep)/nmin
# output matrix
out[1, 1] <- nburn
out[1, 2] <- nkeep + nburn
out[1, 3] <- nmin
out[1, 4] <- signif(iratio, digits = 3)
# print output
print(list(params=c(r=r, s=s, q=q), out=out))
## $params
## r s q
## 0.005 0.950 0.025
##
## $out
## Burnin(M) Required sample size(N) Minimum sample size(Nmin)
## [1,] 30 32295 3746
## Dependence factor(I=(M+N)/Nmin)
## [1,] 8.62