東京大学・株式会社Nospareの菅澤です．今回から前編・後編に分けて，標準的なベイズ推測の枠組みを拡張した 一般化ベイズ(general Bayes)法 について紹介します．

標準的なベイズ推測

$y_1,\ldots,y_n$（以下ではまとめて$Y_n$と表記）をパラメータ$\theta$の確率分布$f(x|\theta)$からのランダムサンプルとします．$\theta$の事前分布を$\pi(\theta)$とすると，$\theta$の事後分布は

\pi(\theta \mid Y_n) \propto \underbrace{\pi(\theta)}_{事前分布} \times \underbrace{\prod_{i=1}^n f(y_i \mid \theta)}_{尤度}

で与えられます．この手順は，事前分布（事前信念）にデータの情報を尤度という形で与え，分布（信念）を更新していると捉えることができます．

このような標準的な事後分布の応用上の制約としては 「$\theta$に対する推測をするのにデータの確率分布を特定化する必要がある」 という点です．問題によってはこの特定化が困難であったり，局外パラメータ(興味のないパラメータ)として省略したい場合があります．(具体例は後述します．)

一方で，頻度論的な方法論では，尤度を用いた推定方法だけでなく，M推定量のように損失関数を自由にデザインすることで様々な推定方法が扱うことが可能です．ただし， 損失関数の形によっては推定量を計算するための最適化が困難であったり，漸近理論による推定結果の不確実性評価が容易でない問題 が頻繁に発生することにも注意が必要です．

ベイズ法は自然な不確実性評価の枠組みを事後分布という形で与えることができ，事後分布からのサンプリング（最適化に依らない方法）によって推測を行うことができます．したがって，尤度に拘らない一般的な損失関数を使ってベイズ推測を実行する枠組みがあれば 「M推定量のような汎用性や利便性を享受しつつ，不確実性評価まで自然な枠組みで実行できるベイズ推測の枠組み」 ができることになります．

これがまさに一般化ベイズ法です．

一般化事後分布

$L(Y_n,\theta)$をパラメータ$\theta$を推定するための一般の損失関数とします．すなわち，$\theta$のM推定量が$\text{argmin}_\theta L(Y_n,\theta)$で与えられるとします．

このとき一般化事後分布は以下のように定義されます．

\pi(\theta\mid Y_n) \propto \underbrace{\pi(\theta)}_{事前分布} \times \underbrace{\exp(-\omega L(Y_n,\theta))}_{\text{損失関数によるデータの情報}}

この形の事後分布はGibbs posteriorやquasi-posteriorと呼ばれることもあります．

ここで，$\omega$はlearning rateと呼ばれるチューニングパラメータです．損失関数$L(Y_n,\theta)$の最適化のみを考える場合は，スケールを変えたとしても（例えば損失関数の値を100倍しても）最適な$\theta$の値は変わりません．しかし，一般化事後分布の中で損失関数を使う場合は，スケールの違いは事後分布のスケール(または不確実性)の度合いに直結するため，それを調整するのが$\omega$の役割です．妥当な$\omega$の設定方法については様々な議論があり，詳細は後編で簡単に紹介します．

負の対数尤度を損失関数にして$\omega=1$とすると，上記の事後分布は通常の事後分布に一致します．したがって，一般化事後分布は通常の事後分布の自然な拡張になっていることがわかります．

このように，対数尤度と損失関数の兼ね合いを考慮すると，上記の一般化事後分布の形は妥当そうな印象を受けますが，実際に理論的な枠組みのもとで一般化事後分布が妥当な形になっていることがBissiri et al. (2016)で示されています．こちらについても後編で紹介します．

具体例

ここでは，一般化ベイズ推測が有用であると思われる例を2つ紹介します．以下の例はBissiri et al. (2016)でも議論されているものなので，具体的な数値例などは論文をご覧ください．

生存時間解析

$y_1,\ldots,y_n$を打ち切りのある生存時間とし，$y_i$に対する説明変数（ベクトル）を$x_i$とします．このような状況での代表的な回帰モデルとして，以下のようなCox比例ハザードモデルがあります．

