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

　Bayes統計学を今回は

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

作者:豊田秀樹
朝倉書店

Amazon

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

4.　Metropolis-Hastings法

4.　Metropolis-Hastings法

　復習としてMetropolis-Hastings法のアルゴリズムを記述する：

　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.に戻る。

以降、特殊なケースにおけるMetropolis-Hastings法(MH法)を取り扱う。

4.7　独立MH法

　提案分布は一般的に条件付き分布の形式で表現する。これを敢えて1時点前の条件付きでない無条件分布を提案分布に用いることもできる。この場合、提案候補は互いに独立になる。このようなMH法を独立MH法という。
　独立MH法では、提案分布 $q(\cdot)$ が無条件分布であるから、

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

と表現できる。すなわち

図表1　独立Metropolis-Hastings法のアルゴリズム

提案分布 $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)})f(a)}{q(a)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.に戻る。

　独立MH法は提案分布の良し悪しに応じて、収束までの成績に大きな違いが生じる。実際のデータ分析では目標分布、更には目標分布の母数ですら不明であるため、不適切な提案分布を選ぶリスクがある。

4.7　ランダムウォークMH法

　独立MH法における「不適切な提案分布を選ぶリスク」を解決するには、ランダムウォークを利用することが効果的である。具体的には候補提案を

$\begin{aligned} a=\theta^{(t)}+\varepsilon,\ \varepsilon\sim F \end{aligned}$

と攪乱項を持たせる。分布 $F$ としては、平均が $0$ であるような正規分布や区間 $[-\delta,\delta],\ \delta\gt0$ の一様分布などの対称な分布を選ぶ。
　提案分布は、対称な分布を選んだ方が便利である。対称な分布を選ぶと提案分布が

$\begin{aligned} q(a|\theta^{(t)})=q(\theta^{(t)}|a) \end{aligned}$

となる。このときランダムウォークMH法における補正係数は

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

と書ける。

図表2　ランダムウォークMetropolis-Hastings法のアルゴリズム

提案分布 $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{f(a)}{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.に戻る。

4.8　生成量・仮説が正しい確率

　一般に期待値は大数の法則により、 $T$ が十分に大きければ

$\begin{aligned} {\hat{\theta}}_{EAP}\approx\displaystyle{\frac{1}{T-B}\sum_{t=B+1}^{T}\theta^{(t)}} \end{aligned}$

が成り立つことがMarkov連鎖Monte Carlo法の動機であった。ここで $B$ はバーンイン期間として控除するサンプル数*1である。
　事後分散は $T$ が十分に大きいので不偏性を持たせるための補正を無視すれば

$\begin{aligned} \hat{V}[\Theta]=\displaystyle{\frac{1}{T-B}\sum_{t=B+1}^{T}\left(\theta^{(t)}-{\hat{\theta}}_{EAP}\right)^2} \end{aligned}$

である。
　一般の関数 $g(X)$ に関する期待値

$\begin{aligned} E[g(X)]=\displaystyle{\int g(X)f(x|\theta) dx} \end{aligned}$

でも、母数とデータを入れ替えた

$\begin{aligned} E[g(\Theta)]=\displaystyle{\int g(\Theta)f(\theta|\boldsymbol{x}) d\theta} \end{aligned}$

について大数の法則から

$\begin{aligned} E[g(\Theta)]\approx \displaystyle{\frac{1}{T-B} \sum_{t=B+1}^{T}g(\theta^{(t)})} \end{aligned}$

と近似できる。このとき $g(\theta^{(t)})$ を生成量という。MCMC法では生成量を大量に与えるため、その標準誤差は事後標準偏差の推定値として利用できるため、それからパーセント値を利用すれば生成量の確信区間を数値的に求めることが出来る。
　特に $g(x)$ として定義関数 $\boldsymbol{1}_{\{\boldsymbol{x}\in U\}}(\boldsymbol{x})$ ( $U$ は確率を求めたい仮説の事象)を用いることで

$\begin{aligned} u^{(t)}=g(\theta^{(t)})=\boldsymbol{1}_{\{\boldsymbol{x}\in U\}}(\theta^{(t)})=\begin{cases} 1,&\ \ \ \theta^{(t)}\in U\\ 0,&\ \ \ \theta^{(t)}\notin U \end{cases} \end{aligned}$

であり、 $\bar{u}^{(t)}=\displaystyle{\frac{1}{T-B}\sum_{t=B+1}^{T}\boldsymbol{1}_{\{\theta^{(t)}\in U\}}(\theta^{(t)})}$ で仮説 $U$ が正しい確率を評価することが出来る。

4.9　事後予測分布の評価方法

　事後予測分布について

$\begin{aligned} f(x^{*}|\boldsymbol{x})&=\displaystyle{\int_{-\infty}^{\infty}f(x^{*}|\theta)f(\theta|\boldsymbol{x})d\theta }\\ &=E_{f(x^{*}|\boldsymbol{x})}\left[f(x^{*}|\theta)\right] \end{aligned}$