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

モデル不確実性と周辺尤度の計算

まず２つのモデル空間が$\mathcal{M}_1$と$\mathcal{M}_2$からなる単純な場合を考えます．ここでそれぞれのモデルの事前確率を$p(\mathcal{M}_1)$，$p(\mathcal{M}_2)$とします．ベイズの定理によって事後オッズ比（posteriro odds ratio）は

\frac{p(\mathcal{M}_1|y)}{p(\mathcal{M}_2|y)} = \frac{p(\mathcal{M}_1)}{p(\mathcal{M}_2)}\times\frac{p(y|\mathcal{M}_1)}{p(y|\mathcal{M}_2)}

と書けます．ここで

p(y|\mathcal{M}_k) = \int_{\Theta_k} p(y|\theta_k,\mathcal{M}_k)p(\theta_k|\mathcal{M}_k)d\theta_k,\quad k=1,2

で，$\theta_k$はモデル$\mathcal{M}_k$の未知メータ，$p(y|\theta_k,\mathcal{M}_k)$と$p(\theta_k|\mathcal{M}_k)$はそれぞれ尤度関数とモデル$M_k$が所与としたときの事前分布を表します．モデル$\mathcal{M}_k$の周辺尤度，データの周辺分布あるいはエビデンスと呼ばれる$p(y|M_k)$の比をベイズファクターと呼びます．また$p(\mathcal{M}_1|y)+p(\mathcal{M}_2|y)=1$です．

モデル選択は決定論の考えに基づいて事後期待損失を最小化するように行われます．具体的には，$p(\mathcal{M}_1|y)/p(\mathcal{M}_2|y)$が第1と第2の過ちによる損失の比を超える場合にモデル$\mathcal{M}_1$が選ばれるようになります．モデル平均は関心のある数量の事後期待値を両方のモデルについて求める際に$p(y|\mathcal{M}_1)$と$p(y|\mathcal{M}_2)$をウェイトとして用いることになります．

いずれにせよ，重要となるのは周辺尤度$p(y|\mathcal{M}_k)$の評価が必要になることです．通常これが解析的に得られる場合というのは非常に限られています．変分ベイズ（VB）ではエビデンスの下限が副産物として得られますが，ここでは周辺尤度を直接評価する方法を取り扱います．

周辺尤度における積分は，$\mathcal{M}_k$について定義されたの事前期待値として見ることができます．事後期待値ではなく事前期待値であることは，

$p(\theta_k|\mathcal{M}_k)$が適切な密度関数である場合にのみ定義されるので，ベイズファクターと事後オッズ比は，無情報（典型的には非正則）な事前分布または客観的な事前分布の元では計算することができない
シミュレーションによって周辺尤度を計算する最も直接的な方法，つまり$p(\theta_k)$から$\theta_k$を生成し，それについて$p(y|\theta_k)$の平均を取る方法は，$p(\theta_k)$が必ずしも尤度関数の値が大きい領域でで高い密度を持つとは限らないため，一般的に不正確なものとなる

という問題があり，特にシミュレーションベースの推定については，精度を向上させるために何らかの方法でデータから情報を得た$\theta_k$の生成を行うことが必要になります．

重点サンプリング（Importance sampling, IS）

周辺尤度を計算するのにISを利用することができます．

w^{(i)}=\frac{p(y|\theta_k^{(i)})p(\theta_k^{(i)})}{q(\theta_k^{(i)}|y)},\quad i=1,\dots, M

とし，提案分布$q(\cdot|y)$から$M$個のシミュレーションを生成するとIS推定量

\hat{p}(y|\mathcal{M}_k)=\frac{1}{M}\sum_{i=1}^M w^{(i)}

が得られます．ここでは$p(y|\theta_k^{(i)},\mathcal{M}_k)$，$p(\theta_k^{(i)})$, $q(\cdot|y)$が既知であり，追加的な正規化が必要ないものとします．

ISの考え方は

提案分布からのシミュレーションを使って加重平均を計算する
ウェイトにある提案分布の密度関数を評価する

という２つの側面があります．Reciprocal IS（RIS）という方法では２つ目の側面を利用する一方でシミュレーション自体は事後分布自身から行います．$p(\theta_k|y)$のサポートに含まれる$q(\cdot|y)$について，$g(\theta_k)=q(\theta_k|y)/[p(y|\theta_k)p(\theta_k)]$と定義すると，

E[g(\theta_k)|y]=[p(y|\mathcal{M}_k)]^{-1}

という関係を得ることができます．これから，$p(\theta_k|y)$からのシミュレーションを元に上の式の逆数を周辺尤度の推定値として用いることができます．調和平均推定量と呼ばれるものは，RISの$q(\theta_k|y)=p(\theta_k)$という特別な場合に相当します．ただし事前分布の裾が事後分布よりも厚い場合には，調和平均はまったく使えないものであるということも知られています．RISの考え方に基づいてブリッジサンプラー（事後分布と提案分布の両方からシミュレーションを行う）や事後分布からのサンプリングを行う際に事後分布を変分近似したものを提案分布として用いることなども提案されています．

MCMC

事後分布$p(\theta_k|y)$の定義から，周辺尤度は

p(y|\mathcal{M}_k)=\frac{p(y|\theta_k)p(\theta_k)}{p(\theta_k|y)}

