雑記: モデルを選択したい話(カステラ本7.4節~7.7節)

私の誤りは私に帰属します。お気付きの点がありましたらお手数ですがご指摘いただけますと幸いです。
テキスト(カステラ本)その他の参考文献まとめ(解釈を含む)
  • 訓練データの被説明変数には独立に同一の分布にしたがう平均 0 のノイズがのっていると考える。と、その訓練データに最適化したモデルはたまたま出たノイズに引っ張られている可能性がある。したがって、訓練データ上の誤差(訓練誤差  \overline{{\rm err}})でモデル間を比較するのは適切ではない。手元に訓練データしかないとしても、あらゆるノイズの出方に対する  \overline{{\rm err}} の期待値をとったもの(訓練標本内誤差  {\rm Err}_{\rm in})で比較するべきである。
  • が、訓練標本内誤差  {\rm Err}_{\rm in} は知り得ないので、あらゆるノイズの出方に対して最適化したときの期待値  E_y \bigl[{\rm Err}_{\rm in}\bigr] を推定することになる。例えば線形回帰モデルで2乗誤差の和を損失としたときは、 E_y \bigl[{\rm Err}_{\rm in}\bigr] = E_y \bigl[\overline{{\rm err}}\bigr] + 2d \sigma_\varepsilon^2/N になることがわかる(d は説明変数の次元数、\sigma_\varepsilon^2 はノイズの分散、N は訓練データの個数)。ので、 E_y \bigl[\overline{{\rm err}}\bigr] を計算可能な  \overline{{\rm err}} で代用した  \overline{{\rm err}} + 2d \sigma_\varepsilon^2/N {\rm Err}_{\rm in} の推定値とすればよい。2d \sigma_\varepsilon^2/N は「手元の訓練データにたまたま出たノイズを学んだせいで誤差が過小評価されている分を補完(?)するペナルティ項」といえる。
  • また、損失を交差エントロピーとしたときはより一般のモデルで漸近的に  E_y \bigl[{\rm Err}_{\rm in}\bigr] \approx E_y \bigl[\overline{{\rm err}}\bigr] + d/N が成り立つ(d はモデルの有効パラメータ数)。この右辺を  \overline{{\rm err}} で代用した  \overline{{\rm err}} + d/N赤池情報量規準AIC)とよばれる。
  • 他方、元より「モデル候補 \mathcal{M}_m \, (m=1, \cdots, M) のうち訓練データ Z に対して最も尤もらしいモデルはどれか」を出発点に、ベイズ的にモデル選択しようとすることもできる。このとき、モデル候補 \mathcal{M}_m のパラメータの最尤推定値を \hat{\theta}_m とすると \log {\rm Pr} (Z|\mathcal{M}_m) = \log {\rm Pr} (Z| \hat{\theta}_m, \mathcal{M}_m) - d_m \log N / 2 + O(1) が成り立つ(d_m はモデル候補 \mathcal{M}_m の自由度)。これに基づいた  \overline{{\rm err}} + (d \log N )/Nベイズ情報量規準(BIC)とよばれる。
  • BICAICはペナルティ項が \log N 倍異なるが、手元のモデル候補のどれが真かを見出したいときはBICが、専ら予測誤差の最小化に興味があるときはAICが適していると考えられる。
キャラクターの原作とは無関係です。
f:id:cookie-box:20180305232608p:plain:w60

モデルを選択するとき、何に注意してどう選択するべきなのかって気になりますよね。そこでカステラ本の7章が「モデルの評価と選択」であるようなんです。面倒なので261ページから読みましょう。

f:id:cookie-box:20180305231302p:plain:w60

いや俺は気になっていないんだけど。とにかくテキストの261ページね…なんか誤差の定義がやたら多くない? 訓練誤差とか汎化誤差とか期待誤差とか…。

f:id:cookie-box:20180305232608p:plain:w60

モデルの誤差には「どの訓練データで学習したのか」「その訓練データにはどんなノイズがのっていたのか」「それをどのテストデータに適用したのか」がすべて関係してきますからね。何を固定して何を固定していないかによって呼び名を変えているようです。以下にまとめてみました。

