雑記：モデルをアンサンブルしたい話（その2―カステラ本7.11節、8.2節、8.4節、8.8節、10.1～10.4節）

私の誤りは私に帰属します。お気付きの点がありましたらお手数ですがご指摘いただけますと幸いです。

テキスト（カステラ本）

統計的学習の基礎 ―データマイニング・推論・予測― | Trevor Hastie, Robert Tibshirani, Jerome Friedman, 杉山将, 井手剛, 神嶌敏弘, 栗田多喜夫, 前田英作, 井尻善久, 岩田具治, 金森敬文, 兼村厚範, 烏山昌幸, 河原吉伸, 木村昭悟, 小西嘉典, 酒井智弥, 鈴木大慈, 竹内一郎, 玉木徹, 出口大輔, 冨岡亮太, 波部斉, 前田新一, 持橋大地, 山田誠 |本 | 通販 | Amazon の 256～257, 326～327ページ。

関連記事

まとめ（解釈を含む）

並行世界の自分が学習したモデルを取り寄せればバリアンスを小さくできる。が、そんなことはできないので自分で並行世界をつくり出す方法がブートストラップ法である（そのモデルたちを平均するのがバギングである）。
カテゴリ値を取るデータの標本 $Z$ からのノンパラメトリックブートストラップ標本は、「標本 $Z$ に対して、 $Z$ 内の各カテゴリの割合を生成したディリクレ分布をベイズ推定して、各カテゴリの割合がその事後分布に（ほぼ）したがうように標本をサンプリングする」ことをしたものであるといえる。ことから、バギングしたモデルは近似的にモデルのベイズ事後分布の平均であるとみなせる（ので普通に訓練したモデル＝MAP推定より2乗誤差が小さくなることが見込める）。
もっともベイズ事後平均を取るのであれば、異なるモデルたちを候補にしてもよい。そのときモデルを足し合わせる重みはBIC（「どのモデルがこの訓練データに対して最も尤もらしいか」という指標だった）を用いてもよいし、直接最適な重みを推定してもよい（スタッキング）。
そもそも足し合わせる各モデルが同じ標本を学ばなくてもよい。1つ前のモデルが誤分類したデータを重点的に学んでいくのがブースティング法の最も基本的なアルゴリズムであるアダブーストM1である。アダブーストM1は前向き段階的加法的モデリングにおいて指数損失を採用しても導出される。というより、前向き段階的加法的モデリングを解く手法全般をブースティングという。

キャラクターの原作とは無関係です。

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

前回は「並行世界の自分が学習したモデルを取り寄せればバリアンスを小さくできる。が、そんなことはできないので自分で並行世界をつくり出す方法がブートストラップ法である」というところまででしたね。

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

だいぶブートストラップ法に対する語弊がないかそれ…。そもそもブートストラップ法って何？自分で並行世界をつくるなんてことしていいの？

テキストでブートストラップ法が最初に出てくるのは286ページの7.11節ですね。前回もいった通りここでのブートストラップ標本とは「訓練データ $Z = \bigl\{ (x_1, y_1), \cdots, (x_N, y_N)\bigr\}$ から $N$ 点を復元抽出する」ことを $B$ 回繰り返して $B$ セットの新しい標本 $Z^{\ast 1}, Z^{\ast 2}, \cdots, Z^{\ast B}$ を得たものです。復元抽出なので元の訓練データの点はある $Z^{\ast b}$ に複数回抽出されていることもありますし、1回も抽出されていないこともあります。カギカッコ内は「訓練データの経験分布から $N$ 点生成する」といっても同じですね。なお、これはあくまでノンパラメトリックブートストラップといわれる方法であって、ブートストラップ標本を得る方法は他にも色々あります（テキスト302ページのノイズを添加する方法など）。より一般に、訓練データ自体を用いて並行世界たちの訓練データを複製する方法であればブートストラップ法のようです。そもそもブートストラップとは靴の後部の小さなつまみのことで、"pull oneself up by one's bootstraps" という慣用句で「他人の助けを借りずに自力で達成する」という意味であるようです。ので、この訓練データだけから新たな標本を生成する方法がブートストラップ法と名付けられたのでしょう。初出は1979年とか。

そういう手法なので、期待値や分散や信頼区間のモンテカルロ推定に利用されるのがメインなのではないでしょうか。7.11節に出てくるブートストラップ法の用法は以下ですね。汎化誤差の期待値の推定は込み入っています。ブートストラップ標本は元の標本の点を全ては含んでいませんから、そのディスアドバンテージ分をどう補ってやるかということになってくるようですね。

