こんにちは，株式会社Nospare・千葉大学の小林です．今回は前回の記事「【論文紹介】ベイズ統計における数値計算の進歩の歴史（前編：1763年から20世紀まで）」に引き続き，Martin, Frazier and Robert (2020)によるベイズ統計における数値計算（Bayesian computation）の進歩の歴史についてまとめた論文`Computing Bayes: Bayesian computation from 1763 to the 21st Century（arXivのリンク）'について紹介します．本記事は中編として，論文の中盤部分にあたる21世紀における発展について紹介します．

21世紀：２回目の革命

なぜ２回目が必要だったか？

20世紀終盤でMCMCによるBayesian computationの手法が台頭し，非常に多くの分野・領域で利用されるようになりました．しかし応用研究における必要性から，これらの初期のMCMCアルゴリズムには制約があることが認識され始めました．

多くの問題において$p(\theta|y)\propto p(y|\theta)p(\theta)$は事後密度の核$p^*(\theta|y)$までしかわからず，直接シミュレーションを行うことができない．
事後密度の核がわかっていても，$p(y|\theta)$と$p(\theta)$は解析的に得られている必要があり，特に$p(y|\theta)$は観測値$y$において評価できる必要がある．

MCMCやISはある意味間接的に$p(\theta|y)$からシミュレーションを行うことで最初の問題を解決しましたが，それでも$p(y|\theta)$の評価できることが必要とされます．よって，この仮定は

データ生成過程が解析的に確率密度を伴わない場合（未知の遷移密度をもつような連続時間モデルや正規化定数が未知であるような空間モデルなど）
$p(y|\theta)$の評価（可能であっても）が$O(n)$の計算コストがかかり，ビッグデータの問題に対してスケーラブルではない場合

に制約的になってしまいます．

また，積分が必要な潜在変数の数が多い場合など，未知変数の次元がとても大きい（高次元）とき，通常のMCMCは適用できたとしても有限の計算時間で十分に精確な推定値を得ることが難しいことがあるため，やはり未知変数の次元に対して必ずしもスケーラブルではないと言えます．以下で紹介する方法は，何らかの近似を導入することによってこれらの課題のいくつかを克服しようという試みとなっています．

Pseudo-marginal（PM）法

データ拡大法では，パラメータ$\theta$と潜在変数$z$の同時事後分布$p(\theta,z|y)$からのシミュレーションを（ギブスサンプラーなどで）行うことで，$p(\theta|y)$からのシミュレーションを得ることができるというものでした．しかし，問題によっては$(\theta,z)$の空間をMCMCが遷移するのが遅いという欠点もあります．

$z$を生成するときに使用する確率変数を$u\in \mathcal{U}$とします．PMのアイディアは，もし$u$を使うことで，尤度関数$p(y|\theta)$の不変推定量を得ることができるのであれば，$(\theta,u)$の空間でのMCMCは不変分布$p(\theta|y)$からのシミュレーションが行える，というものです．$h(u)$を$z$を生成するのに用いられる確率分布とし，$h(y|\theta,u)$を$p(y|\theta)$の不変推定量とします：$E_u[h(y|\theta,u)]=p(y|\theta)$．このとき

h(\theta|y)\propto \int_\mathcal{U}h(y|\theta,u)p(\theta)h(u)du = p(\theta)E_u[h(y|\theta,u)] = p(\theta)p(y|\theta)\propto p(\theta|y)

となります．MHアルゴリズムにおいて$h(y|\theta,u)$を使うと採択確率は

\min\left\{\frac{h(y|\theta^c,u^c)p(\theta^c)q(\theta^{(i)}|\theta^c,y)}{h(y|\theta^{(i)},u^{(i)})p(\theta^{(i)})q(\theta^{c}|\theta^{(i)},y)},1\right\}

と書け，ここで$\theta^c$は$q(\theta^c|\theta^{(i)},y)$からの提案で，$u^{(i)}$と$u^{c}$は$h(u)$から独立に生成されます．この尤度の不変推定量によって置き換えられた採択確率をもつMHアルゴリズムをPMMHアルゴリズムと呼びます．

このPMの考え方は状態空間モデルに対するパラメータ推定にも利用でき，パーティクルフィルターによって尤度関数の不変推定値を得て，MHアルゴリズムによってパラメータの事後分布からシミュレーションを行うアルゴリズムをパーティクルMCMCと呼んだりします（「[Python]particlesパッケージで逐次モンテカルロ法を適用してみた（中編）」では，pythonでの適用を扱ってます）．

PMの利点は3つ挙げられます：

ギブスサンプラーではパラメータと潜在変数の相関が高い（$p(\theta|z,y)$と$p(z|\theta,y)$の）のに対し，PMはもっと効率的であると期待できる
潜在変数に対して前向きのシミュレーションしかできない場合には$p(y,z|\theta)$が評価できないのに対して，その場合においてもPMMHは適用可能（例えばブートストラップフィルターでは$p(z_t|z_{t-1},\theta）$からのシミュレーションは可能）
$y$の次元が大きいときには，$p(y|\theta)$の不変推定値はデータのサブセットからでも構築することができ，PMMHの計算負荷のコントロールができる

この方法での事後期待値の推定値はexactとなります（シミュレーション誤差はあります）．

近似ベイズ推論（approximate Bayesian inference）

シミュレーションに基づくBayesian computationにおいても近似を用いることで，上述の課題を克服しようという試みがなされています．以下で紹介するapproximate Bayesian computation（ABC）とBayesian synthetic likelihood（BSL）は$p(y|\theta)$を直接評価することを避けることができ，変分ベイズ（variational Bayes, VB）とintegrated nested Laplace approximation（INLA）は$p(y|\theta)$の評価は必要であるものの，シミュレーションをある種の最適化問題に置き換えることで計算負荷の軽減することが可能です．

ABC

ABCは，$p(y|\theta)$を直接評価することが困難であっても，$p(y|\theta)$（確率モデル）からのデータのシミュレーションを行える場合に$p(\theta|y)$を近似するために適用することができます．最も単純なABCでの棄却法は次のように行います

$p(\theta)$から$\theta^i$を生成し，それらのもとで人工でーた$x^i$を$p(\cdot|\theta^i)$から生成する（$i=1,\dots,M$）
$x^i$をもとにシミュレーションの要約統計量$\eta(x^i)$を作り，観測された要約統計量$\eta(y)$を作る
距離$d(\cdot,\cdot )$と小さな閾値$\epsilon$を使い，$d(\eta(x^i),\eta(y))\leq \epsilon$となるような$\theta^i$を残す．

このように，ABCは（十分統計量ではない）統計量$\eta(y)$を所与とした事後分布$p(\theta|\eta(y))\propto p(\eta(y)|\theta)p(\theta)$をターゲットとして$p(\theta|y)$の近似とします．上記の棄却法における尤度関数は，

p_{\epsilon}(\eta(y)|\theta)=\int_\mathcal{X}p(x|\theta)I(d(\eta(y),\eta(x))\leq \epsilon)dx

と書け，$\theta^i$と$\eta(x^i)$が与えられたとき，これのシミュレーションでの対応は

\hat{p}_\epsilon(\eta(y)|\theta^i)=I(d(\eta(y),\eta(x^i))\leq \epsilon)

と書けるため，一様カーネルを用いたノンパラメトリック推定量と見ることができます．
またこれは，$p_\epsilon(\theta|\eta(y))\propto p_\epsilon(\eta(y)|\theta)p(\theta)$に対するPMにおける尤度関数の推定量としても見ることができます．

BSL

BSLも$\eta(y)$を所与とした$\theta$の事後分布をターゲットとし，$p(y|\theta)$から人工データを生成可能であれば適用することができます．ABCが暗に尤度関数のノンパラメトリックな推定値を使うのに対して，BSLは$p(\eta(y)|\theta)$に対して正規分布によるパラメトリックな近似を行います：

p_a(\eta(y)|\theta)=\mathcal{N}(\eta(y);\mu(\theta),\Sigma(\theta))

ここで$\mu(\theta)=E[\eta(y)]$，$\Sigma(\theta)=Var[\eta(y)]$．パラメトリックなカーネルによって，（理想的な）BSL事後分布は

p_a(\theta|\eta(y))\propto p_a(\eta(y)|\theta)p(\theta)

となります．ここで，$a$は$p(\eta(y)|\theta)$に対する$p_a(\eta(y)|\theta)$による正規近似を表します．一般的には$\eta(y)$の平均と共分散行列は未知であるので，これをシミュレーションに基づいて次のように推定します．$x_j\sim i.i.d. p(\cdot|\theta), \ j=1,\dots,m$を生成し，$\mu(\theta)$と$\Sigma(\theta)$をモンテカルロ近似$\mu_m(\theta)=\frac{1}{m}\sum_{j=1}^m\eta(x_j)$，$\Sigma_m(\theta)=\frac{1}{m-1}\sum_{j=1}^m(\eta(x_j)-\mu_m(\theta))(\eta(x_j)-\mu_m(\theta))'$に置き換えます．これによって尤度関数の近似を

p_{a,m}(\eta(y)|\theta)=\int_\mathcal{X}\mathcal{N}(\eta(y);\mu_m(\theta),\Sigma_m(\theta))\prod_{j=1}^mp(\eta(x_j)|\theta)dx_1,\dots,dx_m

とすると，ターゲットとなるBSL事後分布は

p_{a,m}(\theta|\eta(y))\propto p_{a,m}(\eta(y)|\theta)p(\theta)

となります．

VB

ABCやBSLは（十分統計量とは限らない）$\eta(y)$を所与とした事後分布の近似をターゲットとしているのに対し，VBはシミュレーションを最適化問題にお置き換えることで事後分布自体$p(\theta|y)$に対して直接近似を粉います．

VBでは，ある密度分布のクラス$\mathcal{Q}$を設定し，その中から$p(\theta|y)$をKLダイバージェンスの意味で最もよく近似できるようなものを見つけるということを行います：

q^*(\theta)=\arg\max_{q(\theta)\in\mathcal{Q}} KL [q(\theta)|p(\theta|y)]

ここで，

\begin{split}
KL[q(\theta)|p(\theta|y)]&=\int\log(q(\theta))q(\theta)d\theta)-\int\log(p(\theta|y))q(\theta)d\theta \\
&= E_q[\log q(\theta)]-E_q[\log(\theta,y)]+\log p(y)
\end{split}

正規化定数$p(y)$は通常未知ですが，VBでは不要です．VBではevidence lower bound（ELBO）取り扱います．

ELBO[q(\theta)]=E_q[\log(\theta,y)]-E_q[\log q(\theta)]

$KL[q(\theta)|p(\theta|y)]$は定数項$\log p(y)$以外まで$-ELBO[q(\theta)]$と等しいので，次の最大化問題を考えることができます：

q^*(\theta)=\arg\max_{q(\theta)\in\mathcal{Q}} ELBO[q(\theta)]

VBの長所としては，特定の問題で特定のクラス$\mathcal{Q}$が与えられたとき（特に標準的な分布で定義されたとき）にはこの最大化問題は解析的に得ることができる，あるいは数値計算によって高速に解くことができる点が挙げられます．また，VBは特に$\theta$が高次元でMCMCを効率的に適用できないときに有用です．そのような場合には$\mathcal{Q}$を，事後分布の近似が$\theta$が高次元でも扱いやすいように選択することができます．そしてVBの副産物として，最大化問題を解くことで周辺尤度$p(y)$の下限も得ることができ，$ELBO[q^*(\theta)]$をモデル選択において使用する数量の推定値として用いることが可能です．

INLA

1774年，ラプラスは事後一致性を示すために，特定の事後確率の$n$についての漸近近似を求めました．1986年にTierney and Kadaneは，任意の事後期待値を漸近近似するためにラプラス近似を取り上げ理論化しました．さらに20年後Rue et al. (2009)はこの手法をさらに発展させ，潜在ガウスモデル（LGM）における周辺事後分布の近似に適用しました．一連の入れ子型のラプラス近似と低次元の数値積分を用いたこの手法はINLAと呼ばれます．LGMは一般化線形モデル，非ガウス型状態空間モデル，空間モデル，時空間モデルなど非常に幅広いモデルが含まれているため，INLAの適用範囲はとても広くなります．INLAの具体的な方法についてはNospareの過去の記事「Integrated Nested Laplace Approximationを用いた近似ベイズ解析」で解説されてますのでそちらを参照ください．

ハイブリッドな方法

ここまでで4つのアイディアを取り上げてきました．それぞれの長所・短所をまとめると，

ABCとBSLは$p(y|\theta)$が評価できないときに適用できるが，$\theta$（また$\eta(y)$）が高次元の場合に適用が困難になる
VBは$\theta$が高次元のときに有用であるが，$p(\theta,y)$つまり$p(y|\theta)$の評価が必要になる
PMは$p(y|\theta)$の評価が困難であるときに，それを不変推定量$h(y|\theta,u)$で置き換えることができる

というものでした．最近ではこれら4つのアルゴリズムの長所を組み合わせたハイブリッドな手法（例えば，尤度関数が解析的に得られなくてかつ$\theta$が高次元な場合など）も考えられてきています．

MCMCアルゴリズム再訪

新しい（近似的なものを含む）計算手法が次々と登場しているにもかかわらず，20世紀後半の主役であるMCMCが別の手法にとって変わられたとはいえません．ここでは，最初のアルゴリズムが登場した後に行われたMCMCの主な進歩について簡単に紹介します．これらの進歩の多くは，大規模なデータセットや高次元の未知数を効率的に扱う必要性に駆られたものです．

まず，MCMCアルゴリズムに関する3つの注意点から始めます．

MC-MCアルゴリズムとはまさにマルコフ連鎖モンテカルロアルゴリズムのことで，アルゴリズムの設計上、ターゲットの事後分布の局所的な探索を行います．シミュレーションされた値のパラメータ空間における位置は，前の値の位置に依存します．問題としては，自己相関のの高いMCMCアルゴリズムは，ターゲット事後分布の事後密度が高い領域の探索に時間がかかる可能性があり，この問題は通常パラメータ空間の次元が大きいほど深刻になります．別の味方をすると，$M$個のMCMCシミュレーションについて，それらの（通常は正の）自己相関の度合いが大きいほど，ターゲットからの$M$個のiidのシミュレーションに基づく推定値と比較して，MCMCベースの推定値の（分散が大きくなる意味で）効率性が低下します．
MC-MCアルゴリズムは、マルコフ連鎖モンテカルロアルゴリズムでもあります．すなわち，適切な仮定の下では，連鎖の自己相関の程度に関わらず，事後期待値の$\sqrt{M}$-consistentな推定値を生成します．原理的には，どのようなMCMCアルゴリズムでも$M$を増やすだけで，任意の精度をもつ推定値を得ることができます．しかし$M$を増やすことは計算コストの増加を伴います．この増加の度合いは，提案を行うための（反復ごとの）コストと，MHステップでは採択確率を計算するためのコストに依存します．この2つのコストは各反復においてシミュレーションや評価が必要となる未知数の数に応じてに増加します．
効率性という概念は，マルコフ連鎖が（$M$について漸近的に）不偏である場合にのみ意味があり，これは正しい不変分布からシミュレーションが生成されることに依存します．つまり，精確なMCMCベースの推定値を生成するためには，連鎖の自己相関を減らすことや反復回数を増やすことだけではなく，連鎖が実際にターゲット分布を探索し，それによってバイアスを回避することを保証することも必要になります．

よってすべてのMCMCの進歩は基本的に，アルゴリズムがターゲットとしている事後分布の密度が高い領域を効率的に探索し，事後期待値の推定精度を高めることを目的としています．具体的には，

マルコフ連鎖の自己相関を減らす
反復あたりの計算コストを減らす
バイアスを取り除く

ということを行っています．それらに対するアプローチは

ターゲット事後分布の幾何学的な情報を用いる（HMC，no-U-turnサンプラーなど）
よりよい提案分布を使用する（適応的MCMC，テンパリング，multiple-try MHなど）
並列化，バッチ，サブサンプルを用いる
分散を減らす方法の適用（ラオ・ブラックウェル化など）

などに分類されます．

後もう少しだけ続きます

本記事では21世紀に入ってからのBayesian computationのさらなる発展について紹介しました．変分法は機械学習の文脈などで名前など聞いたことがあるかもしれません（MLPシリーズでも一冊でています）．ABCやBSLは2010年前後に一気に流行って，最初はMCMCやSMCなどのアルゴリズムの発展について研究が行われ，その後事後分布の漸近的性質や$\eta(y)$の選択についてなどの理論的な研究が行われるようになりました．
本当は元論文の最後の内容まで紹介したかったのですが，qiitaで数式を打つのに思ったより時間がかかったのでまた次回後編として記事を書きたいと思います．

データサイエンス研修やります！
株式会社Nospareではベイズ統計学に限らず統計学の様々な分野を専門とする研究者が所属しており，新たな知見を日々追求しています．統計アドバイザリー，ビジネスデータの分析，データサイエンス研修につきましては弊社までお問い合わせください．インターンや正社員も随時募集しています！

【論文紹介】ベイズ統計における数値計算の進歩の歴史（中編：21世紀）