訓練誤差
 \overline{{\rm err}} = \displaystyle \frac{1}{N} \sum_{i=1}^N L\bigl( y_i, \hat{f}(x_i) \bigr)
訓練データ \mathcal{T} = \bigl\{ (x_1, y_1), \cdots, (x_N, y_N) \bigr\} で学習したモデル \hat{f} で訓練データを予測したときの予測誤差の平均値。\hat{f} は訓練データにたまたま出たノイズまで学習しがちなので、誤差の指標としては「あまりに楽観的(261ページ)」である。
訓練標本内誤差
 {\rm Err}_{\rm in} = \displaystyle \frac{1}{N} \sum_{i=1}^N E_{Y^0} \Bigl[ L\bigl( Y_i^0, \hat{f}(x_i) \bigr)  \, \Big| \, \mathcal{T} \Bigr]
訓練データ \mathcal{T} で学習したモデル \hat{f} で訓練データを予測したときの予測誤差の、各データのノイズの出方に関する期待値の、平均値。訓練データ \mathcal{T} を、たまたま出たノイズは学ばずに、きちんと学べたかという指標になる。興味があるのはこの指標ではなく汎化誤差だが、モデル間を相対比較する分にはこの指標が使用しやすい。訓練誤差を訓練標本内誤差にするための調整分としての最善度 {\rm op} \equiv {\rm Err}_{\rm in} - \overline{{\rm err}} は「ノイズを学んでしまった分」となるがこれは直接推定できない。ふつうは訓練データの被説明変数に関する期待値を取った平均最善度  \omega \equiv E_y \bigl[ {\rm op}\bigr] を推定する(「最善度」と名付けられている割に、大きいとよい指標ではなく、ゼロになるべき指標である)。
汎化誤差
 {\rm Err}_{\mathcal{T}} = E_{X^0, Y^0} \Bigl[ L\bigl( Y^0, \hat{f}(X^0) \bigr) \, \Big| \, \mathcal{T} \Bigr]
訓練データ \mathcal{T} で学習したモデル \hat{f} でテストデータ点 (X^0, Y^0) を予測したときの予測誤差の、テストデータ点 (X^0, Y^0) に関する期待値。(おそらく)この指標に最も興味があるが、訓練データ \mathcal{T} しか手元にない以上これを直接推定することはできない。
期待誤差
 {\rm Err} = E_{\mathcal{T}} E_{X^0, Y^0} \Bigl[ L\bigl( Y^0, \hat{f}(X^0) \bigr) \, \Big| \, \mathcal{T} \Bigr]
汎化誤差の、訓練データ \mathcal{T} の取り方に関する期待値。
そして、2乗損失や0/1損失など色々な損失で  \omega = (2/N) \sum_{i=1}^N {\rm Cov}(\hat{y}_i, y_i) が成り立つそうです。2乗誤差の場合でこれを示せというのが演習問題になっていますね。ハヤト、試しに示してみますか?

f:id:cookie-box:20180305231302p:plain:w60

えっ、まずその \omega が何なのかよくわかっていないんだけど…まあどのみち先に {\rm Err}_{\rm in}\overline{{\rm err}} が要るみたいだから、L を2乗誤差にすると、

\begin{split} {\rm Err}_{\rm in} &= \displaystyle \frac{1}{N} \sum_{i=1}^N E_{Y^0} \Bigl[ \bigl( \hat{f}(x_i) - Y_i^0 \bigr)^2  \, \Big| \, \mathcal{T} \Bigr] \\ &= \frac{1}{N} \sum_{i=1}^N \hat{f}(x_i)^2 + \frac{1}{N} \sum_{i=1}^N E_{Y^0} \Bigl[(Y_i^0)^2 \, \Big| \, \mathcal{T} \Bigr] - \frac{2}{N} \sum_{i=1}^N  E_{Y^0} \Bigl[ Y_i^0 \, \Big| \, \mathcal{T} \Bigr] \hat{f}(x_i) \\ &= \frac{1}{N} \sum_{i=1}^N \hat{f}(x_i)^2 + E_{Y^0} \Bigl[ (Y_i^0)^2 \, \Big| \, \mathcal{T} \Bigr] - \frac{2}{N} E_{Y^0} \Bigl[ Y_i^0 \, \Big| \, \mathcal{T} \Bigr] \sum_{i=1}^N \hat{f}(x_i) \\ \overline{{\rm err}} &= \displaystyle \frac{1}{N} \sum_{i=1}^N \bigl(\hat{f}(x_i) - y_i \bigr)^2 \\ &= \frac{1}{N} \sum_{i=1}^N \hat{f}(x_i)^2 + \frac{1}{N} \sum_{i=1}^N y_i^2 - \frac{2}{N} \sum_{i=1}^N  y_i \hat{f}(x_i) \end{split}

