基礎からのベイズ統計学（豊田秀樹著、朝倉書店）からの学習ノート
ーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーー
理解を深める核心的な重要な記述：マルコフ連鎖モンテカルロ法のメトロポリス・ヘイスティング法についての解説

4.1 事後分布からの乱数の発生
ベイズ統計学では、データからの知見を事後分布を経由して入手する。しかしながら、事後分布の解析的な評価は実践的な場面では困難を極める。２項分布モデルやポアソン分布モデルなどの事後分布は、理論的道筋を明解にするために教材として利用されているが、あくまで教材であり、実践の場ではそんなに簡単に当てはめられるようなものはないということ。

そこで発想を転換して、事後分布f(θ|x)は母数の確率分布であるので、その確率分布から乱数を発生させる。
伝統的な統計学では、母数θが固定され、観測変数は確率分布に従い、観測変数が実現値として観察される。しかし、ここでは逆で、観測変数が固定され、確率変数である母数θが事後分布に従い、母数の実現値が観測される。

しかし、そう簡単に事後分布に従う乱数列を手にすることはできない。
たとえば、正規分布に従う正規乱数は、計算機言語に実装されているf(x|μ,σ）である。一方で、ここでは、f(μ,σ|x)が求められる。

加えて、ベイズの定理　f(θ|X) = f(x|θ)f(θ)/f(x)では、 事後分布f(θ|X) ∝ カーネル：f(x|θ)f(θ)
として、分母を正規化定数としてカットして利用してきたが、事後分布を求めるためには、分母f(x)が必要であり、f(x) = ∫f(x|θ)f(θ)dθの積分値であり、この正規化定数が不明であれば確率分布は求められない。

そこでどうやって、事後分布のカーネルだけを用いて、事後分布f(θ|x)に従う乱数を発生させるか、が課題となる。

4.2 マルコフ連鎖
ネクタイ問題：ある高校の制服では、三種類のネクタイA B Cが用意され、100人の生徒が入学初日には0.6, 0.25, 0.15の比率でそれらのネクタイを着用する。それ以後は、前日に締めたネクタイに応じてだけ次の日のネクタイが選ばれる。
ネクタイ変更の条件付き確率を以下のように規程する。遷移核という。
　　　 　　当日(j)
　　 　 　　　A　 B　 C
　　　　　　 A 　　 0.3 0.3 0.4
前日(i)　　 B 　　0.1 0.5 0.4
　 C 　　0.2 0.6 0.2

二日目は、Aのネクタイは、0.3 x 0.6 + 0.1 x 0.25 + 0.3 x 0. 15 = 0.235
同様にBのネクタイは　0.395、Cのネクタイは0.37と計算できる。
計算を続けると、1/6:1/2:1/3で落ち着くことがわかる。マルコフ連鎖の定常分布stationary distribution、もしくは不変分布invariant distributionという。
ここでは初期条件{0.6:0.25:0.15}を変化させて、{0.3:0.3:0.4}や、{0.1:0.1:0.8}としても同じ比率に収束することがわかる。

サンプリングしたい分布（目標分布）に対して、それを定常分布とするマルコフ連鎖を構成する方法をマルコフ連鎖モンテカルロ法（MCMC)という。MCMCでは、定常分布が既知であり、遷移核が未知である。

4.4 詳細釣り合い条件detailed balance condition
マルコフ連鎖が定常分布に収束するための十分条件：
p(X=i|X’=j)p(X’=j) = p(X=j|X’=i)p(X’=i)

この式の両辺を添え字iに対して和をとり
Σ(p(X=i|X’=j)p(X’=j)) = Σ(p(X’=j|X=i)p(X=i))
左辺にΣp(B|A) = 1を適応すれば
p(X’=j)) = Σ(p(X’=j|X=i)p(X=i))

p(X’=j|X=i)=遷移核なので、

p(X’=j)) = Σ(遷移核 * p(X=i))

母数θは連続型の確率変数であることがほとんどなので、詳細釣り合い条件を確率密度関数で置き換えれば
f(θ|θ’)f(θ’) = f(θ’|θ)f(θ)
目標分布がf(θ)であり、遷移核がf(θ’|θ)
逆に 目標分布がf(θ’)であり、遷移核がf(θ|θ’)と言っても良い。

f(θ) = f(θ|θ’)f(θ’)/f(θ’|θ) = ∫f(θ|θ’)f(θ’)θ’
発射地点θ’からθに飛んでくる確率密度のあらゆる発射地点に関する平均確率密度（右辺）がθの確率密度そのもの（左辺）であることを意味している。

4.5 メトロポリス・ヘイスティング法
前章では
f(θ|θ’)f(θ’) = f(θ’|θ)f(θ)
目標分布がf(θ)であり、遷移核がf(θ’|θ) であったが、
既知な事後分布f(θ)に対して、上記式を満たすような遷移核f(|)を見つけることは困難である。
そこで遷移核f(|)の代わりに適当な遷移核q(|)を提案分布(proposal distribution)として利用する。
適当なので
q(θ|θ’)f(θ’) > q(θ’|θ)f(θ)
というようにもなる。この式を釣り合い条件に向けて確率補正する方法がメトロポリス・ヘイスティング法である。

ここで遷移確率との補正を行うために、符号が正の未知の補正係数cとc’を導入する
f(θ|θ’) = c q(θ|θ’)
f(θ’|θ) = c’ q(θ’|θ)
すると補正後は
c q(θ|θ’)f(θ’) = c’ q(θ’|θ)f(θ)
となる。
ここでの数学的問題は、一つの方程式にcとc’という2つの変数が含まれていることと、補正係数が0〜1の範囲内で収まらないと確率的行動が取りにくいことである。そのために両辺をc’で除算し、c/c’で方程式を解く。

r =c/c’ = q(θ’|θ)f(θ) / q(θ|θ’)f(θ’)
r’ = c’/c’ = 1

r q(θ|θ’)f(θ’) = r’ q(θ’|θ)f(θ)

rq(θ|θ’)とr’q(θ’|θ)を遷移核として採用すれば、詳細釣り合い条件を満たす遷移が実現できる。

MHアルゴリズム
1) 提案分布q(|θ(t))を利用し、乱数αを発生する。
2) 以下の命題を判定する。
q(α|θ(t))f(θ(t)) > q(θ(t)|α)f(α)

[真] 不等号条件により、θ=α, θ’=θ(t)の状態と判定される。この場合、提案分布はq(θ|θ’)であるから、
r = q(θ(t)|α)f(α) / q(α|θ(t))f(θ(t))
を計算し、確率rでαを許容し、θ(t+1)=αとする。確率1-rでαを破棄し、θ(t+1)=θ(t)とする。

[偽]　不等号条件により、θ’=α、θ=θ(t)の状態と判定される。この場合、提案分布はq(θ’|θ)であるから、係数r’を使って確率的補正をする。確率１(=r’)で必ずαを受容し、θ(t+1)=αとする。

39 t = t + 1として、1)へ戻る

まとめると、MH法は、提案された候補点αを確率min(1,r)で受容(θ(t+1)=α)とし、さもなくばその場に留まる(θ(t+1) =　θ(t))ことを繰り返す。

ここで重要なこと：
r = q(θ(t)|α)f(α) / q(α|θ(t))f(θ(t))
分母、分子はそもそも事後分布の式となっている。要するに提案分布と事後分布の比を取ることで、
式の中には積分値である正規化定数が消去される。z乱数の発生は困難であるが値の評価は容易な事後分布のカーネルと、乱数発生が容易な提案分布を利用して、事後分布から母数の実現値の列（擬似乱数）列を入手する方法である。