« 読了：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パッケージの使用例を紹介... 以下、機能の紹介。]
モデルの事前分布について。以下をご用意しております。
- 一様分布。bms(data, mprior="uniform")と指定する。
- 二項モデル事前分布。モデルサイズ(モデルに使う説明変数の数)について考えよう。モデルサイズの事後分布の期待値は、ある説明変数がモデルに含まれる確率(以下PIPと略記)を全変数で合計したものだ。いっぽう、モデルサイズの事前分布は、モデルの事前分布が一様だとすると、問答無用で、$K/2$を中心にした左右対称の山形になっている[←なるほど...]。モデルサイズの事前分布を思い通りに指定できないか？ そこで登場するのが「二項モデル事前分布」。モデル$M_\gamma$のモデルサイズを$k_\gamma$として、
$p(M_\gamma) = \theta^{k_\gamma} (1-\theta)^{K-k_\gamma}$
こうすると、モデルサイズの事前分布の期待値が$\bar{m}=K\theta$になる。mprior="fixed", mprior.size=$\bar{m}$とする。 - PIPの直接指定。mprior="pip", mprior.size = c(0.01, ...)という風に指定すると、最初の変数のPIPが0.01になる。
- 二項ベータモデル事前分布。二項モデル事前分布を使うと、モデルサイズの事後分布は$\bar{m}$のあたりでぐわっと高くことが知られている。そうじゃなくて、$\theta$をベータ分布からドローしようという提案がある[ややこしい...]。モデルサイズの事後分布がもっとなだらかになる。mprior="random", mprior.size=$\bar{m}$とする。
- ほかにもいろいろあるよ。
変数の数$K$が大きくて、$2^K$本のモデルを数え上げるだけで寿命が尽きる。そんなあなたのためにMCMCサンプラーをご用意しました。
- birth-deathサンプラー。まずランダムにひとつ変数を選ぶ。もしそれが現在のモデルに含まれていたら落とし、もし含まれていなかったら含める。mcmc="bd"と指定する。
- reversible-jumpサンプラー。50%の確率でbirth-deathサンプリング。残りの50%の確率で、いま含まれている変数と含まれていない変数をひとつづつランダム選びスワップ。mcmc="rev.jump"と指定する。
- MCMCじゃなくて全部数え上げる。mcmc="enumerate"と指定する。
$K$が14までならenumerate、それ以上ならbdがデフォルト。なお、MCMCのときには初期モデルをうまいこと選ぶ機能もつけてあるよ。ほかにもこんな工夫をしているよ[めんどくさいのでパス]。
パラメータの事前分布について。[このパッケージでは、どうやらg事前分布を使うというのはマストらしい]
- unit information prior。ハイパーパラメータ$g$をサンプルサイズ$N$にする。g="UIP"と指定。これがデフォルト。
- Bayesian Risk Information Criterion。$g=max(N, K^2)$とする。g="BRIC"と指定する。
- お好きなgを指定する。たとえばg=5とか。
- モデルごとに$g$を変える。経験ベイズ法(まずOLSで推定しちゃって、そのF値でgを決める)と、$g/(1+g)$がベータ分布に従うと考えてハイパーパラメータを指定する方法を搭載。[気が狂いそうなので大幅に省略。先生方、ちょっと凝りすぎじゃないですか...]
そのほか、予測分布の出し方とか、グラフィック機能の話とか。
... 読みながらパッケージをいじってたんだけど、ちょっと意外なことに気がついた。
ためしに、変数がすごく多くてサンプルサイズがすごく小さい訓練データ(N=15, K=64)でBMAを試してみたら、過学習がすさまじい。それはしょうがないんだけど、意外にも、モデルサイズの事前分布をいじって、事後分布の平均を2くらいまで下げても、やっぱり過学習するのである。つまりモデルサイズにペナルティをかけたってダメだということだ。回帰係数のハイパーパラメータgを下げれば、さすがに過学習は起きなくなるけど... 結局、素直に小さなPLS回帰モデルをひとつつくったほうが、学習性能も般化性能もよかった。
うーん、モデル・アベレージングは過学習に強いと聞いていたのだが、なぜうまくいかないのだろうか。共線性が強すぎるからだろうか、それともMCMCのやりかたがまずいのか...
論文：データ解析(2015-) - 読了：Zeugner & Feldkirchner (2015) RのBMAパッケージ