こうなるよな。これを {\rm op} \equiv {\rm Err}_{\rm in} - \overline{{\rm err}} に代入して、

\begin{split} {\rm op} &= E_{Y^0} \Bigl[ (Y_i^0)^2 \, \Big| \, \mathcal{T} \Bigr] - \frac{2}{N} E_{Y^0} \Bigl[ Y_i^0 \, \Big| \, \mathcal{T} \Bigr] \sum_{i=1}^N \hat{f}(x_i) - \frac{1}{N} \sum_{i=1}^N y_i^2 + \frac{2}{N} \sum_{i=1}^N  y_i \hat{f}(x_i) \end{split}

さらにこれを訓練データの被説明変数について期待値をとったのが \omega だから、訓練データのノイズが変わったら変わりうるところに E_{Y^0}[ \, \cdot \, | \, \mathcal{T}] をかぶせればいいわけだから、

\begin{split} \omega &= E_{Y^0} \Bigl[ (Y_i^0)^2 \, \Big| \, \mathcal{T} \Bigr] - \frac{2}{N} E_{Y^0} \Bigl[ Y_i^0 \, \Big| \, \mathcal{T} \Bigr] \sum_{i=1}^N E_{Y^0} \Bigl[ \hat{f}(x_i)\, \Big| \, \mathcal{T} \Bigr] \\& \quad \;- \frac{1}{N} \sum_{i=1}^N E_{Y^0} \Bigl[ (Y_i^0)^2 \, \Big| \, \mathcal{T} \Bigr] + \frac{2}{N} \sum_{i=1}^N E_{Y^0} \Bigl[ Y_i^0 \hat{f}(x_i) \, \Big| \, \mathcal{T} \Bigr] \\ &= - \frac{2}{N} E_{Y^0} \Bigl[ Y_i^0 \, \Big| \, \mathcal{T} \Bigr] \sum_{i=1}^N E_{Y^0} \Bigl[ \hat{f}(x_i)\, \Big| \, \mathcal{T} \Bigr] + \frac{2}{N} \sum_{i=1}^N E_{Y^0} \Bigl[ Y_i^0 \hat{f}(x_i) \, \Big| \, \mathcal{T} \Bigr] \\ &= \frac{2}{N} \sum_{i=1}^N \biggl\{ E_{Y^0} \Bigl[ Y_i^0 \hat{f}(x_i) \, \Big| \, \mathcal{T} \Bigr]- E_{Y^0} \Bigl[ Y_i^0 \, \Big| \, \mathcal{T} \Bigr] E_{Y^0} \Bigl[ \hat{f}(x_i)\, \Big| \, \mathcal{T} \Bigr] \biggr\} \end{split}

これって共分散の公式の形だ! だったら  \omega = (2/N) \sum_{i=1}^N {\rm Cov}(\hat{y}_i, y_i) だ(\hat{y}_i = \hat{f}(x_i) とする)。この {\rm Cov}(\cdot, \cdot) は訓練データの説明変数は固定した下での被説明変数に関する共分散(ノイズの出方が色々変わったときの共分散)だな。どう、これでいい?

f:id:cookie-box:20180305232608p:plain:w60

間違っていないとは思いますが、カステラ本には演習の解答がないのでなんとも。

f:id:cookie-box:20180305231302p:plain:w60

ええ解答ないの!?

f:id:cookie-box:20180305232608p:plain:w60

