Sisson, S.A., Fan, Y. (2021) Likelihood-Free MCMC. in Brooks, S., et al.(eds) “Handbook of Markov Chain Monte Carlo,”, CRC Press.
尤度フリーMCMCについての解説。いま仕事の関係で悩んでいることについて、ちょっとした手がかりでも得られるかも、と思って目を通した。高い本を買っちゃったから無理にでもモトを取らなくちゃ、というのがもう一つの動機である。
1. イントロダクション
パラメータ\(\theta\)の事後分布\(\pi(\theta|y)\)がわからないどころか、尤度\(\pi(y|\theta)\)さえ評価できない。それでもなんとかして事後分布について調べたい。そんな辛い状況下でも、大丈夫、まだ生きる道はある。それが尤度フリーMCMCである。近似ベイジアン計算ともいう。
一番シンプルな例として、尤度フリー棄却サンプリングを紹介しよう。こういうアルゴリズムである。
- 事前分布\(\pi(\theta)\)から\(\theta’\)をドローする。
- モデル\(\pi(x|\theta’)\)からデータセット\(x\)を生成する。
- もし\(x \approx y\)ならば\(\theta’\)を受容する。
2. 尤度フリー理論・手法のレビュー
2.1 尤度フリーの基礎
[ここ、何度読んでも頭が混乱したので、思い余ってしばらく逐語訳]
困難な状況下でサンプラーの効率を高めるためのお決まりの方法は、目標事後分布を補足モデルのなかに埋め込むことである。この状況ではモデルに補足パラメータが導入されるが、その目的は単に計算を楽にすることにある。たとえばsimulated temperingとか、焼きなまし法がそうだ。
尤度フリー推論も似たアプローチに沿っている。目標事後分布\(\pi(\theta|y) \propto \pi(y|\theta) \pi(\theta)\)に補足パラメータ\(x\)を導入して$$ \pi_{LF}(\theta, x|y) \propto \pi(y|x, \theta)\pi(x|\theta)\pi(\theta)$$ とする。\(x\)は\(\pi(x|\theta)\)からドローされたデータセットで、\(y \in \mathcal{Y}\)と同じ空間上にある。後述するが、分布\(\pi(y|x, \theta)\)は、事後分布\(\pi(\theta|x)\)が、\(x\)と\(y\)が似ている領域では高い密度で重みづけられるように選んだ分布である。
[あ！ 逐語訳してやっと意味が分かった。さっきの尤度フリー棄却サンプリングでいうと、\(x \approx y\)のとき1、そうでないとき0になるような関数が\(\pi(y|x, \theta)\)に相当するわけね。さっきは\(\theta\)を使わずに値を決めたけど、使ってもいいわけだ]
密度\(\pi(y|x, \theta)\)は点\(x=y\)では\(\theta\)に関して一定だと仮定する。すると、\(c\)をなんらかの正の定数として\(\pi(y|y, \theta)=c\)となり、その結果、\(x = y\)においては目標事後分布が復元されることになる、つまり\(\pi_{LF}(\theta, y|y) \propto \pi(y|\theta) \pi(\theta)\)となる。
究極的な関心はふつうは周辺事後分布$$ \pi_{LF}(\theta|y) \propto \pi(\theta) \int_\mathcal{Y} \pi(y|x, \theta)\pi(x|\theta)dx$$にある。実践的には、この式の積分は、同時事後分布\(\pi_{LF}(\theta, x|y)\)を目標にしたなんらかのサンプラーのアウトプットから、補足データの実現値を単に捨てることによって、数値的になされる。いっぽう後述するように、\(\pi_{LF}(\theta| y)\)を直接に目標とするサンプラーを作るという路線もある。
2.2 事後分布近似の性質
\(\pi(y|x, \theta)\)が、\(x = y\)のときのみ密度を持ち、そうでなかったら0であるような密度関数だったら、\(\pi_{LF}(\theta|y)\)は目標事後分布\(\pi(\theta|y)\)を正確に復元できる… わけだけど、もちろんそんな密度関数ではうまくいかない(すごく低次元で離散の時ならともかく)。
妥協の方向は2つある。
- 標準的なカーネル密度関数\(K\)を使って$$ \pi_\epsilon (y|x, \theta) = \frac{1}{\epsilon}K\left(\frac{|x-y|}{\epsilon}\right)$$
- 一度低次元な要約統計量\(T(x)\)に落として$$ \pi_\epsilon (y|x, \theta) = \frac{1}{\epsilon}K \left(\frac{|T(x) – T(y)|}{\epsilon}\right)$$ これは\(y\)が不完全でもいけるというメリットがある。
たとえば、カーネルとして一様密度を考え、要約統計量\(T(x)\)が\(T(x)\)を中心とした半径\(\epsilon\)の球に一様にちらばるとすると[←???]、2番目のアイデアはこう書ける。$$ \pi_\epsilon(y|x, \theta) \propto \begin{cases} 1, & \mathrm{if} \ \rho(T(x), T(y)) \leq \epsilon \\ 0, & \mathrm{otherwise} \end{cases} $$ ほかにもいろいろある。
[最後の例はよくわからんかったが、要するに、要約統計量とカーネル密度をいいかんじに決めればよい、という話であろう]
2.3 単純な例
[理解を確かめるために目を通したので、途中で頭が混乱してしまった… 逐語訳する]
例示として、ある単純な例における目標事後密度から尤度フリー近似を導出してみよう。
\(\pi(\theta|y)\)が単変量密度\(N(0,1)\)だという場合について考える。この事後分布を尤度フリーの状況で実現するために、尤度を\(x \sim N(\theta, 1)\)とし、\(\theta\)の十分統計量(標本平均)を\(T(x) = x\)と定義する。観察データは\(y=0\)とする。簡略のため、事前分布は\(\pi(\theta) \propto 1\)とする。
[ちょっと待って！ ここまでの話を私にもわかるように書き換えよう。分析者は\(y=0\)を持っている。分析者は\(y \sim N(\theta, 1)\)だと知っているが、\(\theta\)については全く未知で、事前分布は\(\pi(\theta) \propto 1\)としかいいようがない。事後分布\(\pi(\theta|y\)\)からサンプリングしたい。
そこで尤度フリー近似を行う。まず\(\theta’\)を一様分布からドローする。次に\(x\)を\(N(\theta’, 1)\)からドローする。で、\(x\)をある重み\(\pi_\epsilon(y|x,\theta)\)で受理する。
これからその受理関数というか重みづけ関数について定義します。さあ再開だ]
もし重みづけカーネル\(\pi_\epsilon(y|x,\theta)\)が2.2節の最後の式で表され、\(\rho(T(x), T(y))=|x-y|\)ならば、標準正規累積分布関数を\(\Phi\)として$$ \pi_{LF}(\theta, y) \propto \frac{\Phi(\epsilon-\theta) – \Phi(\epsilon-\theta)}{2\epsilon} $$ となる。
[なにをゆうとるのかというと、\(|x-0| \leq \epsilon\)なら受理します、とゆうとるのである。つまり確率\(Prob(x \leq \epsilon) – Prob(x \leq -\epsilon) = \Phi(\epsilon – \theta) – \Phi (-\epsilon -\theta)\)で受理しますというとるのである。あっれーーー？ 分母の\(2\epsilon\)ってなんなの??? よくわかんなくなってきてしまった…]
いっぽう、もし重みづけカーネル\(\pi_\epsilon(y|x,\theta)\)がガウス密度\(y \sim N(x, \epsilon^2/3\)ならば、$$ \pi_{LF}(\theta, y) = N\left(0, 1+\frac{\epsilon^2}{3}\right)$$ となる。
実際に近似してみると、\(\epsilon\)を小さくすればよい近似になる。2つの方法の間にはあまり差がない。
[後略]
3. 尤度フリーMCMCサンプラー
尤度フリーなMetroplis-Hastingsサンプラーについて。
提案分布を$$ q[(\theta, x), (\theta’, x’)] = q(\theta, \theta’)\pi(x’|\theta’) $$ と分解する。つまり、現在状態を\((\theta, x)\)として、提案分布\(q(\theta, \theta’)\)から\(\theta’\)をドローし、\(\pi(x|\theta’)\)から\(x’\)をドローする。定常分布\(\pi_{LF}(\theta, x|y)\)を持つマルコフ連鎖を作るため、詳細釣り合い条件 $$ \pi_{LF}(\theta, x|y)P[(\theta, x), (\theta’, x’)] = \pi_{LF}(\theta’, x’|y)P[(\theta’, x’), (\theta, x)]$$ を課したい。
これを満たすのがこういうアルゴリズムである。まず\(\theta_0, x_0\)を初期化して、
- 提案分布\(q(\theta_t, \theta)\)から\(\theta’\)をドローする。
- \(\theta’\)のもとでモデル\(\pi(x|\theta’)\)から\(x’\)をドローする。
- 受理するかどうか決める。受理確率は$$ \min\left( \frac{\pi_\epsilon(y|x’, \theta’) \pi(\theta’) q(\theta’, \theta_t)}{\pi_\epsilon(y|x, \theta) \pi(\theta) q(\theta, \theta’)} \right) $$ とする。受理したら\((\theta_{t+1}, x_{t+1}) = (\theta’, x’)\), しなかったら\((\theta_{t+1}, x_{t+1}) = (\theta, x)\)。
- \(t = t+ 1\)として最初に戻る。
このアルゴリズムは詳細釣り合い条件を満たす。なぜなら…[略]
[…中略…]
これが基礎的なLF-MCMCアルゴリズムである。ここからはそのバリエーションを紹介しよう。
3.1 周辺状態サンプラー
最初に述べたように、究極的に関心があるのは、ふつう周辺事後分布$$ \pi_{LF}(\theta|y) \propto \pi(\theta) \int_\mathcal{Y} \pi(y | x, \theta) \pi(x | \theta) dx $$ ですわね。いま\(x\)の独立ドローが\(S\)個あったら、モンテカルロ積分$$ \pi_{LF} \approx \frac{\pi(\theta)}{S} \sum_{s=1}^S \pi_\epsilon(y, x^s, \theta) $$は点ごとの不偏推定になる[\(x^s\)の\(s\)は累乗じゃなくてインデクス]。これを使ってMCMCサンプラーをつくれば、\(\pi_{LF}(\theta|y)\)を直接に目標にできる。[…略…]
[ああなるほど。\(S=1\)なら基礎的なLF-MCMCになるわけだ。この話、まだまだ続くけど、疲れたのでパス]
3.2 誤差分散で補足したサンプラー
\(\epsilon\)は小さいほうがいいんだけれど、そうすると受容率が下がって遅くなる。そこで、\(\epsilon\)をテンパリング・パラメータとみたシミュレーテッド・テンパリングを考える。$$ \pi_{LF}(\theta, x, \epsilon) \propto \pi_\epsilon(y|x, \theta)\pi(x|\theta)\pi(\theta)\pi(\epsilon)$$ と分解することになりますね。
[なるほど。発想は理解した、ということで、パス]
3.3 MCMCサンプラーのさまざまな代替案
ほかにもいろんなアイデアがあるぞ。
- 周辺状態サンプラーのモンテカルロドロー数\(S\)を適応的に変える。
- \(\epsilon\)が異なる複数のマルコフ連鎖をカップリングする。
- \(\epsilon\)じゃなくて、\(T(\cdot)\)を求める際のデータ点の数をテンパリングする。[ふえええ。そんなのありなのか]
4. LF-MCMCのための実践的ガイド
[ここからシミュレーション。とても勉強になりそうな内容ではあるが、当面必要でないのでパス]
5. 考察
本章ではMCMC連鎖に基づく尤度フリーなシミュレーションに焦点をあてたが、\pi_{LF}(\theta|y)\)からサンプリングする方法はほかにもある。\(T(x)\)と\(\theta\)の間の関係を標準的な重回帰で推定する棄却サンプラーとか。尤度フリーな逐次モンテカルロとか、尤度フリーな逐次重点サンプリングとか。
いろんな課題も山積みである。\(T(\cdot)\)の選び方とか、事後分布シミュレーションをもっと効率的にやる方法とか…[略]。
最後に、尤度フリー推論を使ったモデル選択というトレンドもある。云々。
———
話の半分くらいしかわかってないような気がするし、仕事の役には立たなかったけど、尤度フリーという言葉の意味が何となく分かったような気がするから、よしとしよう。
要するに尤度を使わないというだけだから、たとえばM-Hアルゴリズムの受容関数だけ書き換えた尤度フリーMCMC、というようなのがありうるわけだ。きっとSanborn & Griffiths (2007)も尤度フリーMCMCの一種だといえるんでしょうね(Barker受容関数を使っているからではなくて、それを尤度抜きで評価しているから)。…という理解で、あってますでしょうか。