« 読了：Dahan & Mendelson (2001) コンセプト・テストの極値モデル | メイン | 読了：Zeugner & Feldkirchner (2015) RのBMAパッケージ »
2016年2月 6日 (土)
Hoeting, J.A., Madigan, D., Raftery, A.E., Volinsky, C.T. (1999) Bayesian Model Averating: A Tutorial. Statistical Science, 14(4), 382-417.
ベイジアン・モデル・アベレージング(BMA)についての解説。仕事の都合で必要に迫られ急遽目を通した。オフィスの本棚には解説書もあるはずなのに、そっちは手に取らず、題名に"A Tutorial"なんて書いてある奴をふらふらと読んでしまう、という...
BMAとはなにか。いきなりフォーマルに書いちゃうと気持ちが萎えますが、いまデータ$D$から関心ある量$\Delta$を予測するモデル$M_1, M_2, \ldots, M_K$があるとして、
$P(\Delta | D) = \sum_k P(\Delta|M_k, D) P(M_k|D)$
を求めましょうよ、という話である。
モデルの事後確率$P(M_k | D)$は、事前確率$P(M_k)$と尤度$P(D|M_k)$の積を、全モデルを通じて合計1になるようにしたものである。
尤度$P(D|M_k)$というのは積分尤度である。つまり、モデルのパラメータを$\theta$として、
$P(D|M_k) = \int P(D | \theta_k, M_k) P(\theta_k | M_k) d\theta_k$
である。[←ここですごく混乱したんだけど、一晩寝て読みなおしたら腑に落ちた。ここでいうモデルの事後確率$P(M_k | D)$とは、「推定されたモデル」が正しい確率ではなくて、なんていえばいいんだろう、モデルの式そのものが正しい確率というか、SASのPROC GLMのMODEL文が正しい確率というか、Rのlm()に与えるformulaが正しい確率というか、そういうものを指しているのだ、きっと。だから、尤度は推定されたモデルの持つ尤度$P(D|\theta_k, M_k)$ではなく、周辺尤度$P(D|M_k)$でなければならないのだ]
BMAには次の問題がある。
- A. モデルの数が多すぎて足し上げきれない。
- B. 積分の計算が難しい。MCMCが登場したいまでも大変。
- C. 事前確率$P(M_k)$をどうするか。
- D. どんなモデルを足し上げるのか。
問題Aについて。これには大きく2つのアプローチがある。
- データによって支持されるモデルだけを足しあげる。これをオッカムの窓という。それぞれのモデルについて$P(M_k|D)$を求め、その最大値と比較する。なにか閾値を決めといて、比が閾値を下回ったら、足し上げせずに捨ててしまう。さらに、「それをもっと単純にしたモデル」が生き残っていて、しかもそいつよりも確率が低い、という図々しいモデルも捨てる。こういうことをやってると、たいていモデルの数はひと桁にまで落ちる。[ああ、なるほど... こないだ読んだ機械学習ユーザ向け箴言集みたいな論文で、BMAなんてのは結局モデルを選んでるだけだ、という悪口が出てきたけど、そういうことか...]
- MCMCモデル合成(MC^3)。えーと、モデル空間$M$を状態空間とするマルコフ連鎖$\{M(t)\}$と、均衡分布$P(M_i|D)$を考えまして... [このへんで頭に霞がかかったようになってきたので中略]。柔軟なやり方だけど、収束するとは限らない。ちょっと似ている方法に確率探索変数選択(SSVS)というのもある。importance samplingを応用するのもある... 云々、云々。
問題Bについて。
- 尤度$P(D|M_k)$がほんとに積分で出せちゃう場合もある(線形回帰の場合とか)。
- 尤度$P(D|M_k)$ を積分なしで近似する方法もいくつかある。BICで近似しちゃうとか。[←ああ、そうだったのか... 各モデルで推定したなにかをBICで重み付けするのってどういう意味だろうかと思っていたのだが、あれは積分尤度の代用品なのか]
- [問題Bからはちょっと話が逸れているんだと思うけど] 事後分布$P(\delta|M_k, D)$のほうを、パラメータの最尤推定$\hat\theta$を使って$P(\Delta | M_k, \hat\theta, D)$で近似しちゃう、というのもある。これをMLE近似という。
問題Cに進む前に、ここで具体的なモデル・クラスを4つ挙げる。
- 例その１、線形回帰。予測子を選択せず、すべての可能な予測子セットを通じて平均する。尤度は閉形式で書けるぞ。俺たちの論文を読め(Raftery, et al., 1997 JASA)。反応変数をBox-Cox変換するとき、そのパラメータを決めずにBMAでどうにかする、という論文も書いたから読め。外れ値をどうこうするっても書いたから読め。[この辺、大幅中略。著者らはどうやらRのBMAパッケージの中の人らしいので、すいません、いずれ勉強いたします]
- 例その２、一般化線形モデル。リンク関数と予測子を変えてBMAする。モデルを$M_k$、ナルモデルを$M_0$として、ベイズファクターは周辺尤度の比
$B_{10} = P(D|M_k) / P(D|M_0)$
事前のオッズを$\alpha_k=P(M_k)/P(M_0)$として、モデルの事後確率は
$P(M_k|D) = \alpha_k B_{k0} / \sum_r \alpha_r B_{r0}$
じゃあ周辺尤度$P(D|M_k)$はどうやって出すのかというと... [なんか超めんどくさい話が続いているので、ええい、省略だ。そういう難しいお話は先生方にお任せしますよ！] - 例その3、生存分析。CoxモデルのBMAはどうやってやるかっているとねえ... [省略！！ 落語の「寝床」の丁稚みたいな気分になってきた]
- 例その4、グラフィカル・モデル。[だめだ、視線が紙の上をつつつーっと滑っていく... 省略]
気を取り直して、問題C、モデルの事前確率をどうやって決めるか。
- まず思いつくのは一様分布。
- 予測子セットを動かすBMAで、各予測子について「モデルにこれが含まれる」事前確率を決めておき、積をとってモデルの事前確率にするという方法がある。George & McCullock (1993 Statist. Sinica)を見よ。[←あー、これは面白いわ！ 予測子の有用性なら関係者の合議で決められそうだ。応用場面が目に浮かぶなあ...]
- グラフィカルモデルで、リンク有無に事前確率を振る。
- 専門家の合議で情報事前分布を決める。Madigan, Gavrin, & Raftery (1995 Comm.Statist.TheoryMethods), Ibrahim & Laud (1994 JASA)をみよ。
予測成績をどうやって測るか。[このくだりは次の事例1のための説明らしい。事例2では使ってない]
データを学習データ$D^B$と検証データ$D^T$に折半します。成績はこうやって測ります。検証データのすべての実現値 $d$ について
$- \sum_{d \in D^T} \log P(d | M, D^B)$
BMAそのものの予測成績は、logの右側を、それぞれのモデルにおける確率をモデルの事後確率で重み付け合計した奴に変えればよい。
[←これ、対数スコアリング・ルールじゃん。びっくり。この手法はプロパー・スコアリング・ルールになっているんだよ、なんて書いてある。ゲーム理論の文脈で出てくる考え方なんだけど、こういう場面でも使うのかー]
事例紹介をふたつ。
ひとつめ、Cox回帰の共変量を選ばずにBMAでやったという話。[難しくはなさそうだけど、面倒なので読み飛ばした]
ふたつめ、体脂肪率を回帰で予測するときにBMAでやったという話。データは$N=251$、説明変数は年齢、身長、体重など13個。全部使った重回帰で$R^2=0.75$。さて、データを折半し、学習データ側で変数選択ないしBMAをやってみた。BMAは、モデル$2^{13}=8192$個、モデルの事前確率は一様、モデルの事後確率の求め方は俺たちの論文を読め。変数選択手法としては、F値によるステップワイズ、Mallowの$C_p$最小化、調整済$R^2$最小化を試す。[←論文には「良く知られている」なんて書いてあるけど、ごめんなさい、後ろの2つがよく理解できない。$C_p$なり調整済$R^2$なりを基準に最良サブセットを悉皆検索したってこと？ 8192本のモデルから？ それって良いやり方なの？]
結果：BMAでは、モデルの事後確率は上位3位がそれぞれ12～14%、上位10位で計57%[←ほんとだ、BMAって事実上のモデル選択だ...]。変数選択のほうは、どの手法でも同一の8変数モデルが選ばれた。予測精度はBMAの勝ち。[←いやあ、悪評高きステップワイズ法なんかに勝っても、ねぇ...]
考察。
- 問題Dについて。モデルのクラス(回帰とか)をいろいろとりまぜたBMAというのもあるよ。Draper(1995)を読め。
- 頻度主義的なモデル・アベレージングもあるよ。AICで重みづけするとかね。ほかに機械学習系の、スタッキングとかバギングとかブースティングとか計算学習理論とか[←なにそれ]、いろいろあるけど、そういうのはたいてい点予測に焦点を当てていて、いっぽうBMAでは予測分布に焦点を当てているよ。
- BMAは「モデルのクラスが事前にわかっている」視点でのアプローチだと思われがちだけど、そんなことないよ。オッカムの窓をみてよ、「モデルのクラスが事前にわかってない」という発想に立ち、 いいモデルがみつかったらしょぼいのは捨てているよ。
- BMAの結果は複雑すぎて理解されないっていわれることあるけど、だったら事後確率だけみせればいいじゃん。p値より全然わかりやすいよ。なんならBMAはパラメータについての一種の感度分析だと割り切って、最良のモデルだけ見せればいいんじゃない？[←ああ、そういう思想なんだ...協調学習とは全然ちがうんだな]
- 今後の課題: 事前分布の選択、真のモデルが含まれていなかった時のパフォーマンス、オッカムの窓のような手法のチューニング、いろんなモデルクラスへの適用、など。
この論文には３人の研究者によるコメントと返事がついているので、そっちも一応目を通してみたのだけど、難しい話が多いわ、眠いわ、疲れたわ、途中でうっかり芋焼酎を飲み始めてしまったわで、ほとんど頭にはいらなかった。また次の機会にということで...
論文：データ解析(2015-) - 読了：Hoeting, et al. (1999) ベイジアン・モデル・アベレジーングへの招待