何にせよ、モデルの  {\rm Cov}(\hat{y}_i, y_i) が正である状況とは、「y_i が大きくなる方向にノイズが変化したら \hat{y}_i = \hat{f}(x_i) も大きくなる」といった状況であるわけです。そんなノイズに反応するようなモデルは望ましくありません。 {\rm Cov}(\hat{y}_i, y_i) が負であっても困ります。 {\rm Cov}(\hat{y}_i, y_i) はゼロであってほしいですよね。

f:id:cookie-box:20180305231302p:plain:w60

それはそうだな。じゃあモデルが  {\rm Cov}(\hat{y}_i, y_i) がどうなるかを調べれば、そのモデルがどれくらいノイズを学習してしまうかわかるんだな。それでさっそく263ページに加法的誤差モデルの例が出てきて…いや、 \sum_{i=1}^N {\rm Cov}(\hat{y}_i, y_i) = d \sigma_\varepsilon^2 ってなんで? というか加法的誤差モデルって何?

f:id:cookie-box:20180305232608p:plain:w60

おそらく25ページがここでいう加法的誤差モデルの定義なんじゃないかと思うんですが、25ページの記述は「加法的モデル」なんですよね…まあ、 \sum_{i=1}^N {\rm Cov}(\hat{y}_i, y_i) = d \sigma_\varepsilon^2 になること自体は線形回帰モデルの最小2乗法で確かめるといいのではないかと。257ページの上の方でもその事実をつかっていますし。少なくとも7章に導出はなさそうですが。イメージとしては、最小2乗回帰したモデルは未知の点への予測値を訓練データを重ね合わせてつくり出していると解釈できますから、それを介して訓練データのノイズを取り込んでしまうわけです。

ただ、おそらくここの主張は線形回帰より広いクラスでもこうなんだというものと思います(以下)。

  • 説明変数の各次元に適用したモデルを足し合わせるようなモデル(Ex. 線形回帰)で損失を2乗誤差(2乗誤差以外でも結構成り立つと思われる)とするときは、
    • 平均最善度(「ノイズを学んでしまった分」の期待値)は \omega = 2 d \sigma_\varepsilon^2 / N になる。
    • よって、 E_y \bigl[ {\rm Err}_{\rm in} \bigr] = E_y \bigl[ \overline{{\rm err}} \bigr] + 2 d \sigma_\varepsilon^2 / N になる(7.24式)。
    • よって、説明変数の次元 d が多いほどノイズを学んでしまいやすい。
ここまで平均最善度がはっきりしていると、「説明変数を1次元増やして(自由に調整できる要素を1つ増やして)誤差が 2 \sigma_\varepsilon^2 / N より減らないなら増やさない方がいい」などといった判断ができそうですよね。

ときにハヤト、いま「標本データ上での誤差を、自由度で調整して、モデルのよさを測る」というようなことをしたわけです。「標本データ上での誤差を、自由度で調整して、モデルのよさを測る」と聞いて、何か思い出しませんか?

f:id:cookie-box:20180305231302p:plain:w60

「自由度で調整」って、もしかして赤池情報量規準AIC)のこと? 何次多項式で回帰するべきかって文脈でよくみるよな。意味はよくわかってないけど…。

f:id:cookie-box:20180305232608p:plain:w60

はい、実際次の264ページに、より一般的なモデルに適用できる {\rm Err}_{\rm in} の推定方法として、AICが出てくるんです。まず、天下りですが N \to \infty で7.27式(以下)の関係式が成り立つことがわかっています(いつ成り立つのか曖昧ですが、AICの原論文はフリーアクセスではなさそうなので一旦置いておきます)。

 -2E\Bigl[ \log {\rm Pr}_{\hat{\theta}}(Y)\Bigr] \approx \displaystyle -\frac{2}{N} E\Bigl[ \sum_{i=1}^N \log {\rm Pr}_{\hat{\theta}}(y_i)\Bigr] + \frac{2d}{N} \tag{7.27}

f:id:cookie-box:20180305231302p:plain:w60

確かにめっちゃ降って湧いてきたな。というかこの式どういう意味?