訓練データから計算される何らかの量 $S(Z)$ の分散を推定するのに、 $B$ 個のブートストラップ標本上の不偏分散でもって推定値とする。
汎化誤差を推定するのに、 $B$ 個のブートストラップ標本で学習したモデル $\hat{f}^{\ast 1}, \cdots, \hat{f}^{\ast B}$ による平均予測誤差の、全訓練データの平均を推定値 $\widehat{\rm Err}^{(1)}$ とする。なお、各 $x_i$ の平均予測誤差を出す際、点 $x_i$ を訓練データとして参照したモデルはとばす。また、この推定値には亜種がある（以下）。

上の推定値では、「個々の $\hat{f}^{\ast 1}, \cdots, \hat{f}^{\ast B}$ は本来の約63.2%のデータしか学べていない（ので汎化誤差を過大評価しているかもしれない）」という問題点がある。これを緩和するために、訓練誤差を36.8%混ぜる（このようになる導出はあるらしい）。これを $\widehat{\rm Err}^{(0.632)}$ とする。
上の推定値でも、「過学習のために訓練誤差が小さくなりすぎているかもしれない（ので訓練誤差を混ぜると逆に過小評価になるかもしれない）」という問題点がある。そこで、「相対過学習率」を、「 $\widehat{\rm Err}^{(1)}$ から訓練誤差を引いたもの」を「でたらめな正解ラベルを付けたときのそのモデルの予測誤差から訓練誤差を引いたもの」で割ったものとして定義する。「相対過学習率」が大きいほど $\widehat{\rm Err}^{(1)}$ を信用する。これを $\widehat{\rm Err}^{(0.632+)}$ とする。

ブートストラップ標本は元の標本の点を全ては含んでいないというのに関連して、元の標本が $N$ 点であったときにあるブートストラップ標本にある点が選ばれる確率は以下ですね。

def arutenga_erabareru_kakuritsu(N):
    p = 1.0 - 1.0 / float(N)  # 1回の抽出で自分以外の点が選ばれる
    p_total = 1.0  # N回の抽出すべてで自分以外の点が選ばれる
    for i in range(N):
        p_total *= p
    return 1.0 - p_total

ある点が選ばれない確率は $N = 1$ で $0$ 、 $N = 2$ で $0.25$ で、 $N \to \infty$ の極限でネイピア数の逆数になります。 $N' = -N$ とおけば $(1-1/N)^N = [ (1+1/{N'})^{-N'} ]^{-1} \xrightarrow[-N' \to - \infty ]{} e^{-1}$ ですから。
ただ7.11節にブートストラップ標本によって推定することがどのように正当かという記述はありませんね。そのような話題は8章にあるようです。以下にまとめます。文字の定義はテキストを参照してください。

真の値のノイズが加法的でガウシアンなときに、パラメトリックブートストラップ標本で学習したモデルの期待値が最尤推定の結果と一致する（8.2節）。
ガウス分布 $N(\theta, 1)$ の期待値 $\theta$ をベイズ推定する状況（ $\theta$ の事前分布も事後分布もガウス分布）では、事前分布の分散を無限大にすると（無情報事前分布）、ある点 $z$ を観測したときの事後分布が $N(z,1)$ となるが、パラメトリックブートストラップ標本はこの事後分布から1点 $\theta'$ をサンプリングして $N(\theta', 1)$ を得てそれの期待値を採用したものとみなせる（8.4節）。

つまり、パラメトリックブートストラップ標本は、「各点 $z$ に対して、 $z$ を生み出したのは $N(z,1)$ なんだろうなとベイズ推測して、その分布から点をサンプリングする」ことを各点に対してしたものであるといえる。

カテゴリ数が $L$ のカテゴリカルなデータの標本を得たとする。各カテゴリが生成される割合 $w$ の事前分布をすべてのカテゴリの集中度母数が0の極限のディリクレ分布とすると、事後分布はすべてのカテゴリの集中度母数が標本内での観測数であるディリクレ分布となるが、ノンパラメトリックブートストラップ標本内の各カテゴリの割合がしたがう分布（標本内の観測割合を母数にもつ多項分布）はこれに近い（平均が等しく分散共分散がほぼ等しい）（8.4節）。

つまり、ノンパラメトリックブートストラップ標本は、「標本 $Z$ に対して、 $Z$ 内の各カテゴリの割合を生み出したのは ${\rm Dir} (N \hat{w})$ なんだろうなとベイズ推測して、各カテゴリの割合がその分布にほぼ近くなる方法で標本をサンプリングする」ことをしたものであるといえる。

