ベイズ・ファクターのBICによる近似についてなんだかよくわかんなくなっちゃったので、ちょっとメモしておく。
出典は
Wagenmakers, E.J. (2007) A practical solution to the pervasive problems of p values. Psychonomic Bulletin & Review, 14(5), 779-804.
この論文は前に目を通したんだけど(なんと4年前だ)、長い論文だし、途中で飽きてしまって適当に読み飛ばしていた。
いわく…
BICの差をベイズ・ファクターの近似に変換できる。さらに、いま検討しているモデルたちの事前の確からしさが等しいという仮定の下で、BICの比較からモデルの事後確率の近似を得ることができる。たとえば、帰無仮説を対立仮説と比較して、\(BIC(H_0) = 343.46, BIC(H_1)=341.58\)だったとしよう。これは\(Pr_{BIC}(H_0|D) \approx 0.28, Pr_{BIC}(H_1|D) \approx 0.72\)と変換できる。
なぜかというと…
所与のモデルなり仮説なりを\(i\)として、その事前予測確率[←事前分布から得られるデータの周辺分布のことね]は $$ Pr(D|H_i) = \int Pr(D|\theta_i, H_i) Pr(\theta_i|H_i) d \theta_i$$ ただし\(D\)は観察データ, \(\theta_i\)はモデル\(H_i\)のパラメータベクトルである。
導出は省くけど、ここから次の近似が得られる。$$ \log Pr(D|H_i) = \log Pr(D|\hat{\theta}_i, H_i) – (k_i / 2) \log n + O(1)$$ ただし、\(\hat{\theta}_i\)は\(\theta_i\)の最尤推定値、\(k_i\)はモデル\(H_1\)のパラメータ数、\(n\)は観察の数。\(O(1)\)というのは、\(n \rightarrow \infty\)でも近似誤差が0に近づかないということを示しているが、\(\log Pr(D|H_i)\)に占める割合としてみれば0に近づく。
ある種の事前分布を使えば近似誤差は\(O(n^{-1/2})\)に減る(\(n \rightarrow \infty\)で近似誤差が0に近づく)。Jeffrey事前分布とか、単位情報事前分布とかがそうだ。単位情報事前分布というのは、ひとつの観察が含んでいるのと平均して同じ量の情報を含んでいる事前分布である。たとえば、\(x^n = (x_1, x_2, \ldots, x_n) \sim N(\mu, 1)\)について\(H_0: \mu = \mu_0, H_1:\mu \neq \mu_0\)の検定をするとき、単位情報事前分布は、平均は標本平均でSDは1の正規分布である。
近似の良さを決めるのは、最大尤度\(Pr(D|\hat{\theta})\)、そしてモデルの複雑さ\((k/2) \log n\)である。罰則項に自由パラメータ数だけではなくて標本サイズがはいっている点に注意。
さて、上の式を-2倍するとBICが得られる。つまり、\(L_i=c Pr(D | \hat{\theta}, H_i) \)(\(c\)は定数)として、$$ BIC(H_i) = -2 \log L_i + k_i \log n$$ である。逆に言うと$$ Pr_{BIC}(D|H_i) = \exp\left( -\frac{BIC(H_i)}{2} \right) $$である。
ご存じのように、$$ \frac{Pr(H_0|D)}{Pr(H_1|D)} = \frac{Pr(D|H_0)}{Pr(D|H_1)} \frac{Pr(H_0)}{Pr(H_1)} $$の右辺第一項をベイズ・ファクターという。ここから、$$ BF_{01} \approx \frac{Pr_{BIC}(D|H_0)}{Pr_{BIC}(D|H_1)} = \exp\left( \frac{BIC(H_1)-BIC(H_0)}{2} \right) $$ である。
話を戻して…
いま、\(BIC(H_0) = 1211.0, BIC(H_1) = 1216.4\)だったとしよう。\(BF_{01} = \exp(5.4/2) \approx 14.9\)である。モデルの事前確率が等しいとすれば、\(H_0\)の事後確率は[\( Pr(H_0|D) \approx 14.9 (1-Pr(H_0|D)) \) より] \(14.9 / 15.9 \approx 0.94\)である。
類似した実験で\(BIC(H_0) = 1532.4, BIC(H_1) = 1534.2\)であったとしよう。今度は\(BF_{01} = \exp(1.8/2) \approx 2.5\)である。2つの実験をあわせると\(BF_{01}(Total) = 14.9 \times 2.5 = 37.25\)である。\(Pr(H_0|D) = 37.25/38.25 \approx 0.97 \)となる。[なるほど]
モデルが\(k\)個ある場合は$$ Pr_{BIC}(H_i|D) = \frac{\exp(-\frac{1}{2}BIC(H_i))}{\sum_{j=0}^{k-1} \exp(-\frac{1}{2}BIC(H_j))} $$ モデルが2個のとき、\(H_0\)の事後確率は $$ Pr_{BIC}(H_0|D) = \frac{1}{1+ \exp(-\frac{1}{2} (BIC(H_1) – BIC(H_0))} $$ つまり帰無仮説の事後確率はBIC差の0.5倍のロジスティック関数である。BIC差が12で事後確率は0.998, とても強い証拠となる。なお、この式は\(H_0\)が\(H_1\)にネストされていようがいまいが成り立つという点に注意。
ベイジアン仮説検定をBICで近似する方法には次の2つの美点がある。
- 事前分布を指定しなくて良い。BICは単位情報事前分布を暗黙に仮定しており、主観がはいらない。もっともこれは欠点でもあって、単位情報事前分布はあまりに広すぎて帰無仮説が有利になりがちだという面もあるし、BICで近似していると実質的情報を統合することができないという見方もできる。
- 計算が楽。
いっぽう限界は、モデルパラメータの関数形を無視している点。パラメータ数が同じモデルでも複雑さは異なることがある。たとえば、精神物理学でいうStevens法則は減速関数と加速関数の両方を扱えるが、Fechner法則は減速法則しか扱えない。でもパラメータ数は同じである。本物のベイジアン分析はパラメータの関数形に敏感だ(パラメータ空間全体を通じて尤度を平均するから)。仮にデータが減速関数に従っていたら、Stevens法則の事前予測分布は、パラメータ空間が加速関数をカバーしているぶんだけ損なわれるわけだ。この関数形の問題は、特に複雑な非線形モデルの場合に重要である。
なお、観察数\(n\)は場合によっては自明でない。たとえば、被験者10人が2条件のそれぞれでデータを提供した場合とか。こういう場合の標準的な選択は\(n\)を被験者数とみなすことである。
最後に、BICとANOVAの関係についてご紹介しよう。実験心理学では正規線型モデルを使うことが多いからね。
誤差の正規性を仮定した線形回帰の場合、BICはこう書くことができる。$$ BIC(H_1) = n \log (1-R^2_i) + k_i \log n$$ ANOVAでいえば、\(H_i\)の誤差平方和を\(SSE_i\)として、\(1-R^2_i = SSE_i / SS_{total}\)ですわね。
\(H_1\)と\(H_0\)を比較する分には、\(SS_{total}\)は尤度比から消える。BIC差はとっても簡単で $$ BIC(H_1) – BIC(H_0) = n \log \left( \frac{SSE_1}{SSE_0} \right) + (k_1 – k_0) \log n $$ となる。
云々。
なるほどねえ…特に、関数形の話が勉強になった。
ときに、たとえば\(y = \alpha + \beta x + e\)というモデルは所与として、よくある\(H_0: \beta = 0\)と\(H_1: \beta \neq 0\)の比較じゃなくて、\(H_0: \beta \leq 0\)と\(H_1: \beta > 0\)を比較したい、ってことがあるじゃないですか。実務的にはすごくニーズがある。
こういう線型回帰なら、\(\beta\)の事前分布をうまく設定すれば事後分布を解析的に求められるだろうけど、もっと複雑なモデルだとそうはいかない。じゃあサンプリングしろって云われそうだけど、そんな面倒なことはやってらんないことが多い。
こういうとき、事後確率\(P(H_1|D)\)をBICで近似することはできるのだろうか？ 仮説の事前確率をどうにかして決められるとして(\(P(H_1)=P(H_0)\)とか？ そこも不思議な感じだけど、それはまあ置いといて)、\(y = \alpha + \beta x + e, \beta \leq 0\)というモデルと\(y = \alpha + \beta x + e, \beta > 0\)というモデルをそれぞれ推定してBICを求めていいわけ？ うーん、きっとどっちかの\(\hat{\beta}\)は0の近くに張り付くよね… なんか直感的に変な感じがするなあ… 気のせいかなあ…