MHアルゴリズムを適用することができない場合のベイズ推定
こんにちは，株式会社Nospareの小林です．前回の記事では，正規化定数が未知パラメータに依存しているため，通常のメトロポリス・ヘイスティングス（MH）アルゴリズムの適用が不可能な場合に利用することができるMCMCアルゴリズムを紹介しました．ここでの設定は，
尤度関数が
p(x|\theta)=\frac{h(x|\theta)}{Z(\theta)},\quad Z(\theta)=\int h(x|\theta)dx
のように書け，正規化定数$Z(\theta)$がパラメータ$\theta$に明示的に依存しており，これが解析的に得られない，もしくは計算的に評価することが困難であるものとします．本記事では，MHアルゴリズムの採択確率の中で扱いの困難である$Z(\theta)$をうまくキャンセルさせる方法の代わりに，これを近似して使う方法について紹介します．
Atchade, Larillot and Robert（ALR）の方法
Atchade, Larillot and Robert (2008)の方法では，$Z(\theta)$を重点サンプリングとstochastic approximationで近似してMHアルゴリズムの採択確率に使用します．
まず，
\begin{split}
Z(\theta)&=\int h(x|\theta)dx=\int h(x|\theta)\frac{Z(\theta^{(0)})}{h(x|\theta^{(0)})}\frac{h(x|\theta^{(0)})}{Z(\theta^{(0)})}dx\\
&=Z(\theta^{(0)})\int\frac{ h(x|\theta)}{h(x|\theta^{(0)})}\frac{h(x|\theta^{(0)})}{Z(\theta^{(0)})}dx
\end{split}
と書けるため，ある$\theta^{(0)}$に固定した確率モデルを提案分布とした重点サンプリングの形になっています．ここでもし$\theta$と$\theta^{(0)}$が遠いと，この近似のパフォーマンスは悪いものとなります．そこでALRでは，複数の値$\theta^{(1)},\dots,\theta^{(d)}$を用意し，重点サンプリングの推定値に$\theta$と$\theta^{(i)}$の距離に応じたウェイト$k(\theta,\theta^{(i)})$を導入します：
Z(\theta)=\sum_{i=1}^d k(\theta,\theta^{(i)})Z(\theta^{(i)})\int\frac{ h(x|\theta)}{h(x|\theta^{(i)})}\frac{h(x|\theta^{(i)})}{Z(\theta^{(i)})}dx
ALRでは重点サンプリングとstochastic approximationを繰り返し行っていき，$n+1$回目の反復にける推定値は
\hat{Z}_{n+1}(\theta)=\sum_{i=1}^dk(\theta,\theta^{(i)})\exp(c_{n+1}^{(i)})\frac{\sum_{k=1}^{n+1}\frac{h(x_k|\theta)}{h(x_k|\theta^{(i)})}I(I_k=i)}{\sum_{k=1}^{n+1}I(I_k=i)}
で与えられます．ここで，$I_n\in\lbrace1,\dots,d\rbrace$は重点サンプリングにどの$\theta^{(i)}$の値を使うかを示す変数，$c_n^{(i)}$は（stocastic approximationによって）$\log Z(\theta^{(i)})$に収束する変数です．
Pseudo-marginal（PM）アプローチ
以前この記事で紹介したPseudo-marginalアプローチ（Andrieu and Roberts, 2009）では，MHアルゴリズムの採択確率の中にある尤度関数$p(x|\theta)$を不偏推定量$\hat{p}(x|\theta)$で置き換えて$\theta$のサンプリングを行います．パラメータに依存する$Z(\theta)$があるここの文脈においては，$1/Z(\theta)$の不偏推定値が必要になります．$Z(\theta)$の不偏推定値$\widehat{Z(\theta)}$は確率モデルからのMCMCサンプルによって得ることはできますが，$1/\widehat{Z(\theta)}$は$1/Z(\theta)$に対して一致性をもちますがバイアスがあるので補正を行います（Lyne etal. 2015）．このアプローチでは$\theta$のMCMCサンプリングの各反復において確率モデルからのMCMCサンプリングとバイアス補正が必要になるので，結果として計算負荷は高くなってしまいます．
Noisy MCMC
マルコフ連鎖の推移核$P$がターゲット分布$p$に対して詳細釣り合い条件を満たすときにはasymptotically exactといいます．一方で，$P$が$\hat{P}$によって近似し，$\hat{P}$を使って生成したサンプルは$\hat{p}$に従うことになります．このようなアルゴリズムのことをnoisy MCMCと呼び，noisy MCMCは上で紹介したALRやPMアプローチなどを含む幅広いクラスのアルゴリズムです．
前回の記事で紹介した交換アルゴリズムでは確率モデルから補助変数$y$を1つだけ生成しましたが，これは$Z(\theta)/Z(\theta')$の1つだけサンプルを使った重点サンプリングによる不偏推定値として見ることができます．ここでもし，1つだけ$y$を生成する代わりに複数の$y_1,\dots,y_N$を確率モデル$h(y|\theta')/Z(\theta')$から生成すると，重点サンプリング推定値の分散は小さくなるでしょう．このnoisy交換アルゴリズムの採択確率は
\min\left\{1,\frac{p(\theta')h(x|\theta')q(\theta|\theta')}{p(\theta)h(x|\theta)q(\theta'|\theta)}\frac{1}{N}\sum_{i=1}^N\frac{h(y_i|\theta)}{h(y_i|\theta')}\right\}
で与えられます．このアルゴリズムは詳細釣り合い条件を満たさないのでaymptotically inexactですが，通常の交換アルゴリズムよりも効率的であることが知られています．
関数エミュレータ
ALRと同様に複数の$\theta^{(i)}$のもとで以下の重点サンプリング推定値を考えます：
\log \hat{Z}(\theta^{(i)})=\log\left(\frac{1}{N}\sum_{l=1}^N\frac{h(x_l|\theta^{(i)})}{h(x_l|\tilde{\theta})}\right)
ここで$\tilde{\theta}$は最尤推定値の近似値であるとします．Park and Haran (2020)は，これら$d$個の推定値をデータとしてガウス過程回帰による関数エミュレーションを提案しました：
\log \hat{Z} = \psi\beta + u
ここで$\psi=(\theta^{(1)},\dots,\theta^{(d)})$, $\beta$は線形回帰係数，$u$はガウス過程に従います．そうすると，新しい$\theta^*$の値のもとでの$\log \hat{Z}(\theta^\star)$の値は，$\log \hat{Z}(\theta^\star)|\log \hat{Z}$の条件付分布（正規分布です）から求めることができます．同じ方法を使って，$\log Z(\theta)$に対するエミュレータの代わりに，尤度関数$p(x|\theta)$全体に対するエミュレータの構築も行うことができます．
おわりに
本記事ではMHアルゴリズムの中に出てくる$Z(\theta)$に対して何らかの近似値を使う方法を紹介しました．$Z(\theta)$をうまくキャンセルさせようとする前回の記事との違いは，（本記事では触れませんでしたが）やはり近似を用いた場合の理論的な正当化が必要になるということです．近似値を用いたほうが計算負荷が軽くなるのかというとそうとは限らず，例えば本記事の解説が依拠していたPark and Haran (2018)の数値実験ではdouble MHが一番高速で，PMアプローチが一番遅い結果となっていました．Park and Haran (2020)では関数エミュレータとDMHの比較がなされており，エミュレータは同等かそれ以上のパフォーマンスを持っていました．エミュレータの場合は，一旦ガウス過程のエミュレータの構築を行ってしまえば，新しい$\theta$の値での$\log Z(\theta)$の値が高速に計算できるという利点があります．
参考文献
- Andrieu, C., and Roberts, G. O. (2009), The Pseudo-Marginal Approach for Efficient Monte Carlo Computations, The Annals of Statistics, 37, 697–725.
- Atchade, Y., Lartillot, N., and Robert, C. P. (2008), Bayesian Computation for Statistical Models With Intractable Normalizing Constants, Brazilian Journal of Probability and Statistics, 27, 416–436.
- Park, J. and Haran, M. (2018), Bayesian Inference in the Presence of Intractable Normalizing Functions, Journal of the American Statistical Association, 113, 1372-1390.
- Park, J. and Haran, M. (2020), A Function Emulation Approach for Doubly Intractable Distributions, Journal of Computational and Graphical Statistics, 29, 66-77.
一緒にお仕事しませんか！
株式会社Nospareでは統計学の様々な分野を専門とする研究者が所属しており，新たな知見を日々追求しています．統計アドバイザリーやビジネスデータの分析につきましては弊社までお問い合わせください．