文字に色を付けた部分の気持ちになれば、天下りでなく演繹的にブートストラップ法にたどり着けるような気がするかもしれません。僕はしませんが。

しないのかよ！？あと、「ほぼ等しい」ってなんか中途半端じゃない？どれくらい等しいの？？

テキストにそうかいてあるものはかいてありますから…まあ確かめますか。 ${\rm Dir} (N \hat{w})$ にしたがう $\tilde{w}$ の分散は $\hat{w}_i (1 - \hat{w}_i) / (N + 1)$ で、 $N$ 倍が ${\rm Multi}(N, \hat{w})$ にしたがう $\tilde{w}$ の分散は $\hat{w}_i (1 - \hat{w}_i) / N$ ですから近いですね。共分散は $-\hat{w}_i \hat{w}_j / (N + 1)$ と $-\hat{w}_i \hat{w}_j / N$ なので近いです。

確かに近いな。じゃあこの状況でのブートストラップ法はほぼベイズ推定なんだ。

はい、テキスト8.8節「モデルの平均と統合」の冒頭にも「バギングしたモデルは近似的なベイズ事後平均になっていることがわかる」のような文があります。対して、普通に訓練したモデルはベイズ事後分布の最頻値をとっていることに対応すると。そして、2乗誤差に対して最適なのはベイズ事後平均です。ということからも、バギングすることで2乗誤差を小さくできることが見込まれます。これが一応の「自分で並行世界をつくるなんてことをしていいのか」に対する答えになるのではないかと。

そういわれるとそうなのか。

さらに8.8節では「平均する個々のモデルを異なるモデルにしたらどうだろうか」という議論を巡らせています。どう思いますか？

えっ、そんな、別々のモデルにしちゃったら、モデルたちはベイズ事後分布をなさないよな…「並行世界の自分にモデルをもらったけど並行世界の自分が学習したのは違うモデルだった」ってなったら、自分のモデルとどっちが尤もらしいかってわからないわけだし…それじゃ平均したら駄目なんじゃないの？

ええ、どちらが尤もらしいかわかりませんね。でも、前々回に、「どのモデルがこの訓練データに対して最も尤もらしい？」という指標を取り扱いませんでしたか？

あ、BIC…。

はい、そのような場合は BIC が利用できるとあります。また、より直接的に、「どのモデルにどれだけ重みを割り振れば誤差が最も小さくなる？」と重みを最適化するのが331ページ下部で紹介されているスタッキングです。ここで、モデル $m$ の $i$ 番目のデータに対する誤差の計算時には $i$ 番目のデータのみを除いた訓練データで学習したモデル $m$ を用意し、それで予測します。明確に記述されていませんが、データ数×モデル数だけ学習することになるので、結構面倒そうな気はします。

それで、バギングの重要例であるランダムフォレストは15章で掘り下げられています（ただ、ランダムフォレストでは決定木の分岐点を決めるときに特徴量もランダムに絞り込んでさらにバリアンスを下げる工夫をしているのでバギングというのかバギングの改造版というのかわかりませんが）。

さておき、ここまでで、「モデルを平均するとよい」「異なるモデルを異なる重みで平均してもよい」というところまできたわけです。であれば、その重みは入力空間内の点によっても異なっていいと思いませんか？いい換えると、ある場所ではモデル1が重く、ある場所ではモデル2が重いといったように、あたかも得意不得意で役割分担するイメージです。

その点で予測誤差が小さいモデルを重くするってこと？さすがにそれは後出しじゃんけんじゃない？というか、その点のノイズに引っ張られまくる気が。

いえ、ブートストラップ標本で学習したモデルでそれをやるということではないです。あるモデルはある場所から重点的にサンプリングした訓練データで学習し、また別のあるモデルは別の場所から重点的にサンプリングした訓練データで学習し…といったイメージです。それも入力空間を適当に区切るわけではなく、まだ上手く学習できていない箇所に重点を置きにいくんです。もはや個々のモデルの訓練データが異なるということです。…もうネタばらししてしまうんですが、10章のブースティング法がその手法です。その最も基本的なアルゴリズムである、アダブーストM1（10.1節）が以下です。ゴールは2クラス分類問題を解くことで、各分類器 $G_m(x)$ は $-1$ か $1$ を出力するという状況です。

