読者です 読者をやめる 読者になる 読者になる

廿TT

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

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

R 確率過程

はじめに;動機など

メトロポリス・ヘイスティング(Metropolis-Hastings)法はマルコフ連鎖モンテカルロ(Marcov Chain Monte Calro; MCMC)法の一種である.

MCMCでは, 独立でも同一でもない分布, マルコフ連鎖により乱数を生成する.

MCMCの代表的な手法(アルゴリズム)はMH法とギブス・サンプラーがある. ここで扱うのはMH法のみである.

MH法には主に, 

  1. ベイズ統計の枠組みで用いる
  2. 乱数生成法として用いる

のふたつの使い方がある.

最近流行っているのは前者だが, この記事は後者の観点から書かれている.

ぼくはいま, マルコフ連鎖がどのくらいで定常過程に収束するのかに興味がある.

アルゴリズム;MH法(連続型)

MH法では, ターゲットとなる密度関数 π に従う乱数を得るために, 別の密度関数 p(\cdot|x) に従う乱数を用意する.

この p(\cdot|x) を提案分布と呼ぶ(transition density ともいう).

x は分布のパラメータである.

また,

 \displaystyle \alpha(x,y) = \min \left( \frac{\pi(y) p(y|x)}{\pi(x) p(x|y)}, 1 \right)

とする.

アルゴリズムは以下の通りである.

  • p(\cdot|x) から初期値  x_0 を発生させる
  • 以下を好きなだけ繰り返す(j = 1,2,3,\ldots, N
    1. p(\cdot|x_{j-1}) から y_j を発生させる
    2. 区間 [0, 1] からの一様乱数 U を発生させる
    3. もし  U \le \alpha(x_{j-1},y_j) ならば x_{j+1}=y_j, さもなくば  x_j = x_{j-1} とする.
  • x_jを返す

(棄却サンプリングの手法を用いている. 
 R による棄却サンプリング法(おれが)入門 - 廿TT

数値実験; 正規分布

平均 0, 分散 1 の正規乱数から, 平均 10, 分散 1 の正規乱数を生成する.

pi<-function(x, log=FALSE) {
  dnorm(x, 10, 1, log=log)
}
p <- function(x,y,log=FALSE){
  dnorm(y, x, log=log)
}

alpha <- function(x, y) {
  return(exp(pi(y,log=TRUE)+ p(y,x,log=TRUE) - 
               pi(x, log=TRUE) - p(x, y, log=TRUE)))
}

MCMC <- function(X0, N) {
  X<-numeric(length=N) 
  Xj <- X0
  for(j in 1:N){
    Yj <-rnorm(1, Xj)
    U <- runif(1)
    if( U < alpha(Xj,Yj) ){
      Xj <- Yj
    }
    X[j] <- Xj
  }
  return( c(X0, X) )
}

N <- 1000
j <- 0:N
X <- MCMC(0, N)
plot(j,X, type="l", ylab=expression(X[j]) ) 

f:id:abrahamcow:20150118132402p:plain

だいたい j = 200 くらいで定常分布に収束しているようだ.

j = 1, ...,200 を捨ててヒストグラムとQ-Qプロットを書いてみる.

こうやって捨てる区切りとなる点 200 のことを burn in(バーンイン;焼きなまし)periodと呼ぶ.

f:id:abrahamcow:20150118132655p:plain

hist(X[-c(1:200)], freq=FALSE)
curve(dnorm(x,10,1),add=TRUE,lty=2)


f:id:abrahamcow:20150118132740p:plain

qqnorm(X[-c(1:200)])
qqline(X[-c(1:200)],col="red3",lwd=2)

補足;ベイズ統計

ベイズ統計の枠組みでは, p を事前密度関数, π を事後密度関数と呼ぶ.

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

今後の課題

定常過程への収束判定はどうやるのだろうか.

参考文献

計算統計学の方法―ブートストラップ・EMアルゴリズム・MCMC (シリーズ予測と発見の科学 5)

計算統計学の方法―ブートストラップ・EMアルゴリズム・MCMC (シリーズ予測と発見の科学 5)