f:id:cookie-box:20180305232608p:plain:w60

まず E[\cdot]Y の真の分布の上での期待値ですね。ここで、真の分布は {\rm Pr}_{\theta^\ast}(Y) で表されることを仮定しています。つまり、真のパラメータ  \theta^\ast さえ突き止めれば真の分布を再現できるという状況です。といっても、手元には有限の N 個のサンプルしかないのでこの理想的な状況でも真の分布の推定が必ずできるわけではありません。限られた手元の標本データで尤度が最大になるようにしたパラメータが  \hat{\theta} です。それで、7.27式の左辺は真の分布と、手元の標本データで最善を尽くした分布の交差エントロピーの2倍ですね。これは一番知りたいものです。しかし、真の分布は知りえませんから、手元の標本データの経験分布で代用すると  -(1/N) \sum_{i=1}^N \log {\rm Pr}_{\hat{\theta}}(y_i) となりますね。この値の期待値(色々標本データを取り直したときの期待値)の2倍が右辺第1項です。しかし、この右辺第1項は、 \hat{\theta} に標本データのクセが反映されている以上、左辺の推定値としては「あまりに楽観的」ですよね。その調整分が右辺第2項であるわけです。なので結局、右辺第1項の期待値を手元のデータでの値で代用した {\rm AIC} = -(2/N) \sum_{i=1}^N \log {\rm Pr}_{\hat{\theta}}(y_i) + 2d/N をモデルのよさの推定値とせよというのがここでの主張ですね。右辺第1項に出てくる \sum_{i=1}^N \log {\rm Pr}_{\hat{\theta}}(y_i) は最大対数尤度(=現在の d で尤度を最大にした \hat{\theta} での尤度)になっていますが、「手元のデータで交差エントロピーを推定したもの」といった方が意味合い的にしっくりくる気がします。同じなんですが。

f:id:cookie-box:20180305231302p:plain:w60

ふーん…あれでも、最小2乗線形回帰モデルだと  E_y \bigl[ {\rm Err}_{\rm in} \bigr] = E_y \bigl[ \overline{{\rm err}} \bigr] + 2 d \sigma_\varepsilon^2 / N だった気がするんだけど、調整分が \sigma_\varepsilon^2 だけずれていない?

f:id:cookie-box:20180305232608p:plain:w60

どちらかというと、最小2乗線形回帰モデルの  E_y \bigl[ {\rm Err}_{\rm in} \bigr] = E_y \bigl[ \overline{{\rm err}} \bigr] + 2 d \sigma_\varepsilon^2 / N が7.29式(7.27式)の両辺に \sigma_\varepsilon^2 をかけたものになっていますね。順を追って説明します。

  • まず、一般にAICを利用するときは、「真の分布からの標本データ \{x_1, \cdots, x_N\} が手元にあって、なるべく真の分布 f(x) に近い予測分布 \hat{f}(x) をつくりたい(近さの基準は交差エントロピー)」ことが大前提になるはずです。
  • 対して7.4節の文脈での当面の目的は「訓練標本内誤差 {\rm Err}_{\rm in} を小さくしたい」です。より詳しくいうと、「訓練データ \bigl\{ (x_1, y_1), \cdots, (x_N, y_N) \bigr\} の説明変数は固定した下で、被説明変数は標本なのだと考えて、真の被説明変数の分布上での期待損失が小さいモデルをつくりたい」です。真の分布 f(x,y) の予測分布 \hat{f}(x,y) がつくりたいのだと考えれば、AICを利用したい状況と同じになるでしょう。ただし、「損失が交差エントロピーであれば」です。
  • 最小2乗線形回帰モデルは損失が2乗誤差なので、AICとは目的が少し違っています。が、ノイズの分散 \sigma_\varepsilon^2 が既知であり、最小2乗線形回帰モデルの予測値を中心に分散 \sigma_\varepsilon^2 をもつガウス分布を予測分布と考えるなら、両者の目的は一致するんです。なぜなら、このときの交差エントロピー  -\sum_i \int f(x_i,y) \log \hat{f}(x_i,y) dy は2乗誤差を 2\sigma_\varepsilon^2 で割ったものに定数バイアスを足したものになるからです(※)。交差エントロピー最小化が2乗誤差最小化と一致するんです。
  • ところでAICは「交差エントロピーの2倍」の推定値として構築されていました。これをいまの状況で「2乗誤差」の推定値にかき換えたいなら、\sigma_\varepsilon^2 をかければいいです。「交差エントロピーの2倍」が2乗誤差を \sigma_\varepsilon^2 で割ったものになっているわけですから。定数バイアスは両辺で打ち消し合いますたぶん。
