こんにちは，株式会社Nospare・千葉大学の小林です．今回はMartin, Frazier and Robert (2020)によるベイズ統計における数値計算（Bayesian computation）の進歩の歴史についてまとめた論文`Computing Bayes: Bayesian computation from 1763 to the 21st Century（arXivのリンク）'について紹介します．本記事は前編として，論文の前半部分にあたる1763年から20世紀までの発展について紹介します．
この論文の狙いはベイズ統計にあまり馴染みのない人に向けて，
- なぜBayesian computationにはたくさんの方法があるのか？
- それらのつながりは？
- いつどの手法を使うのか？
などといった問いに対してインサイトを与えるところにあります．そのために共通の記法を用いて歴史に沿ってBayesian computationについて紹介していきます．
Bayesian computationの方法は大まかに次の３通りありますが，
- 決定論的な積分法
- シミュレーションによる方法
- 近似（漸近近似含む）による方法
この論文は主に2と3に焦点を当てています．
この論文ではパラメトリックモデル（有限個の未知数）に対するシミュレーションや近似に基づく方法のみを取り扱い，ノンパラメトリックモデルは扱いません．また，ベイズ最適化，ベイジアンブートストラップ，ベイジアン経験尤度などもスコープ外としています．一方で最新の理論の展開などについてもたくさん言及しているので参考文献をチェックするのにも役立つと思います．
はじまり
ベイズは「成功確率$\theta$のベルヌーイ試行を$n$回行ったときの$\theta$が$a$と$b$に含まれる確率は何か」というとても単純なモデルを考えていました．これを現在の記法で表記すると$Y_i|\theta\sim Bernoulli(\theta),\ i=1,\dots,n$で尤度関数$p(y|\theta), y=(y_1,\dots,y_n)$，$(0,1)$上の一様事前分布$p(\theta)$と書けます．ベイズは
$$
P(a<\theta<b)|y)=\int_a^b p(\theta|y)d\theta
$$
を計算しようとしていました．ここで$p(\theta|y)$は事後密度
$$
p(\theta|y)=\frac{p(y|\theta)p(\theta)}{p(y)}
$$
で，$p(y)=\int_0^1p(y|\theta)p(\theta)d\theta$は周辺尤度です．$x=\sum_{i=1}^ny_i$とすると，事後分布はベータ分布
$$
p(\theta|y)=\frac{\theta^x(1-\theta)^{n-x}}{B(x+1,n-x+1)}
$$
となります（$B(x+1,n-x+1)=\Gamma(x+1)\Gamma(n-x+1)/\Gamma(n+2)$はベータ関数，$\Gamma(x)$はガンマ関数）．ベイズは上記の事後確率を計算しようとしていたため，不完全ベータ関数の評価が必要でしたが当時は困難でした．こうして事後分布に関する数量を計算する必要が認識されるようになりました．
Bayesian computation
数値計算の問題
ベイズが計算しようとしていた確率は事後期待値$E[I_{[a,b]}|y]=\int_\Theta p(\theta|y)d\theta$の形で書けます（$I_{[a,b]}$は指示関数）．
一般的に，$\theta=(\theta_1,\dots,\theta_p)'\in\Theta$とし，結合事後密度を$p(\theta|y)$とすると，ベイズ推定ではある関数$g(\theta)$の事後期待値に関心があります：
$$
E[g(\theta)|y]=\int_\Theta g(\theta)p(\theta|y)d\theta
$$
事後確率に加えて，事後モーメント$E[\theta|y]=\int_\Theta\theta p(\theta|y)d\theta$， $Var(\theta|y)=\int_\Theta(\theta-E[\theta|y])(\theta-E[\theta|y])' p(\theta|y)d\theta$や，$g(\theta)=p(y_{n+1}^*|\theta,y)$とすると$y_{n+1}$の予測分布，意思決定$d$についての損失関数$L(\theta,d)$について$g(\theta)=L(\theta,d)$とするとベイズ決定論の最小値も含まれます．またモデル$\mathcal{M}$の周辺尤度は$g(\theta)=p(y|\theta,\mathcal{M})$として次の事前分布$p(\theta|\mathcal{M})$に対する期待値で書けます：
$$
E[g(\theta)|\mathcal{M}]=\int_\Theta g(\theta)p(\theta|\mathcal{M})d\theta
$$
周辺尤度の比はベイズファクターとなり，2つのモデルの比較に用いられます．何れにせよ，ベイズ分析（推論，予測，決定論，モデル選択）において必要となる数量は期待値で表されることになります．
早い時期の計算方法
1763年：ベイズの積分
上のベルヌーイ分布の例において，ベイズは$x$あるいは$n-x$が小さい場合に二項展開と項ごとの積分によって解析的に上述の事後確率を求めることができました．しかし，$x$と$n-x$が両方とも大きな場合にはこの方法では計算ができないため，求積法による上限と下限を求めるにとどまりました．
1774年：ラプラスの漸近近似
ラプラスは最初ベルヌーイモデルにおいて，任意の$w$について$P(|x/n-\theta|<w|y)=P(x/n-w<\theta<x/n+w|y)\rightarrow 1$を示そうと試みていました（事後一致性）．しかしベータ事後分布について$P(a<\theta<b|y)$の計算というベイズと同じ問題に躓きました．結局ラプラスは適当な条件のもとでベータ事後分布を正規分布で大標本による漸近近似を用いました：
$$
P(a<\theta<b|y)=\int_a^b p(\theta|y)d\theta=\int_a^b \exp(n f(\theta))d\theta
$$
ここで$f(\theta)=n^{-1}\log p(\theta|y)$．そして$f(\theta)$のモード$\hat{\theta}$周辺で２階のテイラー展開$f(\theta)\approx f(\hat{\theta})+\frac{1}{2}f^{''}(\hat{\theta})(\theta-\hat{\theta})^2$を用いて正規近似を行いました：
P(a<\theta<b|y)\approx\exp(nf(\hat{\theta}))\int_a^b \exp\left(-\frac{1}{2\sigma^2}(\theta-\hat{\theta})^2\right)d\theta\\
=\exp(nf(\hat{\theta}))\sqrt{2\pi\sigma^2}\left(\Phi\left(\frac{b-\hat{\theta}}{\sigma}\right)-\Phi\left(\frac{a-\hat{\theta}}{\sigma}\right)\right)
ここで$\sigma^2=-[nf^{''}(\hat{\theta})]^{-1}$，$\Phi(\cdot)$は標準正規分布の累積確率関数です．ラプラスの考え方は一般的な事後期待値の近似に用いられるようになりました．
ラプラスの後170年の間，積分を計算するための方法に大きな進歩は有りませんでした．それには$p(\theta|y)$から乱数を発生させるという新しい考え方とそれが可能になるためのマシンが必要でした．
1953年：モンテカルロ・シミュレーションとメトロポリスアルゴリズム
メトロポリスらは
$$
E[g(x)]=\int_\mathcal{X}g(x)p(x)dx
$$
という期待値を計算しようとしていました．ここで$p(x)$は$N$この粒子の$R^2$上の集合のボルツマン分布です．ここでの問題は次のようなものでした．
- 積分の次元が$2N$と大きい
- $p(x)$は通常正規化定数が未知である
1.の問題は$L$個のグリッドを各$2N$次元にわたる矩形積分の適用が不可能であり（近似誤差が$O(L^{-1/2N})$）であり，2.の問題は$M$個の$p(x)$からのiidな乱数$x^{(i)}, i=1,\dots,M$に基づく次元によらない近似誤差$O(M^{-1/2})$を持つようなモンテカルロ推定値$\hat{E}_{MC}(g(x))=\sum_i g(x^{(i)})$を得ることができないことを意味しています．
以下ではこれまでのベイズ推定の表記に戻って書きます．彼らは不変分布$p(\theta|y)$を持つマルコフ連鎖$\theta^{(i)}, i=1,\dots,M$をシミュレートすることによって事後期待値を計算することを提案しました．$i+1$回目のイタレーションでの値は$i$回目の値$\theta^{(i)}$をもとに候補の値$\theta^c=\theta^{(i)}+\delta\epsilon$を作り（$\epsilon\sim U(-1,1)$，$\delta$はチューニングパラメータ），確率
\alpha=\min\left\{\frac{p^*(\theta^c|y)}{p^*(\theta^{(i)}|y)},1\right\}
で$\theta^c$を$\theta^{(i+1)}$として採択します（$p^*$は$p$の核）．可逆的なマルコフ連鎖の理論に基づいて，このアルゴリズムは不変分布$p(\theta|y)$からの独立でないサンプルを生成します．これらのサンプルの標本平均によって事後期待値を
\bar{g(\theta)}=\frac{1}{M}\sum_{i=1}^Mg(\theta^{(i)})
によって推定します．この推定量について大数の弱法則や中心極限定理によって$\sqrt{M}$一致性や漸近正規性を示すことができます．
マルコフ連鎖の自己相関によって，メトロポリス推定量は（実行不可能な）$p(\theta|y)$からのiidな乱数に基づいたモンテカルロ推定量よりも分散が大きくなりますが，メトロポリスMCMCアルゴリズムは$p(\theta|y)$の正規化定数がわからなくても使用でき，また$p(\theta|y)$自体からから直接シミュレーションを行わなくてもよいという利点によりこのあと大きく発展していくことになります．
1964年：重点サンプリング
重点サンプリング（importance sampling, IS）は当初は積分をシミュレーションに基づいて推定する推定量の分散を小さくする目的で考案されました．重点サンプリングについては米倉さんの記事でも解説が行われていますので参照ください．
$p(\theta|y)$をと似ている提案密度$q(\theta|y)$が与えられたとき，IS推定量は$q(\theta|y)$からのiidなサンプルをもとに$\bar{g(\theta)}_{IS}=M^{-1}\sum_i g(\theta^{(i)})w(\theta^{(i)})$の形で書けます（ここで重み$w(\theta^{(i)})=p(\theta^{(i)}|y)/q(\theta^{(i)}|y)$）．特に$p(\theta^{(i)}|y)$の正規化定数がわからないときには$w(\theta^{(i)})$を評価することができないので，推定量は次のように修正されます：
\bar{g(\theta)}_{IS}=\frac{\sum_{i=1}^M g(\theta^{(i)})w(\theta^{(i)})}{\sum_{j=1}^Mw(\theta^{(j)})}
ここで$w(\theta^{(j)})=p^* (\theta^{(i)}|y)/ q^* (\theta^{(i)}|y)$で
$p*$と$q*$はそれぞれ$ p^* $と$q^*$の核です．重点サンプリングの利点はやはりメトロポリスアルゴリズムと同様に$p(\theta|y)$から直接シミュレートすることができないときに利用できる，というところにあります．さらに独立なサンプルを用いるという点は，重みのリウェイティングと組み合わせてパーティクルフィルターといった逐次モンテカルロ法の発展につながっていき，パーティクルMCMCなどに用いられるようになります．
1970年：ヘイスティングスによるメトロポリスアルゴリズムの一般化
ヘイスティングスらは$p(\theta|y)$の核だけ評価することができればよいというメトロポリスアルゴリズムの利点を認識しており，メトロポリスアルゴリズムの採択確率を次のように一般化しました：
\alpha=\min\left\{\frac{p^*(\theta^c|y)q(\theta^{(i)}|\theta^c,y)}{p^*(\theta^{(i)}|y)q(\theta^c|\theta^{(i)},y)},1\right\}
この採択確率は$q(\theta|y)$が$\theta^c$と$\theta^{(i)}$について大将であるときもとのランダムウォークメトロポリスアルゴリズムになります．この一般化によって，チューニングの手間を減らし潜在的には自己相関つまり$E[g(\theta)|y]$の推定値の分散を減らすように$q(\theta|y)$を選択することができるようになりました．$\theta$の次元が大きいときには$\theta$をブロックにわけて他のブロックを条件づけてあるブロックを更新することを提案しました．
20世紀終盤：ギブスサンプリングとMCMC革命
20世紀終盤には
- コンピュータの速度アップ，特にデスクトップのパソコン
- 結合事後分布$p(\theta|y)$からのMCMCのシミュレーションは，低次元でしばしば標準的な条件付事後分布からの反復的なサンプリングによって構成することができる
ということから，シミュレーションに基づいたBayesian computationはめざましく発展しました．事後分布の拡大（augmentation）とMCMCアルゴリズムの理論的な解明によって2.はギブスサンプリング（メトロポリスヘイスティングス（MH）副連鎖を含む場合もある）とつながり，1990年代のBayesian computationの中心的な役割を担っていくことになりました．
MHアルゴリズムは上記の採択確率によって$p(\theta|y)$に収束するようなマルコフ連鎖が構築されることにより機能します．もっと言えば，提案密度$q(\theta|y)$と採択確率は不変分布$p(\theta|y)$を持つ遷移核を構成することになります．ギブスサンプラーも同様に$p(\theta|y)$が不変分布となるようなマルコフ連鎖を構築しますが，遷移核は完全条件付事後分布の積で定義されます．例えば単純な二次元ベクトル$\theta=(\theta_1,\theta_2)'$の場合にはギブスサンプラーは，まず$\theta_2$の初期値$\theta^{(0)}_2$を設定し，$i=1,\dots,M$について$\theta_1^{(i)}$を$p_1(\theta^{(i)}_1|\theta^{(i-1)}_2,y)$からサンプリングし，$\theta_2^{(i)}$を$p_1(\theta^{(i)}_2|\theta^{(i)}_1,y)$からサンプリングするというステップを繰り返します．
ギブスサンプリングの考え方は，未知数$\theta$を潜在データ$z=(z_1,\dots,z_n)'$に拡大して条件付分布$p(\theta|z,y)$と$p(z|\theta,y)$を得るという方法と関連していました．これにより潜在データの導入により$p(\theta|z,y)$が標準的な形にすることができシミュレーションを行いやすくなり，拡大された空間上で扱いやすい条件付事後分布にすることで，高次元の複雑なモデルをベイズの枠組みで取り扱えるようになりました．
後編に続く
本記事では論文の前半部分にあたる20世紀までのBayesian computationの発展の歴史について紹介しました．後編では論文の後半部分にあたる21世紀に入ってからのさらなるMCMCの革命と最新の動向について紹介します．
データサイエンス研修やります！
株式会社Nospareではベイズ統計学に限らず統計学の様々な分野を専門とする研究者が所属しており，新たな知見を日々追求しています．統計アドバイザリー，ビジネスデータの分析，データサイエンス研修につきましては弊社までお問い合わせください．インターンや正社員も随時募集しています！