各データの重み $\{w_i\}_{i=1}^N$ を均等に初期化する。
重み $\{w_i\}_{i=1}^N$ で分類器 $G_1(x)$ を学習する。
$G_1(x)$ の重み付き誤分類率 ${\rm err}_1$ を計算し、正解しやすさの対数オッズ $\alpha_1 = \log \bigl( (1 - {\rm err}_1)/{\rm err}_1 \bigr)$ を計算する。
各データの重みを更新する。 $G_1(x)$ が正しく予測できたデータの重みは更新しない。 $G_1(x)$ が正しく予測できなかったデータの重みには $e^{\alpha_1}$ をかける。
2. に戻る（更新した重み $\{w_i\}_{i=1}^N$ で分類器 $G_2(x)$ を学習する）。
最終的な予測モデルを $G(x) = {\rm sign} \bigl[ \sum_{m=1}^M \alpha_m G_m(x) \bigr]$ とする。

「正しく予測できなかったデータの重みを増幅しておく（ことで次以降のモデルに託す）」ことをやっているのがわかると思います。また、重みの増幅率やモデルを足し合わせる係数に対数オッズが出てくるのにも理由があります。

強引に10章にとんだな…。

ブースティングの話に触れておきたかったので。モデルのアンサンブルという文脈でバギングとブースティングが対比される場合が多いと思うんです。というよりは、ランダムフォレストと勾配ブースティング木が対比されるんでしょうか…しかし、テキスト10章の冒頭にもあるように両者は「本質的に異なる手法」なんです。どちらも決定木がたくさんあるくらいで。なので話を上手く10章にもっていけませんでした。なんていうか、一度「モデルを組み合わせよう」という気持ちを捨てた方がいい気がしますね。

うん、いきなりアダブーストM1っていわれてもブートストラップ標本が出てきたときより突拍子がないんだけど、ブートストラップ標本がベイズ推定から導出されるなら、そのアダブーストM1は何から導出されるの。

それが10.3節の「前向き段階的加法的モデリング」ですね。いま学習したいモデルが10.3式のように他のモデルの線形和でかけるとします。便宜上、足し合わせるモデルを弱学習器とよび、足し合わせたモデルを強学習器とよびます。10.3式のようなモデルは全てのパラメータを一気に最適化すべきですが、計算コストなどの問題でそれが上手くできないとき、弱学習器を1つずつ最適化できないかということになります。それが前向き段階的加法的モデリングに他なりません。391ページのアルゴリズム10.2です。

391ページのアルゴリズム10.2ね…って、387ページのアダブーストM1とだいぶ違わない？重みとか出てこないし。それに、アダブーストM1はステップ $m$ で学習した弱学習器の損失をみていたけど、この式だとステップ $m$ までの強学習器の損失って感じだし。

いえ、アルゴリズム10.2で $L$ を指数損失としてみましょう。そうする各ステップですべき最適化が10.9式になります。「ステップ $m$ までの強学習器の損失」が「重み」と「ステップ $m$ で学習した弱学習器の損失」の積になっていますよね。むしろ「重み付き誤分類率」とは「ステップ $m$ までの強学習器の損失」だったんです。だから重み付き誤分類率を最小化する弱学習器 $G_m(x)$ を学習すればいいんです。

なるほど。あれでも、弱学習器 $G_m(x)$ をどれだけの重みで強学習器に足せばいいんだ…って、10.12式か。え、これどうやって出すの？

演習10.1になっていますね…そうですね、以下で出ます。

$\sum_{i=1}^N w_i^{(m)} \exp \bigl( - \beta_m y_i G_m(x_i) \bigr)$ が最も小さくなればいいんですよね。
$\sum_{i \in {\rm OK}_m} w_i^{(m)} \exp \bigl( - \beta_m \bigr) + \sum_{i \in {\rm NG}_m} w_i^{(m)} \exp \bigl( \beta_m \bigr)$ と書き換えます。 ${\rm OK}_m$ は $G_m(x_i) = y_i$ であるインデックスの集合、 ${\rm NG}_m$ は $G_m(x_i) \neq y_i$ であるインデックスの集合とします。
$\beta_m$ についての微分がゼロなら $- \sum_{i \in {\rm OK}_m} w_i^{(m)} \exp \bigl( - \beta_m \bigr) + \sum_{i \in {\rm NG}_m} w_i^{(m)} \exp \bigl( \beta_m \bigr) = 0$ です。
よって $\exp \bigl( 2 \beta_m \bigr) \sum_{i \in {\rm NG}_m} w_i^{(m)} = \sum_{i \in {\rm OK}_m} w_i^{(m)}$ です。
よって $\displaystyle \beta_m =\frac{1}{2} \log \frac{ \sum_{i \in {\rm OK}_m} w_i^{(m)}}{\sum_{i \in {\rm NG}_m} w_i^{(m)}} =\frac{1}{2} \log \frac{ 1 - {\rm err}_m}{{\rm err}_m}$ です。

（その3があれば）つづく