「大人の教養・知識・気付き」を伸ばすブログ

一流の大人（ビジネスマン、政治家、リーダー…）として知っておきたい、教養・社会動向を意外なところから取り上げ学ぶことで“気付く力”を伸ばすブログです。

MENU

【ベイズ統計学】入門　ベイズ統計学（06/X）

統計学

　Bayes統計学を今回は

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

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

Amazon

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

4.　Metropolis-Hastings法

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

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

  1. 提案分布q(\cdot|\theta^{(t)})を利用して乱数aを生成する。
  2. 生成した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とする。
  3. 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法のアルゴリズム

  1. 提案分布q(\cdot|\theta^{(t)})を利用して乱数aを生成する。
  2. 生成した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とする。
  3. 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法のアルゴリズム

  1. 提案分布q(\cdot|\theta^{(t)})を利用して乱数aを生成する。
  2. 生成した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とする。
  3. 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}

が成り立つから

\begin{aligned} f(x^{*}|\boldsymbol{x})\approx \displaystyle{\frac{1}{T-B}\sum_{t=B+1}^{T}f(x^{*}|\theta^{(t)})} \end{aligned}

*1:厳密にはサンプルの最初からの番号である。

プライバシーポリシー お問い合わせ