Metropolis Algorithm の実装例

「入門ランダムウォーク：Metropolis Algorithm」からの続き。


Metropolis algorithm の実装例の前。このアルゴリズム、つまり MCMC 等のランダムウォークが、すべての予測に適切ではないことを述べているのが以下の引用

The cost of introducing the Metropolis algorithm in this simple one-parameter scenario is that you may come away with the false belief that people would actually do this in real data analysis. They would not. For a single parameter with a limited range, grid approximation can be more accurate than MCMC, and grid approximation allows you to estimate p(D) easily.

限定された幅の 1 つのパラメータにおいて、離散型の近似の方が MCMC より正確で、p(D) の見積りも容易。

上図は、同じ観測データと事前確率で、左は Figure 5.2 の連続型（数学的）で求めた事後確率、右は Figure 6.2 の離散型で求めた事後確率。同じ結果となった。

左が、本書提供のスクリプト BernMetropolisTemplate.R で求めたもの。これも上記二つを同じと判断できる。

Metropolis algorithm を実装したこの R スクリプトを元に、実装の詳細を見る。


P_{move：動く確率}

前回の Metropolis algorithm の紹介では、数学的な表記は避けたが、ここで一つだけ取り上げる。

前回の例で 4 の位置から 5 の提案があった場合、θ_proposed = 5, θ_current = 4 とすると、動く確率は

　P_move = min( θ_proposed / θ_current , 1 ) = min( 5/4, 1 ) = 1

つまり「100% 動く」ことになる。これを一般化したのが以下の式

これは 7.1 式で「移動する確率」。

min はその引数の内で小さい方を返すので、「提案された θ の確率」が「現在の θ の確率」より大きい場合は、移動する確率が 1.0 、つまり「必ず移動する」。小さい場合は、その割合が動く確率になる。

次は、事前確率を β 分布、尤度を Bernoulli 関数の場合（"7.3 The metropolis algorithm more generally" より）。

ランダムウォークの実装

BernMetropolisTemplate.R の要点を解説する。

#尤度
likelihood = function( theta , data ) {
　z = sum( data == 1 )
　N = length( data )
　pDataGivenTheta = theta^z* (1-theta)^(N-z) #Bernoulli likelihood
　pDataGivenTheta[ theta > 1 | theta < 0 ] = 0
　return( pDataGivenTheta )
}
#事前確率
prior = function( theta ) {
　prior = rep( 1 , length(theta) ) # uniform density over [0,1]
　prior[ theta > 1 | theta < 0 ] = 0
　return( prior )
}

targetRelProb = function( theta , data ) {
　targetRelProb = likelihood( theta , data ) * prior( theta )
　return( targetRelProb )
}

ここまでが 尤度、事前確率、「尤度×事前確率」の関数の定義。以降は、これらの関数を使ったロジックとなる。


myData = c( 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0 )
proposedJump = rnorm( 1 , mean = 0 , sd = 0.1 )
probAccept = min( 1,
　targetRelProb( currentPosition + proposedJump , myData )
　/ targetRelProb( currentPosition , myData ) )

これは 7.1式 の実装。提案後の θ の尤度（分子）が、提案前の θ の尤度（分母）より大きければ probAccept は 1 となる。今回の prior は常に 1 を返すため、targetRelProb は尤度 θ^z * (1-θ)^(N-z) の Bernoulli Likelihood となる。


if ( runif(1) < probAccept ) {
　trajectory[ t+1 ] = currentPosition + proposedJump
　if ( t > burnIn ) { nAccepted = nAccepted + 1 }
} else {
　trajectory[ t+1 ] = currentPosition
　if ( t > burnIn ) { nRejected = nRejected + 1 }
}

先のロジックで probAccept（提案を受入れる確率）が runif（0 から 1 までの乱数）より大きければ提案の確率を trajectory（「軌道点」）に追加、小さければ現在の θ を追加。


trajLength = 55556
burnIn = ceiling( 0.1 *burnIn = ceiling( 0.1 * trajLength ) 
acceptedTraj = trajectory[ (burnIn+1) : length(trajectory) ]

今回は burn-in の回数を、総「軌道点」数の約 10% としたので、acceptedTraj（有効な「軌道点」）には 10% 以降の約 90% が入る。

以上が Metropolis algorithm のランダムウォークで、ここまでに格納した変数値を使って、左のグラフ（冒頭のグラフと同じ）生成。ただし、p(D) の算出方法は後述。

mean：事後確率分布の平均値

> mean(acceptedTraj)

[1] 0.7489788

N_pro：総軌道点数
> length(acceptedTraj)
[1] 50000

N_acc/N_pro：総軌道点数における提案を受入た割合、すなわち「移動数の割合」。

> signif(nAccepted/length(acceptedTraj))

[1] 0.71804

Evidence p(D) の算出

ここでは面倒な数学的な解説は省きたいので、今回は本書から要素だけを書く。

左は式 10.8 。一行目はベイズルール p(θ | D) = p(D | θ)p(θ) / p(D) を変形しただけ。トリッキーなのは二行目で

　∫ dθ h(θ) = 1

を右辺に掛けている。1 の掛算なので左辺に変化はなく、右辺を 3, 4 行目のように変形することができる。

最後の 4 行目は h(θ) / p(D | θ)p(θ) の平均値で、この p(D) の算出の実装例は

wtdEvid = dbeta( acceptedTraj, a, b ) / (

likelihood( acceptedTraj, myData ) * prior( acceptedTraj ) )

pData = 1 / mean( wtdEvid )

ポイントは分子 h(θ) にベータ関数を使っていること。そうする理由は

• h(θ) は likelihood(θ) × prior(θ) の類似であること
• likelihood(θ) が binominal distribution の場合は、 h(θ) は ベータ分布にすべき

以下は、スクリプト中の「ベータ分布の求め方」

Compute a, b parameters for beta distribution that has the same mean and stdev as the sample from the posterior. This is a useful choice when the likelihood function is Bernoulli.

ベータ分布の平均と標準偏差が、事後確率の標本の平均と標準偏差と同じ、引数 a, b とする。

よって a, b の算出は

meanTraj = mean( acceptedTraj )
sdTraj = sd( acceptedTraj )

a =   meanTraj   * ( (meanTraj*(1-meanTraj)/sdTraj^2) - 1 )

b = (1-meanTraj) * ( (meanTraj*(1-meanTraj)/sdTraj^2) - 1 )

「実習7.1-2 ランダムウォークを観察」に続く。

/** tex **/