Example of Metropolis Hastings Algorithm

ということで初めてのMCMCだよ。パターン認識と機械学習下 - ベイズ理論による統計的予測のP252-253、マルコフ連鎖モンテカルロ法 (統計ライブラリー)のP9-14とここの資料を参考にしながら作りました。

What to do

ガンベル分布(二重指数分布)からのサンプリングが困難な状況を仮定します。そういう状況でサンプリングしてきて、ガンベル分布の平均とか分散を計算してみたい。Rejection SamplingとかImportance Samplingとかあるけど、定数決めるのが難しかったり、そもそも多次元だと適切な提案分布を見つけるのが困難＞＜ということで、そういう状況でも戦えるマルコフ連鎖モンテカルロ法(MCMC)を導入します!!RにはMCMCpackとかいうライブラリはあることは知っているけど、最初なので敢えて自分で書いてみるよ!!

Metropolis Algorithm & Metropolis Hastings Algorithm

両アルゴリズムでも今までのように提案分布とかを使います。この辺の提案分布の決め方はまだよく分かってないんだけど、割りと自由っぽい(?)感じかなあ。この辺は「不偏分布が目標とする分布であるようなエルゴート的マルコフ連鎖を構成し、推移を繰り返せば、やがて目標分布を得ることができるというのがMCMCの概略である」と書いてある付近と関連しているのかしら。どういうものを提案分布として使えるのかはもうちょい勉強しないといけない。

Metropolis Algorithmでは提案分布がパラメータに関して対象ということを仮定しています。Metropolis Hastings Algorithmではその仮定がはずされています。今回はMetropolis Algorithmを使った例として(提案分布が対象という意味で)random walkを、Metropolis Hastings Algorithmではindependent chainというものを使っています。

Rejection Samplingのところとかでは棄却されたサンプルはただただ捨てられていたんだけど、上のアルゴリズムではちょっと違う感じの扱いになっています。

Result of sampling

ガンベル分布を
$f(x)=\frac{1}{2b}\exp(\frac{-|x-\mu|}{b})$
として、今回は $\mu=0$ 、 $b=3$ の時を考えてみます。この分布の平均と分散は解析的に得ることができて、平均 $\mu$ 、分散 $2b^2$ となります。これを知った上でそれぞれの結果を評価してみましょう。

Result of random walk

分散が大きくなるにつれて、採択率が下っていることが確認できます。精度はそこまで大きく違わないなーという印象ですです。

c	平均	分散	採択率
0.25	-0.4	27.65	97%
1	0.56	16.07	88%
5	-0.06	17.11	57%
50	-0.19	20.65	10%
150	0.56	17.72	3%

Result of independent chain

random walkの時とは違い、分散が小さいからといって、必ずしも採択率が大きいというわけではなさそうです。ピークがあるという感じ。あと、cによって分散の精度がやたら違うということが確認できる。

c	平均	分散	採択率
0.25	-0.36	4.4	57%
1	-0.15	3.12	62%
5	-0.03	17.26	76%
50	0.31	18.64	10%
150	-0.25	19.8	3%

Consideration

今回は10000サンプルを取ってきて、最初の1000サンプルを捨ててます。MCMC全般において、初期値の影響を受ける連鎖のはじめのT'までの繰り返しは推定に用いないそうです(この期間をバーンイン(burn-in)期間と呼ぶそうです)。で、本に書いてあるから1000捨てたわけですが、捨てる数をどうやって決定するかとかも考える必要がありそうです。

また、調整母数(tuning parameter)によって近似の精度がやたらと影響を受けていることが上の表からも分かりますが、この辺はどうやっていくのかも問題になりそう。まあ、たぶんここは最適化の時と同じように試行錯誤だろうな。

Program

こんな感じで必要な情報plotしちゃうよ!!

random.walk <- function(c,N){
  init <- 0
  xt <- rep(NA,(N+1))
  xt[1] <- init
  u <- runif(N)
  accept <- 0
  for (i in 1:N){
    x <- rnorm(1,xt[i],c)
    a <- min(1,
             exp(
                 - 1/3 * (abs(x) - abs(xt[i]))
                 )
             )
    if(u[i] < a){
      xt[i+1] <- x
      accept <- accept+1
    }else{
      xt[i+1] <- xt[i]
    }
  }
  return(list(accept=accept/N,
              xt=xt[((N/10)+1):(N+1)]))
}

independent.chain <- function(c,N){
  init <- 0
  xt <- rep(NA,(N+1))
  xt[1] <- init
  x <- rnorm(N,0,c)
  u <- runif(N)
  accept <- 0
  for (i in 1:N){
    a <- min(1,
             exp(
                 - 1/3 * (abs(x[i]) - abs(xt[i]))
                 - 1/2 * ((xt[i]/c)^2 - (x[i]/c)^2)
                 )
             )
    if(u[i] < a){
      xt[i+1] <- x[i]
      accept <- accept+1
    }else{
      xt[i+1] <- xt[i]
    }
  }
  return(list(accept=accept/N,
              xt=xt[((N/10)+1):(N+1)]))
}

plot.metropolis.hastings <- function(data,title){
  accept <- data$accept
  xt <- data$xt
  par(oma=c(0,0,2,0))
  par(mfrow=c(1,2))
  title(main=title)
  plot(density(xt),
       main=paste("密度推定\n",
         "mean=",round(mean(xt),digits=2),
         ",var=",round(var(xt),digits=2)
         )
       )
  plot(xt,type="l")
  title("Iteration")
  mtext(side=3,line=-0.75,outer=T,text=paste(title,"\t(採択率:",round(accept*100),"%)",sep=""),cex=2)
}

plot.new()
r <- random.walk(1,10000)
plot.metropolis.hastings(r,"酔歩連鎖")
r <- random.walk(5,10000)
plot.metropolis.hastings(r,"酔歩連鎖")
r <- random.walk(50,10000)
plot.metropolis.hastings(r,"酔歩連鎖")

r <- independent.chain(1,10000)
plot.metropolis.hastings(r,"独立連鎖")
r <- independent.chain(5,10000)
plot.metropolis.hastings(r,"独立連鎖")
r <- independent.chain(50,10000)
plot.metropolis.hastings(r,"独立連鎖")