というように表すこともできます．Chib (1995)の手法は，上の等式を事後分布の密度が高い値$\theta^\ast_k$で上の式を評価します．分子の尤度関数と事前密度は通常簡単に評価することができるのですが，事後密度である分母の評価は直接の評価は通常不可能です．ただし，$\theta_k^\ast=(\theta^\ast_{k,1},\dots,\theta^\ast_{k,p_k})$と書くと，$\theta_k\ast$で評価された同時事後分布は

p(\theta_k^\ast|y)=p(\theta_{k,1}^\ast|\theta_{k,2}^\ast,\dots,\theta_{k,p_k}^\ast,y)p(\theta_{k,2}^\ast|\theta_{k,3}^\ast,\dots,\theta_{k,p_k}^\ast,y)\cdots p(\theta_{k,p_k}^\ast|y)

と分解できます．右辺の最後の項$p(\theta_{k,p_k}^\ast|y)$は通常のMCMCのアウトプットから評価することができます．次に$p(\theta_{k,p_{k-1}}^\ast|\theta_{k,p_{k}}^\ast,y)$は$\theta_{k,p_{k}}^\ast$を固定したMCMCを適用する（$\theta_{k,p_{k}}$をサンプリングせずに固定する）と評価することができます．これを繰り返して一つずつパラメータを固定したMCMCを適用することで周辺尤度を評価することができます．サンプリングにおいてはギブスだけでなくMHアルゴリズムなども利用できます．

リバーシブルジャンプ（reversible jump）MCMC

上述の手法では周辺尤度を各モデルにおいて個別に計算するということが行われていました．一度各周辺尤度が得られると，モデルの事後確率の計算，モデル選択，モデル平均などができる，という感じです．

一方で，モデルの不確実性を利用して未知パラメータを拡張し，この拡張された空間をターゲットとしてサンプラーを設計することもできます．これは，回帰モデルの変数選択や有限混合モデルのコンポネント数など，様々な論文で採用されている原理でもあります．ここではGreen (1995)のアプローチにのみ焦点を当てます．

Green (1995) は，未知のモデルの問題をどのモデルが使われているかによって（モデル固有の）未知の次元が変化する問題と特徴づけ，異なる次元のパラメータ空間の間を移動できるリバーシブルジャンプMCMC（Reversible jump, RJMCMC）アルゴリズムを提案しました．RJMCMCは基本的にはデータ拡大法（data augmentation）の1つであるとみることができます．

いま候補モデルの集合が$\mathcal{M}=\{ \mathcal{M}_1,\mathcal{M}_2,\dots\}$で与えられ，それらのインデックスを$k=1,2,\dots$で表すとします．各モデルには$p_k$次元の未知パラメータ$\theta_k$があるとします．RJMCMCアルゴリズムは

p(k,\theta_k|y)=\frac{p(y|k,\theta_k)p(\theta_k|k)p(k)}{p(y)}

をターゲットとし，一時的な補助変数を導入して拡張された空間の次元を同じにすることで，任意の2つのパラメータ空間の間を移動することができるようになります．$\{k,\theta_k\}$の拡張された空間からのシミュレーションを用いて任意の事後期待値$E[g(\theta_k)|y]$を計算することができます．また$k$番目のモデルの事後確率$p(k|y)$を含む期待値$E[g(k)|y]$も計算することができるため，RJMCMCの副産物としてモデルの事後確率を計算することが可能になります．

予測

$p(y|\theta)$が真のデータ生成過程と一致しているという仮定のもと，ベイズ統計における予測は

p(y_{n+1}^\ast|y)=\int_\Theta p(y_{n+1}^\ast|\theta,y)p(\theta|y)d\theta

という事後期待値に基づいて行います．分布$p(y_{n+1}^*|y)$は，仮定されたモデル（条件付き予測$p(y_{n+1}^\ast|\theta, y)$）と事後分布の構造に基づいている）と$p(\theta|y)$に情報を与える事前情報の両方を条件として，$y_{n+1}$に関するすべての不確実性を要約しています．$y_{n+1}$の点予測および区間予測，そして実際に他の任意の分布の概要を上の定義から抽出することができます．モデル自体が不確かであり，モデル空間に有限個のモデル$\mathcal{M}_1, \mathcal{M}_2,...,\mathcal{M}_K$が存在すると仮定した場合，モデル平均の原理に基づいて予測値を以下のように得ることができます：

p_{MA}(y_{n+1}^\ast|y)=\sum_{k=1}^Kp(y_{n+1}^\ast|y,\mathcal{M}_k)p(\mathcal{M}_k|y)

未来

いままで約250年に渡ってベイズ統計における数値計算の手法はとても大きく進歩してきましたが，今後もまだまだたくさんやるべきことがあります．国際会議Bayes Comp 2020では「スケーラビリティ」，「大きなデータセット」，「高次元問題」などをキーワードにする報告が多くあり，確率的トピックモデル，医学，脳機能イメージング，遺伝学，生物学，疫学，生態学，心理学，計量経済学など様々な応用分析が取り上げられていました．またモデル誤特定のもとでの数値計算というトピックも挙げられており，

パラメトリックモデルが誤特定されていた場合の示唆は？
従来の尤度ベースのアプローチが使われない場合における数値計算上の示唆は？

といった点についても近年は研究がなされています．

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