h(t \mid x_i)=h_0(t)\exp(x_i^\top \beta), \ \ \ \ i=1,\ldots,n.

ここで，$h(t \mid x_i)$は$x_i$が所与のもとでのハザード関数であり，$h_0(t)$はベースラインハザードと呼ばれる関数です．

応用上，回帰係数$\beta$が興味の対象になることが多く，$h_0(t)$は局外パラメータとして扱われます．古典的な頻度論的アプローチでは，$\beta$の集約尤度である部分尤度(partial likelihood)を用いた推定を行いますが，標準的なベイズ法でこのモデルを扱うためには$h_0(t)$も推定する必要があります．

$h_0(t)$はノンパラメトリックなモデル化が求められるため，計算負荷が大きくなる等の困難が生じる可能性があります．

一方で，一般化ベイズ推測の枠組みでは，部分尤度を用いて事後分布を直接的に定義することができます．$D$を観測データとすると，$\beta$の一般化事後分布は

\pi(\beta \mid D)\propto \pi(\beta)\exp\left\{-\omega \sum_{i=1}^n \ell(\beta, x_i) \right\}, \ \ \ \ \ell(\beta, x_i)=\log \left\{\frac{\exp \left(x_i^\top\beta\right)}{\sum_{j \in R_{i}} \exp \left(x_j^\top \beta\right)}\right\}

となります．ここで，$R_i$はリスク集合（時刻$y_i$においてイベントも発生しておらず打ち切りもされていない個体の集合）を表します．また，$\sum_{i=1}^n \ell(\beta, x_i)$は対数部分尤度に相当します．

このような一般化事後分布を用いることにより，局外パラメータ$h_0(t)$のモデル化を考えずに$\beta$のみのベイズ推測を行うことができます．事後分布は$\beta$の関数として特殊な形をしていますが，random-walk MHアルゴリズム等で比較的簡単に$\beta$の事後サンプルを生成できます．(Stanにおいて対数尤度を自分で定義することで事後サンプルを生成することもできると思います．）

クラスタリング

通常のベイズ推測の枠組みでクラスタリングを行うためには，データに対して混合モデル

f(x \mid C)=\sum_{j=1}^{K} p_{j} \ f_{j}\left(x \mid C_{j}\right)

を適用する必要があります．ここで，$C_j$は$j$番目の要素分布がもつパラメータ，$p_j$は$j$番目のクラスターに割り当てられる事前確率です．各要素分布$f_{j}$としては多変量正規分布等を使います．

しかし，要素分布を特定化したクラスタリングは，要素分布の特定化によって敏感に結果が変わってしまう可能性があります．データのクラスタリングのみに興味がある場合は，要素分布は局外パラメータとなり，混合モデルによるクラスタリングは冗長な方法になってしまいます．

今，データとして$x_1,\ldots,x_n$（以下ではまとめて$X_n$と表記）が得られているとします．また，データにインデックスの集合$\{1,\ldots,n\}$の任意の分割を$S$とします．（分割$S$がデータのクラスタリングを定めることに注意してください．）$S$を推定するために，以下のような損失関数を考えることができます．

L(S, X_n)=\sum_{C_k\in S}\|x_i - \bar{x}_{C_k}\|^2.

ただし，$C_k$は$k$番目のクラスターのインデックスの集合であり，$\bar{x}_{C_k}$は$k$番目のクラスターの平均である．この損失関数を使うことで，$S$の一般化事後分布は

\pi(S \mid X_n)\propto \pi(S) \exp \left\{-\omega L(S, X_n) \right\}

となります．これにより，要素分布を特定化することなく分割$S$に対する事後分布を直接的に定義することができます．

おわりに

今回は一般化ベイズ法について簡単に紹介しました．次回の記事では，Bissiri et al. (2016)で議論されている一般化事後分布の理論的妥当性やlearning rateの選択方法について紹介したいと思います．

株式会社Nospareでは統計学の様々な分野を専門とする研究者が所属しております．統計アドバイザリーやビジネスデータの分析につきましては株式会社Nospareまでお問い合わせください．