評価できない正規化定数がある場合のベイズ推定について【１】

MHアルゴリズムを適用することができない場合がある

こんにちは，株式会社Nospareの小林です．メトロポリス・ヘイスティングス（MH）アルゴリズムは，未知パラメータの事後分布からサンプリングを行うための標準的なアルゴリズムです．ターゲットの事後分布を$p(\theta|x)=p(x|\theta)p(\theta)/p(x)$と表記することにします．ここで，$p(\theta)$は事前分布，$p(x|\theta)$は尤度関数，$p(x)=\int p(x|\theta)p(\theta)d\theta$を周辺尤度とします．マルコフ連鎖の現在の値を$\theta$とすると，MHアルゴリズムは以下のステップを繰り返します．

1. $\theta'$を提案分布$q(\theta'|\theta)$から生成する．
2. マルコフ連鎖の次の値を，確率

\alpha=\min\left\{1,\frac{p(\theta'|x)q(\theta|\theta')}{p(\theta|x)q(\theta'|\theta)}\right\}. 

で$\theta'$とし，確率$1-\alpha$で$\theta$のままとする．通常多くのモデルでは$p(x)$の評価が困難ですが，現在の設定では$p(x)$が評価できなくても，採択確率は

\alpha=\min\left\{1,\frac{\left[p(x|\theta')p(\theta')/p(x)\right]q(\theta|\theta')}{\left[p(x|\theta)p(\theta)/p(x)\right]q(\theta'|\theta)}\right\}=\min\left\{1,\frac{p(x|\theta')p(\theta')q(\theta|\theta')}{p(x|\theta)p(\theta)q(\theta'|\theta)}\right\}

と書くことができるため，尤度関数と事前密度を評価することができればMHアルゴリズムを適用することが可能です．また尤度関数を$p(x|\theta)=h(x|\theta)/Z$と書くことができるとすると（ここで$Z=\int h(x|\theta)dx$は正規化定数），同様に採択確率は

\alpha=\min\left\{1,\frac{h(x|\theta')p(\theta')q(\theta|\theta')}{h(x|\theta)p(\theta)q(\theta'|\theta)}\right\}

となり，事前密度と尤度関数の正規化されていない部分さえ評価することができればMHアルゴリズムは適用可能です．

ところが一部の統計モデルでは，尤度関数が

p(x|\theta)=\frac{h(x|\theta)}{Z(\theta)},\quad Z(\theta)=\int h(x|\theta)dx

のように，パラメータ$\theta$に依存する正規化定数を含む場合があります．本記事の解説が依拠するPark and Haran (2018)では，そんな統計モデルの例をいくつか挙げており，イジングモデル，指数ランダムグラフモデル，空間上の相互効果のある点過程などが挙げられています．例えばイジングモデルは，$m\times n$の格子状で$x_{ij}$は$-1$か$1$の値を取ります．尤度関数は

p(x|\theta)=\frac{1}{Z(\theta)}\exp\left\{\theta S(x)\right\},\quad S(x)=\sum_{i=1}^m\sum_{j=1}^{n-1}x_{ij}x_{i,j+1}+\sum_{i=1}^{m-1}\sum_{j=1}^nx_{ij}x_{i+1,j}

で与えられます．空間上の依存は$S(x)$によって与えられ，大きな$\theta$の値のときには格子状の相互関係が強くなります．ここで問題になるのが$\theta$に依存する正規化定数$Z(\theta)$で，$x$のすべての値の組み合わせについて$S(x)$を足し合わせたものになります．これは$2^{mn}$通り（！）もあるので，通常$Z(\theta)$の評価は計算時間の問題から難しいでしょう．なぜ$Z(\theta)$の評価ができないことが問題になるかというと，$\theta$の事後分布からMHアルゴリズムでサンプリングする場合には，事後分布は

p(\theta)\propto \frac{h(x|\theta)}{Z(\theta)}p(\theta)

となります．そしてMHアルゴリズムの採択確率は

\alpha=\min\left\{1,\frac{h(x|\theta')p(\theta')Z(\theta)q(\theta|\theta')}{h(x|\theta)p(\theta)Z(\theta')q(\theta'|\theta)}\right\}

となり，尤度の正規化されていない部分$h(x|\theta)$だけでなく，正規化定数$Z(\theta)$の評価も必要になることがわかります．つまり上記のモデルのように$Z(\theta)$が評価できない場合には，標準的なMHアルゴリズムを用いたベイズ推測が行えないということになってしまいます．このような状況をdoubly intractableと呼んだりします．

ではどうするのか

直接評価が困難な$Z(\theta)$を含むモデルに対してなんとかしてMHアルゴリズムを適用するようにするには主に次の2つのアプローチが考えられます．

1. 補助変数をうまく導入することで採択確率内の$Z(\theta)$をキャンセルする
2. $Z(\theta)$を推定してその推定値を使う

本記事では1.の補助変数を導入して$Z(\theta)$をキャンセルすることで，採択確率で$Z(\theta)$の直接の評価を避けることを目指したアルゴリズムを3つ紹介したいと思います．

Møllerのアルゴリズム

Møller et al. (2006)は補助変数を用いた最初のアルゴリズムを提案しました．このアルゴリズムでは，もとのターゲットの分布を補助変数$y$によって拡大した

p(\theta,y|x)\propto p(y|\theta,x)p(\theta)\frac{h(x|\theta)}{Z(\theta)}

をターゲットにします．ここで$y$は$x$と同じ空間$\mathcal{X}$を台としており，この拡大されたターゲット分布の周辺分布はもとのターゲット分布$\int_\mathcal{X} p(\theta,y|x)dy=p(\theta|x)$となるものとします．そして新しいターゲット分布に対する結合提案分布$q(\theta',y'|\theta,y)$は

q(\theta',y'|\theta,y)=q(y'|\theta')q(\theta'|\theta)

のように分解できるものとします．ここで，補助変数$y$に対する提案分布を$q(y'|\theta')=h(y'|\theta')/Z(\theta')$とする（補助変数$y$の提案は$x$と同じ確率モデルから乱数を生成）ことで，採択確率は

\alpha=\min\left\{1,\frac{p(y'|\theta',x)h(x|\theta')p(\theta')Z(\theta)\overbrace{h(y|\theta)/Z(\theta)}^{q(y|\theta)}q(\theta|\theta')}{p(y|\theta,x)h(x|\theta)p(\theta)Z(\theta')\underbrace{h(y'|\theta')/Z(\theta')}_{q(y'|\theta')}q(\theta'|\theta)}\right\}\\
=\min\left\{1,\frac{p(y'|\theta',x)h(x|\theta')p(\theta')h(y|\theta)q(\theta|\theta')}{p(y|\theta,x)h(x|\theta)p(\theta)h(y'|\theta')q(\theta'|\theta)}\right\}

となり$Z(\theta)$が含まれておらず，条件付き密度$p(y|\theta,x)$を評価することができればこの確率を評価することが可能になります．ではこの$p(y|\theta,x)$はどういったものにすればよかというと，単純なチョイスとして$p(y|\theta,x)=h(y|\hat{\theta})/Z(\hat{\theta})$が挙げられます．ここで$\hat{\theta}$は最尤推定値，あるいは最尤推定値の近似値を表します．そうすると採択確率の最初の部分の比は

\frac{p(y'|\theta',x)}{p(y|\theta,x)}=\frac{h(y'|\hat{\theta})/Z(\hat{\theta})}{h(y|\hat{\theta})/Z(\hat{\theta})} = \frac{h(y'|\hat{\theta})}{h(y|\hat{\theta})} 

と評価することができます．

最後に，このアルゴリズムは次のステップを繰り返します：

1. $\theta'$を提案分布$q(\cdot|\theta)$から生成する．
2. $\theta'$が与えられたもとで$y'$をモデルからexactな方法（パーフェクトサンプリング）で生成する．
3. $(\theta',y')$を上記の採択確率に基づいて採択・棄却する．

交換（exchange）アルゴリズム

Møllerのアルゴリズムと同様にMurray et al. (2006)は補助変数を導入して拡大された事後分布をターゲットにしたアルゴリズムを提案しました．ここで補助変数$y$は$h(y|\theta')/Z(\theta')$に従い，$\theta$を所与とした$\theta'$の条件付き分布を$q(\theta'|\theta)$とします（ここで$\theta$はデータ$x$を生成した$\theta$）．このとき結合分布は次のように表されます：

p(\theta,\theta',y|x)\propto p(\theta)\frac{h(x|\theta)}{Z(\theta)}q(\theta'|\theta)\frac{h(y|\theta')}{Z(\theta')}. 

次にアルゴリズムの名前にもなっているように，$\theta$と$\theta'$の位置を交換することを考えます．最初は$x$と$\theta$，$y$と$\theta'$の組み合わせになっていましたが，交換したあとは$x$と$\theta'$，$y$と$\theta$の組み合わせになります．この交換を提案する提案分布を

s(\theta^*,\theta'^{*}|\theta,\theta')=\delta(\theta^*-\theta')\delta(\theta'^{*}-\theta)

と書きます．すると採択確率は

\alpha=\min\left\{1,\frac{s(\theta,\theta'|\theta^*,\theta'^{*}) p(\theta',\theta,y|x)}{s(\theta^*,\theta'^*|\theta,\theta') p(\theta,\theta',y|x)}\right\}=\min\left\{1,\frac{p(\theta',\theta,y|x)}{p(\theta,\theta',y|x)}\right\}\\
=\min\left\{1,\frac{p(\theta')\frac{h(x|\theta')}{Z(\theta')}q(\theta|\theta')\frac{h(y|\theta)}{Z(\theta)}}{p(\theta)\frac{h(x|\theta)}{Z(\theta)}q(\theta'|\theta)\frac{h(y|\theta')}{Z(\theta')}}\right\}=\min\left\{1,\frac{p(\theta')h(x|\theta')q(\theta|\theta')h(y|\theta)}{p(\theta)h(x|\theta)q(\theta'|\theta)h(y|\theta')}\right\}

となり，$Z(\theta)$がキャンセルされます．交換アルゴリズムは次のステップを繰り返します．

1. $\theta'$を提案分布$q(\cdot|\theta)$から生成する．
2. $\theta'$の値のもとで補助変数$y$をモデル$h(\cdot|\theta')/Z(\theta')$からパーフェクトサンプリングで生成する．
3. $\theta'$を上記の採択確率に基づいて採択・棄却する．

ダブルMH（DMH）アルゴリズム

ここまで強調してきませんでしたが，上述の2つのアルゴリズムではパーフェクトサンプリングでexactに補助変数$y$を生成する必要があります．ところがパーフェクトサンプリングは複雑なモデルに対しては適用が困難である場合が多くあります．Liang (2010)が提案したダブルHMアルゴリズムはパーフェクトサンプリングの代わりにMHアルゴリズムを使って$y$を生成します．つまり$\theta$をサンプリングするMHアルゴリズム（外側のサンプラー）の中で$y$をサンプリングするMHアルゴリズム（内側のサンプラー）を実行するので，ダブルMHアルゴリズムと呼ばれます．

まず$T^m_{\theta'}(y|x)$を$\theta'$のもとで$x$から$y$まで$m$ステップで動くMHカーネルとします（不変分布は$h(y|\theta')/Z(\theta')$）．これは詳細釣り合い条件を満たします：

\frac{h(x|\theta')}{Z(\theta')}T_{\theta'}^m(y|x)=\frac{h(y|\theta')}{Z(\theta')}T_{\theta'}^m(x|y) \Longleftrightarrow
h(x|\theta')T_{\theta'}^m(y|x)=h(y|\theta')T_{\theta'}^m(x|y)\\
 \Longleftrightarrow \frac{h(x|\theta')}{h(y|\theta')}=\frac{T_{\theta'}^m(x|y)}{T_{\theta'}^m(y|x)}

これを交換アルゴリズムの採択確率に代入すると，DMHアルゴリズムの採択確率は

\min\left\{1,\frac{p(\theta')T_{\theta'}^m(x|y)h(y|\theta)q(\theta|\theta')}{p(\theta)T_{\theta'}^m(y|x)h(x|\theta)q(\theta'|\theta)}\right\}

となります．DMHアルゴリズムでは，パーフェクトサンプリングを使って$y$をexactに生成する代わりに，$m$ステップの内側のMHアルゴリズムを使います．具体的にはDMHは次のステップを繰り返すことになります：

1. $\theta'$を提案分布$q(\cdot|\theta)$から生成する．
2. $\theta'$のもとで，モデル$h(\cdot|\theta')/Z(\theta')$からMHアルゴリズムの更新を$m$回行って$y$を生成する（$y\sim T_{\theta'}^m(\cdot|x)$）．
3. $\theta'$を上記の採択確率に基づいて採択・棄却する．

3つのアルゴリズムの共通点・相違点・使いやすさ

• 上述のアルゴリズム3つすべての共通点は，補助変数$y$を導入して$Z(\theta)$を採択確率の中でキャンセルすることを目指しています．そして$y$は（exactであれ，inexactであれ）$x$と同じモデルから生成されます．
• 3つのうち，Møllerのアルゴリズムと交換アルゴリズムは$y$の生成にパーフェクトサンプリングが必要になります．一方で，DMHは$y$の生成にMHアルゴリズム（exactではない）を使用します．
• Møllerのアルゴリズムと交換アルゴリズムは詳細釣り合い条件を満たすので漸近的にexactな手法になりますが，DMHは内側のサンプラーのステップ数$m$が無限に行かない（実際の適用では有限）と詳細釣り合い条件が満たされないので漸近的にexactな手法ではありません．
• 3つの中で一番適用が簡単なのは，MHアルゴリズムだけで構成されるDMHです．その次がパーフェクトサンプリングが必要となる交換アルゴリズムです．そしてMøllerのアルゴリズムで，これは最尤推定値あるいはそれの近似値が必要になるからです．

おわりに

Doubly intractableな状況に直面するのはあまり頻繁ではないかもしれませんが，空間モデルやネットワークモデルなど特定のモデルに対するベイズ推定を行う際には本記事で紹介したアルゴリズムは有用となるでしょう．各アルゴリズムに対する細かいチューニングの方法は本記事のスコープ外となりますので，Park and Haran (2018)や各アルゴリズムの元論文を参照してください．今後アルゴリズムの実装や$Z(\theta)$を推定するアプローチについての解説なども行っていければと思います．　

参考文献

1. Liang, F. (2010), A Double Metropolis–Hastings Sampler for Spatial Mod- els With Intractable Normalizing Constants, Journal of Statistical Computation and Simulation, 80, 1007–1022.
2. Møller, J., Pettitt, A. N., Reeves, R., and Berthelsen, K. K. (2006), An Efficient Markov Chain Monte Carlo Method for Distributions With Intractable Normalising Constants, Biometrika, 93, 451-458.
3. Murray, I., Ghahramani, Z., and MacKay, D. J. C. (2006), MCMC for Doubly-Intractable Distributions, in Proceedings of the 22nd Annual Conference on Uncertainty in Artificial Intelligence (UAI-06), AUAI Press, 359–366.
4. Park, J. and Haran, M. (2018), Bayesian Inference in the Presence of Intractable Normalizing Functions, Journal of the American Statistical Association, 113, 1372-1390.

---

一緒にお仕事しませんか！
株式会社Nospareでは統計学の様々な分野を専門とする研究者が所属しており，新たな知見を日々追求しています．統計アドバイザリーやビジネスデータの分析につきましては弊社までお問い合わせください．