マルコフ連鎖モンテカルロ法

ベイジアン統計解析の実際(丹後俊郎、Taeko Becque著、朝倉書店)からの学習ノート
——————————————————-
マルコフ連鎖モンテカルロ法
3.3 Metropolis-Hastingアルゴリズム
解析的に表現できない事後分布がπ(θ)を定常分布にもつ

π(θ) = p(θ|y) = p(y|θ)p(θ) / p(y)

このマルコフ連鎖 p(・|・)をどのように構成するか?
これがMetropolisらがその方法を提示し、Hastingがアルゴリズムを完成させた以下のMetropolis-Hastingアルゴリズムである。

Metropolis-Hastingアルゴリズム:
1) まず、解析者が現時点の値(状態)に基いて、次の時点にどの値(状態)に推移するかを示す条件付き分布を提案する。これをサンプラー(sampler)と呼ぶ
2)次に、このサンプラーから一つのサンプル(乱数)を発生させて、次に進むべき候補を示す。
3)この値が採択されればその値に進むが、採択するかどうかは、カギとなる確率式(後で示す採択確率)で与えられる。

4)この値が棄却されれば、現時点と同じ値にとどまる。
5)このプロセスを数多く繰り返せば、発生された値に関するヒストグラムが生成され、繰り返しを増やすことでヒストグラムがなめらかな曲線に収束して定常分布、つまり求める事後分布π(θ)の密度関数となる。
*このアルゴリズムの重要な性質は、これから得られる定常分布が、提案されるサンプラーにかかわらず、π(θ)に一致することである。

採択確率acceptance probability:

α(θ(t), z) = min (1, π(z)q(θ(t)|z) / π(θ(t)q(z|θ(t)))

θ(t+1) = z, with probability α(θ(t), z)
θ(t+1) = θ(t), with probability 1 – α(θ(t), z)

ちょっとわかりにくすぎるので、θ(t)=aとして、採択される場合θ(t+1) z = b を考えてみると

min(1, π(θ(t+1))q(θ|θ(t+1)) / π(θ(t))q(θ(t+1)|θ(t))

min(1, q(a|b)π(b)) / q(b|a)π(a))

事後分布p(θ|y)の計算には積分p(y) = ∫p(θ|y)p(θ)dθが含まれているが、この積分は計算が困難。
しかしながら、q(a|b)π(b)) / q(b|a)π(a)で分母分子でキャンセルされる。

π(θ) = p(θ|y) = p(y|θ)p(θ)/p(y)
p(y) = ∫p(a|y)p(a)da
p(y) = ∫p(b|y)p(b)db

事後分布の核(kernel):π(θ) ∝ p(y|θ)p(θ)だけを指定すればよい。ベイジアン統計解析では、事後分布の核の関数形は既知であるのでMHアルゴリズムはたいへん便利となる。

とにかくよく解らないが、感覚的には、以下のように理解しておくこととする。
山登りをしていて、帰りの道に迷った。夜間で景色は見えない。少し歩いてみて、高いところに行けばダメとし、低いところに行けば良しとする。どちらの方向に行くかは色々な方法がある。例えばサイコロで決めるとか、ぐるっと回って止まったところの方向で進むとか。でも方法にこだわらずに、とにかく進んで低い所を繰り返せば、どんな方法で進む方向を決めようとも、そのうちに低地の町にでれるのではないか。。。。というようなことで理解しておくとする。

具体的な導入では、
1) sampler 1(・|・)の候補はたくさんある。
2) 早く収束(rapid mixing)するものがよいが、定常分布π(・)との関係に依存する。

条件付き確率として代表的なものは以下の2つである。

1)対称サンプラー Metropolis sampler
q(z|θ) = q(θ|z)
この場合、採択確率は、min(1, π(z)/π(θ))
応用では、正規分布 q(z|θ) = N(θ, σ^2)
ランダムウォークサンプラー: q(z|θ) = q(|z-θ|)も特殊な例としてある。

2) 独立サンプラー Independence sampler
q(z|θ) = q(z)
で状態に依存しないサンプラー
この場合、採択率はα(θ, z) = min(1, w(z)/w(θ)), w(・)=π(・)/q(・)

このあと
3.4 単一成分Metropolis-Hastingアルゴリズム
3.5 Gibbsサンプリング
3.6 欠損データの補完
3.7 収束診断
と続くが、これ以上は本書の内容では、とても理解できないのであきらめる。