なので調整項に \sigma_\varepsilon^2 のずれが発生しているわけです。ちなみに(※)はこの記事の最下部でプロデューサーさんが適当に計算しただけなので適当です。

f:id:cookie-box:20180305231302p:plain:w60

あー、何を推定しているのかがずれていたのか。一瞬「AIC を使えば \sigma_\varepsilon^2 が要らないの?」って思ったけど、AIC にしたらしたで右辺第1項の計算時に \sigma_\varepsilon^2 が要るな。

f:id:cookie-box:20180305232608p:plain:w60

ノイズ \sigma_\varepsilon^2 が大きいときほどノイズを学習してしまう危険性は増すわけですからね。極端な話、ノイズが全くないデータなら {\rm Err}_{\rm in}\overline{{\rm err}} は一致して最善度は常にゼロですし。まあそれでも訓練データの隙間や外挿部分の正解は知りえませんが。

f:id:cookie-box:20180305231302p:plain:w60

ところでジュン、265ページの「もし基底関数が適応的に選ばれるなら式 (7.23) はもはや成り立たない」ってどういうこと?

f:id:cookie-box:20180305232608p:plain:w60

かいてある通りなんですが、例えば…ハヤトはいま1次元の入力で1次元の出力を回帰するモデルとして \hat{f}(x) = a_0 + a_1 x + a_2 x^2 + \cdots のような多項式を候補に考えていて、何次多項式にするかを決めたいとします。平均 \hat{f}(x)、分散 \sigma^2 であるようなガウス分布を予測分布としましょう。切片と分散も推定対象とすると、m多項式モデルのパラメータの次元数は d =m + 2 になりますね。これでAICを適用すればいいです。しかしハヤトは考えました。「AICの調整項が面倒だな…そうだ、xm 次の項までのうち一番予測に有用な1個だけを使うことにしよう。これならパラメータの次元数は常に3になるから調整項が要らない。10次くらいまで調べて単に尤度が最大になる次数を選べばいい。このアイデアはモテるだろ!」と。しかし、それでは全然駄目なんです。

f:id:cookie-box:20180305231302p:plain:w60

俺そんなこといわない! 俺のキャラクターの把握が雑! そもそもAICの調整項 2d/N の計算に面倒な要素ないだろ!! …まあともかく、なんか駄目なのはわかるよ。なんか後出しじゃんけんっぽいし。見かけのパラメータ数でAICを計算したら駄目ってことだよな。でも、それならどうやって d を決めればいいんだ?

f:id:cookie-box:20180305232608p:plain:w60

次の7.6節に出てきますが、訓練データへの予測値 \hat{y} を訓練データの被説明変数のベクトル y\hat{y} = Sy と表し、S \in \mathbb{R}^{N \times N} を「訓練データの説明変数のベクトル x には依存するが y には依存しない」行列としたとき、{\rm tr}(S) が有効パラメータ数になるということです。

f:id:cookie-box:20180305231302p:plain:w60

へ、行列のトレース?

f:id:cookie-box:20180305232608p:plain:w60

最小2乗線形回帰 \hat{f}(x) = x^\top (X^\top X)^{-1} X^\top Y なら S = X(X^\top X)^{-1} X^\top になると思います。それならトレースの公式から {\rm tr} \bigl( X(X^\top X)^{-1} X^\top \bigr) = {\rm tr} \bigl( (X^\top X)^{-1} X^\top X \bigr) = {\rm tr} \bigl( I_d \bigr) = d となって成り立っていますね。ここで dx の次元数です。先ほどの多項式の例だと S がどうなるかというのは示せないんですが、S が「y には依存しない」という条件がある以上、適応的に次数を選んだなら調べた次数を全て取り込まないといけない気がします。まずは全ての次数でパラメータを最適化してみなければならないでしょう。その上で、ある特定の次数の項だけを選ぶような強烈な正則化をするのではないかと思うんですが…。

