廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

独立メトロポリス・ヘイスティングス法を用いたベイズ推測の簡単な例題

なぜ MCMC が必要か

ベイズ統計学では母数 θ を確率変数と考える。

ベイズ統計のための準備, ベイズの定理, 事前分布と事後分布 - 廿TT

ベイズ推測はデータ \boldsymbol{x} が得られたもとでの θ の確率分布  f(\theta| \boldsymbol{x}) 事後分布を通じて行われる。

母数 θ の点推定には事後期待値(expected a posteriori ; EAP)推定量などが用いられる。

θ の期待値を求めるためには以下の積分計算が必要。

ベイズの定理と期待値の定義より、

 \displaystyle \int^{\infty}_{-\infty} f(\theta| \boldsymbol{x})d\theta= \int^{\infty}_{-\infty} \theta \frac{f( \boldsymbol{x} | \theta) f (\theta)}{ f(\boldsymbol{x})} d\theta

複数パラメータの場合は多重積分になる。

この難しい積分を避けるために事後分布=母数の確率分布から乱数を発生させて分布を調べる。

メトロポリスヘイスティングス

メトロポリスヘイスティングス法では, ターゲットとなる事後分布に従う乱数を、提案分布 q から発生させる。

R によるメトロポリス・ヘイスティング(MH)法入門 - 廿TT

そのアルゴリズムは大雑把に述べると、以下の通り。

「提案分布を利用し発生させた乱数 a を確率 \min(1,r) で受容(\theta_{\rm new} = a)し、さもなくば棄却(\theta_{\rm new} = \theta_{\rm old})する。」

ここで r は 尤度 L と提案分布(事前分布)の積

 \displaystyle \frac{q(\theta_{\rm old})L(a)}{q(a)L(\theta_{\rm old})}

である。

R で実験

事前分布を区間 [0, 10] の一様分布として、指数分布のパラメータの EAP 推定値を求める。

1万回サンプリングした乱数のトレースラインを描いてみる。

f:id:abrahamcow:20150918151341p:plain

グジグジの帯が真横に移動しているのが収束のサイン。

はじめの1000個をバーンイン期間と見積もり破棄してヒストグラムを描いてみる。

f:id:abrahamcow:20150918151746p:plain

青い線が EAP 推定値。赤い線が真値。

X <-rexp(10,2) #データ。パラメータ 2 の指数分布

1/mean(X) #最尤推定値

lik <- function(par){ #尤度
  prod(dexp(X,par))
}
proposal <- function(par){ #提案分布
  dunif(par,0,10)
}

M <- 10000 #試行回数
chain <- numeric(M)
acc <- 0
chain[1] <- runif(1,min=0,max=10)
lik(0)
for (i in 2:M){
  a <-runif(1,0,10)
  r <- (proposal(chain[i-1])*lik(a))/(proposal(a)*lik(chain[i-1]))
  tmp <- runif(1)
  if (tmp < r){
  chain[i] <-a
  acc<-acc+1
  }else{
  chain[i]<-chain[i-1]
  }
}

plot(chain,type="l",col="grey30")
abline(h=2,col="red3",lwd=3)

hist(chain[1001:10000],freq = FALSE)
abline(v=mean(chain[1001:10000]),col="blue2",lwd=2)

summary(chain[1001:10000])
acc/M #採択率

参考文献

基礎からのベイズ統計学: ハミルトニアンモンテカルロ法による実践的入門

基礎からのベイズ統計学: ハミルトニアンモンテカルロ法による実践的入門

追記

複数パラメータの分布(ガンマ分布)を推定する場合はこんな感じになる。

X <-rgamma(100,2,4) #データ

lik <- function(par){ #尤度
  sum(dgamma(X,par[1],par[2],log=TRUE))
}

optim(c(2,4),lik,control = list(fnscale=-1))

proposal <- function(par){ #提案分布
  dunif(par,0,10,log=TRUE)
}

M <- 10000 #試行回数
chain <- matrix(,M,2)
acc1 <- 0
acc2 <- 0
chain[1,] <- runif(2,min=0,max=10)
lik(0)
for (i in 2:M){
  a <-runif(2,0,10)
#  r1 <- exp((proposal(chain[i-1,1])+lik(a))-(proposal(a[1])+lik(chain[i-1,])))
  r1 <- exp((proposal(chain[i-1,1])+lik(a))-(proposal(a[1])+lik(chain[i-1,])))
  r2 <- exp((proposal(chain[i-1,2])+lik(a))-(proposal(a[2])+lik(chain[i-1,])))
  if (runif(1) < r1){
    chain[i,1] <-a[1]
    acc1<-acc1+1
  }else{
    chain[i,1]<-chain[i-1,1]
  }
  if (runif(1) < r2){
    chain[i,2] <-a[2]
    acc2<-acc2+1
  }else{
    chain[i,2]<-chain[i-1,2]
  }
}
head(chain)
plot(chain[,1],type="l",col="grey30")
abline(h=2,col="red3",lwd=3)
hist(chain[1001:10000,1],freq = FALSE)
abline(v=mean(chain[1001:10000,1]),col="blue2",lwd=2)
abline(v=2,col="red2",lwd=2)

plot(chain[,2],type="l",col="grey30")
hist(chain[1001:10000,2],freq = FALSE)
abline(v=mean(chain[1001:10000,2]),col="blue2",lwd=2)
abline(v=4,col="red2",lwd=2)

summary(chain)
acc1/M #採択率
acc2/M #採択率