Goncalves, F.B, Latuszynski, K., Roberts, G.O. (2017) Barker’s algorithm for Bayesian inference with intractable likelihoods. Brazilian Journal of Probability and Statistics. 31(4). 732-745.
本来私なんかが読むような論文じゃないんだけど、ちょっと事情があって前半部分のみ目を通した。すいません、すいません…
えーと、MCMCではイテレーションのたびに移動を受理するかどうか確率的に決めるじゃないですか。受理の確率を求めるために、MHアルゴリズムの場合だと密度比とかいう式を使いますわね。でもあの式の代わりに別の式を使おうという話があって、そのひとつとしてBarkerのアルゴリズムというのがあるのだそうだ。この論文はその解説。
2023/04/18追記: arXivに上がってたドラフトを読んでたんだけど、公刊されているのをみつけたので、そちらを読み直した。
1. イントロダクション
尤度関数の点ごとの評価が不可能だったり、計算的に極端に大変だったりする問題のことを、intractableな尤度という。intractabilityは、データの欠損とか、モデルの複雑性とか、データが足りないとかのせいで起きる。
データセット\(y\)の下での、状態空間\(\mathcal{X}\)におけるパラメータ\(\theta\)の事後分布を\(\pi(\theta) = \pi(\theta|y)\)とする。\(\pi\)を目標とするMCMCを設計するという問題を考えよう。\(\theta\)から提案される遷移の密度を\(q(\theta, \cdot)\)とする。
標準的なMHアルゴリズムでは、\(\theta\)から\(\phi\)への移動を確率$$ \alpha_{MH} = 1 \land \frac{\pi(\phi) q(\phi, \theta)}{\pi(\theta) q(\theta, \phi)} $$ で受容する。しかし、その計算には少なくとも\(\pi(\phi)\)の評価が必要になる。
近年では、確率そのものを計算することなく上の受容をシミュレートしようとするretrospectiveシミュレーションという方法に関心が集まっている。こうした手法は、確率\(\pi(\phi)\)を持つであろう出来事からのシミュレーションを利用する。しかし、たとえ確率\(\pi(\theta), \pi(\phi)\)を持つ出来事のretrospectiveシミュレーションが可能でありまた効率的であるとしても、上の式で表された確率を持つ出来事からのシミュレーションという問題への解にはつながらない。上の式は\(\pi(\theta), \pi(\phi)\)の非線形な表現だからである。[なにいってんだかよくわからん]
この問題は、ベルヌーイ工場問題という古典的な確率計算問題のひとつである。つまり、我々には評価できない確率 \(p\)があり、しかしその確率を持つ出来事はシミュレートできるとき、確率\(f(p)\)を持つ出来事をどうやってシミュレートするか、という問題である。
ここでは、関数\(f\)が\(f(p_1, p_2) = 1 \land \frac{c_1 p_1}{c_2 p_2}\)という形式を持っている。多くのモデリング文脈では \(f(p_3) = 1 \land c_3 p_3\)という形式にも書き換えられるだろう。
近年ではベルヌーイ工場問題に対する効率的な解が提案されている。\(f(p_3) = 1 \land c_3 p_3\)という形式のときの解は一般には存在しないことが知られているが、幸いなことに、MHアルゴリズムの柔軟性のおかげで、この問題を回避することができる。
本論文では、intractableな尤度に対するMCMCをBarker(1965)の受容確率を使って実装するという枠組みについて紹介する。Barker受容確率は効率的ではないが、単純なベルヌーイ工場型のアルゴリズム(2-coinアルゴリズム)をつかって簡単にシミュレートできる。
2. MHアルゴリズムに対するBarkerの代替案
提案カーネル\(q\)に対して、マルコフ連鎖が定常分布\(\pi\)を持つように受容確率\(\alpha(\theta, \phi)\)を定めるとき、その解はたくさんある。そのうち MHアルゴリズムの受容確率 \(\alpha_{MH}(\theta, \phi)\)は、すべての可能な遷移 \(\theta \rightarrow \phi\)の受容確率を最大化する方法である。このことは、MHアルゴリズムの受容確率が、\(L^2(\pi)\)における関数の積分の推定における漸近的モンテカルロ分散を最小化する受容確率であるということを示している。ここで\(L^2(\pi)\)とは\(\pi\)に関してsquare integrableな関数のヒルベルト空間を指す。[ああなるほど、square integrableな関数のヒルベルト空間ね… あれ酒のつまみにあうよね…。 すいません、途中から全く理解できなくなり写経しました]
しかし、だいたい同じくらい良い性質を持つ受容確率はほかにもたくさんある。ここではBarker(1965)の提案に焦点をあてよう。$$ \alpha_B(\theta, \phi) = \frac{\pi(\phi) q(\phi, \theta)}{\pi(\theta) q(\theta, \phi) + \pi(\phi) q(\phi, \theta)} $$ という受容確率である。$$ \frac{\alpha_{MH}(\theta, \phi)}{2} \leq \alpha(\theta, \phi) \leq \alpha_{MH}(\theta, \phi)$$ という性質がある。つまりどちらも性能は似ている。
[Barkerのアルゴリズムという言葉は\(\frac{\pi(\phi)}{\pi(\theta) + \pi(\phi)}\)のことを指すことにして(Metroplisアルゴリズム\(\frac{\pi(\phi)}{\pi(\theta)}\)との対比で)、上の式はBarker-Hastingsアルゴリズムだということにすると話が簡単なのだが… 歴史的にはBarker(1965)がすでにこの式だったのかもしれない]
厳密にCLTと漸近分散を比べるとこういう性質があることが知られている。\(f \in L^2(\pi)\)とする。そのiidモンテカルロ誤差を \(\sigma^2_\pi = \mathrm{Var}_\pi(f)\)とする。いま、\(f\)とMH連鎖について平方根中心極限定理が成立するならば、つまりMH分散のCLT漸近誤差を\(\sigma^2_{MH}\)として$$ \frac{\sum_i^N f(\theta_i) – \pi(f)}{\sqrt{N}} \Rightarrow N(0, \sigma^2_{MH})$$ ならば、\(f\)とBarker連鎖についてもCLTが成り立ち、Baker連鎖のCLT漸近分散を\(\sigma_B\)として $$ \sigma^2_{MH} \leq \sigma^2_B \leq 2\sigma^2_{MH} + \sigma^2_\pi $$ が成り立つ。
[難しくて何言ってんだかわかんないけど、Barkerの受容確率を使ったMCMC推定の誤差は通常のMHアルゴリズムよりは大きくなる、しかしその上界もわかっている、という話なんでしょうね]
いっぽう、\(\pi(\theta)\)と比例する出来事から効率よくシミュレーションしたいというとき、以下の2-coinアルゴリズムによって、確率\(\alpha_B(\theta, \phi)\)を持つ出来事から効率よくシミュレーションすることができる。
- \(C_1 \sim Ber \left( \frac{c_1}{c_1 + c_2} \right) \)をドローする。
- もし\(C_1 = 1\)なら、\(C_2 \sim Ber(p_1)\)をドローし、もし1だったら1を出力し、0だったら最初に戻る。
- もし\(C_0 = 0\)なら、\(C_2 \sim Ber(p_2)\)をドローし、もし1だったら0を出力し、0だったら最初に戻る。
このアルゴリズムから1が出力される確率は\(\frac{c_1 p_1}{c_1 p_1 + c_2 p_2}\)、0が出力される確率は\(\frac{c_2 p_2}{c_1 p_1 + c_2 p_2}\)になっている。平均実行時間は\(\frac{c_1 + c_2}{c_1 p_1 + c_2 p_2}\)に比例する。
このアプローチは、擬似周辺MCMCというアプローチと似ているんだけど、[…後略]。
3. かんたんな例
[てめえ、ほんとにかんたんなんだろうな、違ってたら殺すぞ… という気持ちです]
事後分布を $$ \pi(\theta) = E_{\eta \sim h} [\pi(\theta | \eta)] $$ だとしましょう。つまり、\(\pi\)は条件付き密度\(\pi(\cdot | \eta)\)の混合であり、混合を表す測度\(h\)はわかっている、でもそれ以上のことはわからない、としましょう。
もうちょっと具体的に、\(\theta | \eta \sim \mathrm{Poisson}(\eta) \)としましょう。すると$$ \pi(\theta | \eta) = \frac{\exp(-\eta) \eta^\theta}{\theta!} $$ となります[これはポアソン分布の定義ね]。\(d(\theta) = \frac{\exp(-\theta) \theta^\theta}{\theta!}\)とすると、\(\pi(\theta | \eta) \leq d(\theta)\)です[\(\pi(\theta | \eta)\)は「所与の時間中に平均で\(\eta\)回生じるできごとが\(\theta\)回生じる確率」で、\(d(\theta)\)は「所与の時間中に平均で\(\theta\)回生じるできごとが\(\theta\)回生じる確率」だから、当然そうなる]。
以下では\(\pi(\theta) = d(\theta) \frac{\pi(\theta)}{d(\theta)} = d(\theta) p(\theta)\)と書くことにします。
\(p(\theta)\)のシミュレーションはかんたんです。まず \(\eta \sim h\)をドローします。つぎに標準一様分布から\(U\)をドローします。で、出来事 \( \{U \leq \pi(\theta | \eta) / d(\theta) \} \) に注目します。その確率は\(p(\theta)\)です。
では、Barkerアルゴリズムについて考えてみましょう。提案分布は対称的ランダムウォークから得られる、すなわち\(q(\theta, \phi) = q(\phi, \theta)\)としましょう。受容確率を $$ \alpha_B(\theta, \phi) = \frac{d(\phi) p(\phi)}{d(\phi) p(\phi) + d(\theta) p(\theta)} $$ とすれば、2-coinアルゴリズムを使うことができます。
実装例として、\(\eta \sim \mathrm{Gamma}(100, 5)\)の場合を示します。提案分布は\( \{\theta-10, \theta-9, \ldots, \theta+10\}\)上の一様分布にしました。受容率は0.367となりました。\(\pi(\theta)\)の厳密な分布は平均20, 分散24の負の二項分布になります。シミュレーションの結果を示します。良い近似が得られています。
4. Barker MCMCによるWright-Fisher族拡散の正確推論
[ここからはBarker受容関数によるMCMCの具体的なご利益を示す本論。パス]
——-
感想: 私には難しくて、読み返してもいまいち理解できない…