f:id:cookie-box:20180305231302p:plain:w60

歯切れ悪いな…まあそれで次の7.7節は、{\rm BIC} = -2 \sum_{i=1}^N \log {\rm Pr}_{\hat{\theta}}(y_i) + d \log N っていうのが出てきた。これ、AICN 倍になっているというわけでもないな。調整項が 2d だったのが d \log N になっているし。でも、AICの調整項は最小2乗線形回帰の例でも裏付けられていたはずで…どうなってんの??

f:id:cookie-box:20180305232608p:plain:w60

AICとは「測っているもの」が異なるということでしょう。Schwarz の論文(BICの原論文)は以下で閲覧できました。

イントロダクション中に、AICはこうだがBICはこうだとありますね。
  • choosing the model for which \log M_j (X_1, \cdots, X_n) - k_j is largest
  • Choose the model for which \log M_j (X_1, \cdots, X_n) - \frac{1}{2} k_j \log n is largest
イントロダクションに「ある自由度のモデルの最尤推定量は、ベイズ定量のサンプル数を大きくした極限として得られる」とあり、だからベイズ的なアプローチでモデル選択できるということだと思います。というかテキスト7.7節内にも結構かいてありますね。BICの出発点は「モデル候補 \mathcal{M}_m \, (m=1, \cdots, M) のうち最もよいモデルはどれか」であるようです。そして、BICは以下が成り立つことに基づいています。

 \displaystyle \log {\rm Pr} (Z|\mathcal{M}_m) = \log {\rm Pr} (Z| \hat{\theta}_m, \mathcal{M}_m) -\frac{d_m}{2} \log N + O(1) \tag{7.40}

この式を日本語訳するならば、「モデル候補 m の下で訓練データ Z が得られる対数尤度は、そのモデルのパラメータを最尤推定 \hat{\theta}_m にした下での対数尤度から d_m \log N /2 を差し引いたものである(d_m はモデル候補 m の自由度)」といったところでしょうか。そしてこれは、AICのベースとなっている7.27式とは左辺が既に違います。あちらは「尤度最大のパラメータでの、訓練標本内交差エントロピーは?」と問うているのに対して、こちらは「どのモデルがこの訓練データに対して最も尤もらしい?」なのですから。

ところで、テキストのAICBICを比べやすくするために、右辺第1項をそろえてみましょう。

\begin{split} \frac{\rm AIC}{2} &= - \frac{1}{N} \sum_{i=1}^N \log {\rm Pr}_{\hat{\theta}}(y_i) + \frac{d}{N}\\ \frac{\rm BIC}{2N} &= - \frac{1}{N} \sum_{i=1}^N \log {\rm Pr}_{\hat{\theta}}(y_i)  + \frac{d \log N}{N} \end{split}

どちらも N \to \infty で調整項がゼロに収束するのは同じです。ただ、収束のスピードが異なります。相対的に、AICの調整項は甘く、BICの調整項は厳しいことになります。つまり、テキストにもありますが、相対的に、AICは複雑な、BICは単純なモデルを選択する傾向があるということですね。そして、「モデル選択が目的のとき、AICBICのどちらを使うべきかがはっきりしているわけではない(269ページ)」と。

f:id:cookie-box:20180305231302p:plain:w60

ええどっちがいいとかないの? そういえば、269ページに「真のモデルを含むモデル集合が与えられたときにBICが正しいモデルを選ぶ確率は、標本数 N \to \infty のとき1に近づく」(漸近一致性をもつ)ってあるけど、逆にAICはこれが成り立たないの? 成り立ってほしいような気がするんだけど…。

f:id:cookie-box:20180305232608p:plain:w60

