【ベイズ統計学】入門　ベイズ統計学（05/X） - 「大人の教養・知識・気付き」を伸ばすブログ

　Bayes統計学を今回は

基礎からのベイズ統計学: ハミルトニアンモンテカルロ法による実践的入門

作者:豊田秀樹
朝倉書店

Amazon

を主に参照しつつ学んでいく。

4.　Metropolis-Hastings法

4.　Metropolis-Hastings法

4.1　事後分布からの乱数発生

　Bayes統計学では、データからの知見を事後分布経由で得る。しかし事後分布の解析的評価は不可能な局面が殆どである。たとえばEAP推定量

$\begin{aligned} \hat{\theta}_{EAP}=\displaystyle{\int\theta f(\theta|\boldsymbol{x})}d\theta \end{aligned}$

は母数に関する積分が必要となるが、積分は評価できるものと出来ないものが存在する。特に多変量への積分になれば、計算負荷が大きくなり、数値的にも不可能という状態になり得る。そのため発想を変え、事後分布が母数の確率分布であることを踏まえて、事後分布から乱数を発生させることにする。すなわち観測変数が固定され、確率変数である母数 $\boldsymbol{\theta}$ が事後分布に従い、母数の実現値も観測されているという状況を考える。
　ただし2点問題がある。

仮に確率分布が評価できても、必ずしも当該確率分布に従う乱数を発生させられるとは限らない。
事後分布に含まれっる正規化定数は評価できない場合が多い。

これらを考慮して、カーネルのみを用いて事後分布に従う乱数を発生させることを考える。

4.2　Monte Carlo積分 *1

　積分の計算方法には様々なものがあるが、1つにはMonte Carlo法を用いるものがある。
　いま $\boldsymbol{\theta}$ の確率密度関数を $s(\boldsymbol{\theta})$ として以下の積分を考える。

$\begin{aligned} \gamma=\displaystyle{\int h(\boldsymbol{\theta})s(\boldsymbol{\theta})d\boldsymbol{\theta}} \end{aligned}$

　仮に $s(\boldsymbol{\theta})$ から独立に $\boldsymbol{\theta}^{(j)},j=1,\cdots,L$ が得られたとすると、　

$\begin{aligned} \hat{\gamma}=\displaystyle{\frac{1}{L}\sum_{j=1}^{L}h\left(\boldsymbol{\theta})^{(j)}\right)} \end{aligned}$

は大数の強法則から、これは $L\rightarrow\infty$ で $L$ に収束する。これをMonte Carlo積分という。
　Monte Carlo積分を用いてさまざまな統計量を（近似的に）計算する。

(1)事後平均

$\begin{aligned} \bar{\boldsymbol{\theta}}_n=\displaystyle{\int\boldsymbol{\theta}\pi(\boldsymbol{\theta}|\boldsymbol{X}_n)}d\boldsymbol{\theta} \end{aligned}$

を近似すると、

$\begin{aligned} \displaystyle{\frac{1}{L}\sum_{j=1}^{L}\boldsymbol{\theta}^{(j)}} \end{aligned}$

が該当する。

(2)事後モード

$\begin{aligned} \hat{\boldsymbol{\theta}}_n=\displaystyle{\arg\max_{\boldsymbol{\theta}}}{\pi(\boldsymbol{\theta}|\boldsymbol{X}_n)} \end{aligned}$

を近似すると、

$\begin{aligned} \displaystyle{\arg\max_{j}}{\ \pi(\boldsymbol{\theta^{(j)}|\boldsymbol{X}_n )}} \end{aligned}$

が該当する。

(3)特定の領域 $Q$ に入る確率

$\begin{aligned} \pi(\boldsymbol{\theta}\in Q|\boldsymbol{X}_n) \end{aligned}$

を近似すると、

$\begin{aligned} \displaystyle{\frac{1}{L}\sum_{j=1}^{L}\boldsymbol{1}\left(\boldsymbol{\theta}^{(j)}\in Q\right)} \end{aligned}$

が該当する。

(4)周辺事後分布

