elsur.jpn.org >

« 読了:Hoeting, et al. (1999) ベイジアン・モデル・アベレジーングへの招待 | メイン | 読了:Xiao, Cao, & Zu (2015) 協調PLS回帰 (Rのenplsパッケージ) »

2016年2月 8日 (月)

Zeugner, S. & Feldkirchner, M. (2015) Bayesian Model Averaging Employing Fixed and Flexible Priors: The BMS Package for R. Journal of Statistical Software, 68(4)。
 Rのベイジアン・モデル・アベレジーング(BMA)のパッケージBMSの紹介。実戦投入前の参詣のような意味合いで目を通した。
 RのBMA用パッケージとしては、他に先週読んだ論文の著者HoetingsさんたちによるBMAパッケージとか、その論文にコメントを寄せていたClydeさんによるBASパッケージとかがあり、優劣は良くわかんないんだけど。

 BMA概論。ここでは線形回帰のBMAについてだけ考える[BMSはそういうパッケージであるようだ]。
 説明変数の候補が$K$個あり、その行列を$X$としよう。そこから$X_\gamma$を選んで、
 $y = \alpha_\gamma + X_\gamma \beta_\gamma + \epsilon, \ \ \ \epsilon \sim N(0, \sigma^2 I)$
を推定したい。でも全部入れたモデルを検討するにはサンプルサイズが足りない。
 そこでBMAだ。ありうる$2^K$個のモデルについて片っ端から推定し、それぞれにベイズ定理で事後確率を与えて平均する。モデルの事後確率(以下PMPと略記する)は
 $p(M_\gamma | y, X) = p(y | M_\gamma, X) p(M_\gamma) / \sum_s p(y|M_s, X) p(M_s)$
分母は分子を全モデルで合計しているだけである。分子のひとつめの$p(y | M_\gamma, X)$はモデルの周辺尤度 [←最大尤度とかじゃないのである。うんうん、ここはHoetingさんたちのおかげで理解できたぞ]。ふたつめの$p(M_\gamma)$はモデルの事前確率。

 統計量$\theta$(たとえば回帰係数$\beta$)について、次のようにして事後分布を得ることができる。
 $p(\theta | y, X) = \sum_\gamma p(\theta | M_\gamma, y, X) p(M_\gamma | X, y) p(M_\gamma) / \sum_s p(M_s | y, X) p(M_s)$
 [←そ、そうなの??? モデルで推定された$\theta$の事後分布$p(\theta | M_\gamma, y, X)$にモデルのPMP$p(M_\gamma | y, X)$を掛けて足しあげればいいんじゃないの? なんでここでまたモデルの事前分布$p(M_\gamma)$が出てくるの?
  探してみると、この論文のドラフトであろう BMSパッケージのVignetteでは、素直に
 $p(\theta | y, X) = \sum_\gamma p(\theta | M_\gamma, y, X) p(M_\gamma | X, y)$
となっている。これはこの論文の誤記ではなかろうか?]

 さて、モデルの周辺尤度$p(y | M_\gamma, X)$、パラメータの事後確率$p(\theta | M_\gamma, y, X)$をどう求めるか。ベイジアン回帰で求めます。[←ここがHoetingさんたちと違うところだ。彼らは個別のモデルの推定は最尤法でいいやと思ってたんじゃないかしら]
 パラメータの事前分布は、切片は$p(\alpha_\gamma) \propto 1$、誤差分散は$p(\sigma) \propto \sigma^{-1}$でよしとしよう。問題は$\beta_\gamma$の事前分布。ゼルナーのg事前分布を使うことにします。
 [...ここでg事前分布についての定義があるけど、パス。ダレダレの事前分布だなんて、そういう難しい話、できれば一生関わらずに過ごしたい...]
 まあとにかく、g事前分布を使うってことは、係数はゼロだよねと思っているということだ。ハイパーパラメータ$g$はその確信の度合いを表す(小さな$g$は「まず間違いなくゼロだよ」に相当)。
 g事前分布を使うと、$\beta_\gamma$の事後分布は$t$分布になり、その期待値は、モデル$M_\gamma$のOLS推定量を$\hat\beta_\gamma$として、$\frac{g}{1+g}\hat\beta_\gamma$になる。要するに、$0$と$\hat\beta_\gamma$を$g$で重みづけて足しているわけだ。[←へー!そりゃ面白いなあ。偉いぞゼルナーくん]。
 g事前分布を使えば、モデルの周辺尤度もすごく簡単に求められる。[数式省略]

 [ここでBMSパッケージの使用例を紹介... 以下、機能の紹介。]

 モデルの事前分布について。以下をご用意しております。

 変数の数$K$が大きくて、$2^K$本のモデルを数え上げるだけで寿命が尽きる。そんなあなたのためにMCMCサンプラーをご用意しました。

$K$が14までならenumerate、それ以上ならbdがデフォルト。なお、MCMCのときには初期モデルをうまいこと選ぶ機能もつけてあるよ。ほかにもこんな工夫をしているよ[めんどくさいのでパス]。
 
 パラメータの事前分布について。[このパッケージでは、どうやらg事前分布を使うというのはマストらしい]

 そのほか、予測分布の出し方とか、グラフィック機能の話とか。

 ... 読みながらパッケージをいじってたんだけど、ちょっと意外なことに気がついた。
 ためしに、変数がすごく多くてサンプルサイズがすごく小さい訓練データ(N=15, K=64)でBMAを試してみたら、過学習がすさまじい。それはしょうがないんだけど、意外にも、モデルサイズの事前分布をいじって、事後分布の平均を2くらいまで下げても、やっぱり過学習するのである。つまりモデルサイズにペナルティをかけたってダメだということだ。回帰係数のハイパーパラメータgを下げれば、さすがに過学習は起きなくなるけど... 結局、素直に小さなPLS回帰モデルをひとつつくったほうが、学習性能も般化性能もよかった。
 うーん、モデル・アベレージングは過学習に強いと聞いていたのだが、なぜうまくいかないのだろうか。共線性が強すぎるからだろうか、それともMCMCのやりかたがまずいのか...

論文:データ解析 - 読了:Zeugner & Feldkirchner (2015) RのBMAパッケージ