実習7.1-2 ランダムウォークを観察

「Metropolis Algorithm の実装例」からの続き。


実習7.1 本書提供のスクリプト BernMetropolisTemplate.R の proposedJump の乱数生成で、正規分布関数の標準偏差を 0.1, 0.001, 100.0 で、事後確率を比較する。

変更するスクリプトの箇所は引数 sd（標準偏差）

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

この変更の影響は、currentPosition + proposedJump が sd の大小に比例すること。例えば、sd = 100.0 の場合は 0.1 よりも提案される値が大きく、0.001 の場合は小さくなる。

以下が実行例

SD = 0.001 の「提案受入れ率」は 0.996 で、SD = 0.1 の 0.712 より高い。

これは、小さい SD の場合は現在位置からの変化が少なく、左の式のように P_move の値が高くなる。つまり、提案を受入れやすい値ばかりでは target distribution は反映しにくくなるということ。

逆に SD = 100.0 と大きい場合は、P_move の分子、つまり提案値が分母（現在値）より小さくなるケースが増え、P_move の値が小さくなることで提案を拒絶されがちになる。提案の受入率が 0.0024 という低さでは、到底 target distribution は反映できない。


実習7.2 スクリプト BernMetropolisTemplate.R に以下の変更をした結果を分析。

• trajlength = 100：ランダムウォーク数を極端に減らす
• burnin = ceiling(0.01 * trajlength)：burn-in の期間を極端に短縮
• trajectory[1] = 0.001：ランダムウォーク開始点を極端に小さく

変更の結果、ランダムウォークは左端 0.01 から始まり、burn-in は 1 ステップで提案数は 99 、この結果が左図。

θ の極端に小さい値からスタートし、且つ標準偏差が 0.1 のため、0.5 以上に達するまでにステップ数を多く要している。

ランダムウォークの様子を知るには、ヒストグラムよりも左のプロットが有効。

> plot(acceptedTraj,1:length(acceptedTraj),
 type="o",ylab="Step",xlab="Parameter Value")


θ = 0.01 から徐々に 0.5 以上の地点に達し、それ以降は 0.5 を下回っていない。

本実習で分かるのは、乱数の開始地点を的確に設定するのは難しいが、burn-in 期間を設けることでその困難を補えるということ。

「実習7.3-4 複数山はムリ、h(θ) は分母の近似」に続く。