$\begin{aligned} \pi(\boldsymbol{\theta}_1|\boldsymbol{X}n)=\displaystyle{\int\pi(\boldsymbol{\theta}_1,\boldsymbol{\theta}_2|\boldsymbol{X}_n)}d\boldsymbol{\theta}_2 \end{aligned}$

を近似すると、

$\begin{aligned} \displaystyle{\frac{1}{L}\sum_{j=1}^{L}\pi(\boldsymbol{\theta}_1,\boldsymbol{\theta}^{(j)}_2|\boldsymbol{X}_n)} \end{aligned}$

が該当する。

(5)予測分布

$\begin{aligned} f(\boldsymbol{z}|\boldsymbol{X}_n)=\displaystyle{\int f(\boldsymbol{z}|\boldsymbol{\theta})\pi(\boldsymbol{\theta}|\boldsymbol{X}_n)}d\boldsymbol{\theta} \end{aligned}$

を近似すると、

$\begin{aligned} \displaystyle{\frac{1}{L}\sum_{j=1}^{L}f(\boldsymbol{z}|\boldsymbol{\theta})^{(j)})} \end{aligned}$

が該当する。

4.3　Markov連鎖

　カーネルのみを用いて事後分布に従う乱数を発生させるために、Markov連鎖の利用を考える。
　一般の確率過程に対して、

$\begin{aligned} P(X^{(t)}|X^{(t-1)},\cdots,X^{(1)})=P(X^{(t)}|X^{(t-1)}) \end{aligned}$

すなわち条件付き確率が直前の時点における状況のみの影響を受けるような確率過程をMarkov連鎖という。Markov連鎖を規定する条件付き確率を遷移核(transition kernel)という。

4.4　定常分布への収束

　Markov連鎖を用いて確率を遷移させていった結果として変化しなくなった確率分布 $\boldsymbol{p}$ をそのMarkov連鎖の定常分布(または不変分布)という。この挙動を定常分布への収束といい、収束までの期間をバーンイン期間 $B$ という。
　Markov連鎖は遷移カーネルと初期状態が、

既約的：有限回の推移で状態空間の要素全てが互いに到達可能であること
正再帰的：状態空間の任意の要素は限りなく何度も遷移し得ること
非周期的：連鎖の状態が一定の周期性を持たないこと

を有するとき、定常分布に収束することが保証される。
　定常分布には以下が知られている：

Markov連鎖 $\{X_n\}_{n=1}^{\infty}$ で
$\begin{aligned}p(\theta^{\prime})=\displaystyle{\int _{-\infty}^{\infty}p(\theta)f(\theta|\theta^{\prime})}d\theta\end{aligned}$
を満たすような定常分布 $p(\theta)$ が存在する。
Markov連鎖 $\{X_n\}_{n=1}^{\infty}$ で

4.5　Markov連鎖Monte Carlo法

　サンプリングしたい分布を定常分布とするようなMarkov連鎖を構成する方法をMarkov連鎖Monte Carlo法と総称する。Markov連鎖Monte Carlo法ではサンプリングしたい分布を目標分布と呼ぶ。このときはすなわち、定常分布が既知である中で遷移核を知ることが問題になる。
　Markov連鎖Monte Carlo法は以下の3点から爆発的に普及した：

MCMC法の発想は極めて単純でコンピュータ上にて容易に実行可能である
従来の手法では扱えなかった広範囲の複雑なモンデるに対しても適用できる拡張性を有している
古典統計学の手法では扱えない一方でMCMC法に依るBayes分析であれば扱えるモデルが存在する

4.6　詳細釣り合い条件

　Markov連鎖が定常分布に収束するための十分条件として詳細釣り合い条件がある。

　標本空間のすべての事象の組 $i,j$ に対し
$\begin{aligned} p(X=i|X^{\prime}=j)p(X^{\prime}=j)=p(X^{\prime}=j|X=i)p(X=i) \end{aligned}$
が成り立つこと。

4.5.1　詳細釣り合い条件の意味

　詳細釣り合い条件において $i$ に関する総和を両辺について取ると

$\begin{aligned} &\displaystyle{\sum_{i=1}^{n}p(X=i|X^{\prime}=j)p(X^{\prime}=j)}=\displaystyle{\sum_{i=1}^{n}p(X^{\prime}=j|X=i)p(X=i)}\\ \Leftrightarrow&p(X^{\prime}=j)=\displaystyle{\sum_{i=1}^{n}p(X^{\prime}=j|X=i)p(X=i)} \end{aligned}$

