慶應義塾大学・株式会社Nospareの菅澤です．
今回はベイズ的に複数モデルを組み合わせる方法として，Bayesian model averaging (BMA) とスタッキング (Stacking) について解説します．
はじめに
データ分析の際には複数のモデル候補が存在するケースが多々あります．例えば，回帰モデルにおいて複数の変数の組み合わせでモデルを立てることを考えると自然とそのような状況になります．
その際に，大きく2つのアプローチが考えられます．
- 何らかのモデル評価規準 (AICやBICなど) を用いてベストなモデルを1つ選ぶ
- 候補モデルを組み合わせて1つのモデルを作る
後者のアプローチは一般にモデル平均と呼ばれており，前者のようなモデル選択に基づく方法よりも予測精度が高くなる傾向が知られています．
モデル平均自体にも様々な方法が提案されておりますが，今回はベイズ的なアプローチに焦点を当て，代表的な方法であるBMAとスタッキングについて解説していきます．
Bayesian model averaging (BMA)
ベイズ的にモデルを平均化する最も代表的な方法として，ベイズモデル平均化 (Bayesian model averaging, BMA) があります．この方法は「モデルが選ばれる確率」を考えることで，モデルが選ばれる事後確率に基づく重み付き平均を考える方法です．
観測データ$y=(y_1,\ldots,y_n)$に対して$K$個のモデル$M_1,\ldots,M_K$を想定しているとします．さらに，モデル$M_k$に対して以下を定義しておきます．
- $\theta_k$: モデル$M_k$のパラメータ (モデル$M_k$ごとに次元が異なる可能性がある)
- $p(\theta_k|M_k)$: パラメータ$\theta_k$の事前分布
- $f(y|\theta_k, M_k)$: モデル$M_k$の尤度関数
$M_k$が所与の (モデルを1つ決めた) もとで，$\theta_k$の事後分布は
\pi(\theta_k|y, M_k)=\frac{\pi(\theta_k|M_k)f(y|\theta_k, M_k)}{\int \pi(\theta_k|M_k)f(y|\theta_k, M_k)d\theta_k}, \ \ \ \ \ k=1,\ldots,K
となります．
次に，$M_k$の事後分布（各モデルが選ばれる事後確率）を考えます．$p(M_k)$をモデル$M_k$が選ばれる事前確率（分布）とすると，事後確率は以下のように表現できます．
p(M_k|y)=\frac{p(y|M_k)p(M_k)}{\sum_{k'=1}^K p(y|M_{k'})p(M_{k'})}.
ここで，$p(y|M_k)$はモデル$M_k$の周辺尤度 (marginal likelihood)
p(y|M_k)=\int p(y|\theta_k, M_k)p(\theta_k|M_k)d\theta_k
です．ベイズ統計学において周辺尤度はエビデンスとも呼ばれており，モデル評価をする際に重要な量です．例えば，2つのモデルの周辺尤度の比 $p(y|M_k)/p(y|M_{k'})$はベイズファクターと呼ばれます．モデルの事前確率として一様分布$p(M_k)\propto 1$を用いた場合，事後確率は相対的な周辺尤度の大きさで決まることがわかりますので，事後確率はベイズファクターと似たような量を考えていることがわかります．
この事後確率を使って，モデル$M_1,\ldots,M_K$を統合して予測することを考えます．ます，モデル$M_k$のもとで，未観測データ$\tilde{y}$に対する事後予測分布は
p(\tilde{y}|y,M_k)= \int p(\tilde{y}|\theta_k,M_k)\pi(\theta_k|y, M_k)d\theta_k
となります．これは，各モデルごとに（MCMCなどを使って）事後分布$\pi(\theta_k|y, M_k)$からのサンプルを生成することができれば，モデル$p(\tilde{y}|\theta_k,M_k)$に従って$\tilde{y}$を生成することで簡単に事後予測分布からのサンプルを生成することができます．
BMAは$\tilde{y}$に対する予測分布を
p(\tilde{y}|y)=\sum_{k=1}^K p(\tilde{y}|y,M_k)p(M_k|y)
の形で与えます．この形は，複数ある予測分布をそれぞれの事後確率をウエイトにした重み付き平均を考えていると解釈できます．
これにより，モデルを1つ選ぶのではなく，確率的な枠組みのもとでモデル選択の不確実性を考慮してモデルを平均化することが可能になります．
BMAの性質として，候補モデル$M_1,\ldots,M_K$のどれかが真のデータ生成過程に一致している (どれかはわからない) 状況では漸近的に妥当な方法として知られています．一方で，そうでない状況（いわゆるM-openと呼ばれる状況）では漸近的にKullback-Leiblerダイバージェンスが最も小さくなるモデルに対するウエイトが1になる (それ以外のモデルは0になってしまう) ことが知られています．
これは応用上とても困る問題です．複数のモデルをバランスよく組み合わせて予測を安定化させるのがモデル平均の基本的なモチベーションですが，BMAでは極端なウエイトが返ってきてしまう可能性がある（本質的にモデル選択をしているのと変わらない結果になる）ことを意味しています．
また，実用上の問題として周辺尤度の計算が簡単ではないという問題もありますので，ベイズ的なモデル平均の概念としては自然なBMAですが，計算コストや実際の数値的なパフォーマンスとしては一般的にあまり使いやすい手法ではありません．（もちろん周辺尤度が簡単に計算できる線形回帰モデルなどではまた話が変わってきます．）
スタッキング (Stacking)
スタッキング (Stacking) は本来複数モデルからの推定値を統合するアルゴリズムですが，Yao et al. (2018) ではスタッキングを用いた予測分布の統合手法が提案されています．
BMAと同様に，各モデルごとの事後予測分布$p(|y, M_k)$を凸結合した予測分布$\sum_{k=1}^K w_kp(\tilde{y}|y, M_k)$を考えます．(BMAでは$w_k$として事後モデル確率を用いていました．)
スタッキングでは，$w_k$を交差検証の考え方を用いて求めます．$y_{-i}=(y_1,\ldots,y_{i-1},y_{i+1},\ldots,y_n)$をデータ$y$から$y_i$のみを除外したデータセットとし，Leave-one-out (LOO) 予測分布を
p(y_i|y_{-i}, M_k)= \int p(y_i|\theta_k, M_k) p(\theta_k|y_{-i}, M_k)d\theta_k
と定めます．これは，$y_{-i}$のみから得られる$\theta_k$の事後分布から$y_i$の事後予測分布を構成したものです．
このとき，ウエイト$w_1,\ldots,w_K$を用いて構成した統合予測分布のLOO予測がうまくいっているかを測る尺度として
\sum_{i=1}^n \log \bigg\{\sum_{k=1}^K w_kp(y_i|y_{-i}, M_k)\bigg\} \tag{1}
を考えます．これは予測がうまくいっているほど$p(y_i|y_{-i}, M_k)$（テストデータ$y_i$の尤度）の値が高くなるため，大きい値ほど予測がうまくいっていることを表します．
そこで，最適なウエイト$w_1,\ldots,w_K$は上記の関数を制約条件$\sum_{k=1}^Kw_k=1$のもとで最大化することで得られます．
この方法を単純に実行する場合は計算上の大きな問題点が生じます．それは，(1)式を計算するのにLOO予測分布$p(y_i|y_{-i}, M_k)$を$n$回計算する必要がある点です．これは，各$y_{-i}$ごとに$\theta_k$の事後分布を（マルコフ連鎖モンテカルロ法などで）計算する必要があることを意味します．
Yao et al. (2018)では，この問題に対して重点サンプリングによる効率的な近似計算手法を与えています．
全データによる事後分布$p(\theta_k|y, M_k)$を用いて，LOO予測分布は
p(y_i|y_{-i}, M_k)= \int p(y_i|\theta_k, M_k) \frac{p(\theta_k|y_{-i}, M_k)}{p(\theta_k|y, M_k)}p(\theta_k|y, M_k)d\theta_k
と表すことができるため
\frac{p(\theta_k|y_{-i}, M_k)}{p(\theta_k|y, M_k)}\propto \frac{1}{p(y_i|\theta_k, M_k)}
を重点ウエイトとした重点サンプリングによって計算することができます．
スタッキングの理論的な性質としては漸近的なウエイトの最適性が知られています．具体的には，(1)式を最適にするウエイトは，漸近的に重み付き密度と真の分布のKullback-Leiblerダイバージェンスが最小になるようなウエイトに一致します．
BMAが単に事後確率に関する重み付けだったのに対して，スタッキングは重みの選択に関する最適性を保持するため，BMAが得意でないケース (M-openのケース) などで効果的であると考えられます．
おわりに
株式会社Nospareでは統計学の様々な分野を専門とする研究者が所属しております．統計アドバイザリーやビジネスデータの分析につきましては株式会社Nospareまでお問い合わせください．