そのようですね。つまり、真のモデルより複雑なモデルがノイズを学んでしまう分に釣り合ったペナルティを d/N では課し切れないということでしょう。…しかし、AICの目的は「真のモデルはどれだろうか」だったでしょうか? 確かにBICの目的は「真のモデルはどれだろうか」でした。しかし、AICの目的は「そのモデルを尤度最大にチューニングした下で真の分布とどれだけ離れているか」なんです。真のモデルを外していようと、誤差を小さくすれば正義なんです。

例えば以下の文献にはAICは漸近有効性をもつとあります。つまり、予測誤差を最小にするのはAICであるということです。

そしてこの文献は、AICBICは何を最適としているかが違うと指摘しています。
  • In fact the conflict is easily resolved once it is acknowledged that ‘‘asymptotically optimal’’ can have several meanings. Asymptotic efficiency and (asymptotic) consistency are different kinds of optimality.(633ページ)
この文献によると、McQuarrie and Tsai (1998) は、AICを筆頭とする予測誤差を最小にするための指標を「efficient estimators」、BICを筆頭とするモデルを確認するための指標を「consistent estimators」と区別したそうです。AICの眷属にはテキスト7.5節にも出てきた Mallows’ C_p や交差検証(cross validation)が、BICの眷属には Hannan and Quinn や Geweke and Meese があるようですね。

f:id:cookie-box:20180305231302p:plain:w60

えっと、モデル選択の指標が2つの派閥(?)に分かれるのとか、両者はなんか目的が違うっぽいのはわかったけど、そこまでいうならどっちを使うべきってないの?

f:id:cookie-box:20180305232608p:plain:w60

先の文献には「真のモデルが非常に複雑で(自由度が非常に大きく)真のモデルを当てられる見込みがなく、モデル候補の自由度をそれに合わせなくてもよく、モデル候補のパラメータが真のモデルのパラメータを包含していなくてもよいし、余分なパラメータを含んでいてもよい」というときは前者を、「真のモデルがシンプルで、モデル候補のうち1つが真のモデルであると期待されるとき」は後者を利用するべきだとあります。確かに前者の状況でBICを使う意味は薄いと思うんですよね。モデル候補が真のモデルを含んでいる見込みが薄いなら。でも、ではこの状況でAICは上手く機能するのかはきちんと導出を追っていないのでわからないんですよね。いいんだと思うんですが。

おわり

※ 交差エントロピーが2乗誤差を 2\sigma_\varepsilon^2 で割ったものに定数バイアスを足したものになっている。
 \begin{split} \displaystyle -\int f(x_i,y) \log \hat{f}(x_i,y) dy &=-\frac{1}{\sqrt{2 \pi \sigma_\varepsilon^2}} \int  \exp \left( -\frac{(y - y^\ast)^2}{2 \sigma_\varepsilon^2} \right) \log \left[ \frac{1}{\sqrt{2 \pi \sigma_\varepsilon^2}} \exp \left( -\frac{(y - \hat{y})^2}{2 \sigma_\varepsilon^2} \right) \right] dy \\ &=- \frac{1}{\sqrt{2 \pi \sigma_\varepsilon^2}} \int   \left( -\frac{(y - \hat{y})^2}{2 \sigma_\varepsilon^2} \right) \exp \left( -\frac{(y - y^\ast)^2}{2 \sigma_\varepsilon^2} \right)dy + \log \sqrt{2 \pi \sigma_\varepsilon^2} \\ &= -\frac{1}{\sqrt{2 \pi \sigma_\varepsilon^2}} \int   \left( -\frac{(y - y^\ast)^2 +(y^\ast -\hat{y})^2 - 2(y - y^\ast)(y^\ast -\hat{y})}{2 \sigma_\varepsilon^2} \right) \exp \left( -\frac{(y - y^\ast)^2}{2 \sigma_\varepsilon^2} \right)dy  + \log \sqrt{2 \pi \sigma_\varepsilon^2} \\ &= \frac{\sigma_\varepsilon^2}{2 \sigma_\varepsilon^2} + \frac{(y^\ast -\hat{y})^2}{2 \sigma_\varepsilon^2} + \log \sqrt{2 \pi \sigma_\varepsilon^2} \end{split}