が成り立つ。これは分布が変化しないように制約を入れた全確率の公式である。
　これを連続な確率変数に関して書けば

$\begin{aligned} f(\theta|\theta^{\prime})f(\theta^{\prime})=f(\theta^{\prime}|\theta)f(\theta) \end{aligned}$

である。ここで $f(\theta)$ が目標分布であり、 $f(\theta^{\prime}|\theta)$ が遷移核である。目標分布が $f(\theta^{\prime})$ であり遷移核が $f(\theta|\theta^{\prime})$ であると言っても同じである。
　この式はあらゆる点において成り立つから、これらを積分しても成り立つ：

$\begin{aligned} &\displaystyle{\int f(\theta|\theta^{\prime})f(\theta^{\prime})d\theta^{\prime}}=\displaystyle{\int f(\theta^{\prime}|\theta)f(\theta)d\theta^{\prime}},\\ \Leftrightarrow&f(\theta)=\displaystyle{\int f(\theta|\theta^{\prime})f(\theta^{\prime})d\theta^{\prime}} \end{aligned}$

これは $\theta^{\prime}$ から $\theta$ に遷移する確率密度のあらゆる開始地点に関する平均確率密度（右辺）が $\theta$ の確率密度そのものであると言っている。したがって詳細釣り合い条件が満たされているならば、初期状態を真の値からデタラメに遠くに取ったとしてもその周辺へと乱数列が急速に近づく。

4.6　Metropolis-Hastings法

　既知である事後分布 $f(\theta)$ に対して詳細釣り合い条件

$\begin{aligned} f(\theta|\theta^{\prime})f(\theta^{\prime})=f(\theta^{\prime}|\theta)f(\theta) \end{aligned}$

を満たすような遷移核 $f(\theta^{\prime}|\theta)$ を見つけるのは困難である。そこで代わりに適当な遷移核 $q(\cdot|\cdot)$ を提案分布として用いる。
　とはいえ適当に選択する以上、詳細釣り合い条件を満たすとは限らない。そこでそれを満たすような方向へ確率補正を掛けるのがMetropolis-Hastings法である。
　確率補正の方法として適当な正の定数

$\begin{aligned} r=\displaystyle{\frac{q(\theta^{\prime}|\theta)f(\theta)}{q(\theta|\theta^{\prime}f(\theta^{\prime}))}} \end{aligned}$

を導入し

$\begin{aligned} r q(\theta|\theta^{\prime})f(\theta^{\prime})=q(\theta|\theta^{\prime})f(\theta) \end{aligned}$

と補正する。
　これを基にMetropolis-Hastings法のアルゴリズムを定式化すると「提案された候補点 $a$ を確率 $\min(1,r)$ で受容( $\theta^{(t+1)}=a$ )するか、その場に留まる( $\theta^{(t+1)}=\theta^{(t)}$ )ことを繰り返す」：

提案分布 $q(\cdot|\theta^{(t)})$ を利用して乱数 $a$ を生成する。
生成した $a$ を用いて
$\begin{aligned}q(a|\theta^{(t)})f(\theta^{(t)})\gt q(\theta^{(t)}|a)f(a)\end{aligned}$
ならば、 $\theta=a$ 、 $\theta^{\prime}=\theta^{(t)}$ と見なし、
$\begin{aligned}r=\displaystyle{\frac{q(\theta^{(t)}|a)f(a)}{q(a|\theta^{(t)})f(\theta^{(t)})}}\end{aligned}$
を計算し、確率 $r$ で $a$ を受容し、 $\theta^{(t+1)}=a$ とする。確率 $1-r$ で $a$ を破棄し、 $\theta^{(t+1)}=\theta^{(t)}$ とする。
$\begin{aligned}q(a|\theta^{(t)})f(\theta^{(t)})\leq q(\theta^{(t)}|a)f(a)\end{aligned}$
ならば、確率 $1$ で $a$ を受容し $\theta^{(t+1)}=a$ とする。
$t=t+1$ として1.に戻る。