統計学 学習ノート

統計検定の4級から1級まで、出題範囲の順に記事と練習問題で学ぶ

1級数理-13 統計数理

モデル評価基準

  • カルバック・ライブラー情報量
  • 情報量規準AIC
  • クロスバリデーション

はじめに

「AICの罰則は、なんで『2 × 母数の個数』なの?」

準1級では、AICの式を使ってモデルを選びました。過去問を丸暗記した受験生の話を覚えていますか? 手元のデータの成績はいつも甘く出るので、本番の成績を見積もるためにハンデを足す、という話でした。

でも、そこで気になりませんでしたか? ハンデが「母数の個数の2倍」なのは、たまたまなのでしょうか。1.5倍でも3倍でもなく、ぴったり2倍である理由は何でしょうか?

実は、この2は適当に決めた数ではありません。手元の成績がどれだけ甘いかを、期待値として計算した結果です。計算は、テイラー展開と行列のトレースでできています。

そこで登場するのが、今回の主役 情報量規準AIC です。AICの土台である カルバック・ライブラー情報量 と、計算でなくデータの分割で本番を再現する クロスバリデーション も扱います。

前回の記事では、1つのモデルの中で推定量のよさを立証しました。今回の記事では、モデルどうしのよさを測る基準として、カルバック・ライブラー情報量の性質、AICの導出、クロスバリデーションの理論について解説します。

カルバック・ライブラー情報量:性質を証明する

はじめに、モデルのよさの物差しであるカルバック・ライブラー情報量について解説します。

真の分布の密度を gg、モデルの密度を ff とします。カルバック・ライブラー情報量(KL情報量)は次の量でした。

D(g ∥ f)=Eg[log⁡g(Z)f(Z)]=∫g(z)log⁡g(z)f(z) dzD(g \,\|\, f) = E_g\left[\log\frac{g(Z)}{f(Z)}\right] = \int g(z)\log\frac{g(z)}{f(z)}\,dz

準1級では、非負で非対称な「分布のズレ」として使いました。論述では、非負性を証明できるようにしましょう。

【非負性の証明】

  1. すべての x>0x > 0 で log⁡x≤x−1\log x \le x - 1 が成り立ち、等号は x=1x = 1 のときだけである
  2. x=f(z)/g(z)x = f(z)/g(z) を代入して gg の下で期待値をとると、次の式を得る
−D(g ∥ f)=Eg[log⁡f(Z)g(Z)]≤Eg[f(Z)g(Z)−1]=∫g>0f(z) dz−1≤0-D(g \,\|\, f) = E_g\left[\log\frac{f(Z)}{g(Z)}\right] \le E_g\left[\frac{f(Z)}{g(Z)} - 1\right] = \int_{g > 0} f(z)\,dz - 1 \le 0
  1. 等号は、確率1で f(Z)/g(Z)=1f(Z)/g(Z) = 1 のときだけ、すなわち f=gf = g のときだけである

イェンセンの不等式を使っても同じですが、log⁡x≤x−1\log x \le x - 1 なら1行で済みます。

次に計算の練習です。2つの正規分布 g=N(μ1,σ12)g = N(\mu_1, \sigma_1^2)、f=N(μ2,σ22)f = N(\mu_2, \sigma_2^2) のKL情報量を求めます。log⁡(g/f)=log⁡(σ2/σ1)−(z−μ1)2/(2σ12)+(z−μ2)2/(2σ22)\log(g/f) = \log(\sigma_2/\sigma_1) - (z - \mu_1)^2/(2\sigma_1^2) + (z - \mu_2)^2/(2\sigma_2^2) の期待値を、gg の下でとります。Eg[(Z−μ2)2]=σ12+(μ1−μ2)2E_g[(Z - \mu_2)^2] = \sigma_1^2 + (\mu_1 - \mu_2)^2 なので、次のようになります。

D(g ∥ f)=log⁡σ2σ1+σ12+(μ1−μ2)22σ22−12D(g \,\|\, f) = \log\frac{\sigma_2}{\sigma_1} + \frac{\sigma_1^2 + (\mu_1 - \mu_2)^2}{2\sigma_2^2} - \frac{1}{2}

例えば、g=N(0,1)g = N(0, 1)、f=N(1,4)f = N(1, 4) なら、log⁡2+2/8−1/2≈0.443\log 2 + 2/8 - 1/2 \approx 0.443 です。逆向きは log⁡(1/2)+5/2−1/2≈1.307\log(1/2) + 5/2 - 1/2 \approx 1.307 で、3倍近く違います。「細い分布から見た太い分布」より、「太い分布から見た細い分布」のほうが、ズレは大きいのです。太い分布のデータは、細いモデルがほとんど想定しない値をたくさん出すからです。

せっかくなので、別の例も見てみましょう。平均 λ1\lambda_1 と λ2\lambda_2 のポアソン分布では、同じ計算で D=λ1log⁡(λ1/λ2)+λ2−λ1D = \lambda_1\log(\lambda_1/\lambda_2) + \lambda_2 - \lambda_1 です。λ1=2\lambda_1 = 2、λ2=3\lambda_2 = 3 なら約0.189です。

ここで、KL情報量を2つに分けます。

D(g ∥ fθ)=Eg[log⁡g(Z)]−Eg[log⁡f(Z;θ)]D(g \,\|\, f_\theta) = E_g[\log g(Z)] - E_g[\log f(Z;\theta)]

第1項はモデルによらないので、モデルの比較では第2項だけを見れば足ります。第2項の η(θ)=Eg[log⁡f(Z;θ)]\eta(\theta) = E_g[\log f(Z;\theta)] を 平均対数尤度 と呼びます。これが「本番の成績」にあたります。ZZ は、手元のデータとは別に、新しく出てくる本番の問題です。

η(θ)\eta(\theta) を最大にする θ∗\theta^* を、擬似真値 と呼びます。モデルが間違っていても、モデルの中で真の分布に最も近い点です。

例えば、真の分布が平均1の指数分布なのに、正規モデル N(μ,σ2)N(\mu, \sigma^2) を使うとしましょう。η=−12log⁡(2πσ2)−Eg[(Z−μ)2]/(2σ2)\eta = -\frac{1}{2}\log(2\pi\sigma^2) - E_g[(Z - \mu)^2]/(2\sigma^2) を最大にすると、μ∗=Eg[Z]=1\mu^* = E_g[Z] = 1、σ∗2=Vg[Z]=1\sigma^{*2} = V_g[Z] = 1 です。つまり、平均と分散をそろえた正規分布です。このときのKL情報量は、Eg[log⁡g]=−1E_g[\log g] = -1 より、12log⁡(2π)−12≈0.419\frac{1}{2}\log(2\pi) - \frac{1}{2} \approx 0.419 となります。

以上が、カルバック・ライブラー情報量についての解説です。

情報量規準AIC:罰則を導く

次に、この記事の中心である情報量規準AICの導出を説明します。

最尤推定量 θ^\hat{\theta} を使ったモデルの本番の成績は、n η(θ^)n\,\eta(\hat{\theta}) です(nn 問分に換算しました)。でも、gg が未知なので η\eta は計算できません。そこで、手元の成績である最大対数尤度 ℓ(θ^)\ell(\hat{\theta}) で代用します。問題は、この2つの差です。

b=E[ℓ(θ^)−n η(θ^)]b = E\big[\ell(\hat{\theta}) - n\,\eta(\hat{\theta})\big]

過去問で鍛えた受験生の、過去問の点と本番の点の差の平均です。これを計算します。

準備として、記号を2つ置きます。擬似真値 θ∗\theta^* で評価した、次の2つの p×pp \times p 行列です(pp は母数の個数)。

J=−Eg[∂2log⁡f(Z;θ∗)∂θ ∂θ⊤],K=Eg[∂log⁡f(Z;θ∗)∂θ∂log⁡f(Z;θ∗)∂θ⊤]J = -E_g\left[\frac{\partial^2 \log f(Z;\theta^*)}{\partial\theta\,\partial\theta^\top}\right], \qquad K = E_g\left[\frac{\partial\log f(Z;\theta^*)}{\partial\theta}\frac{\partial\log f(Z;\theta^*)}{\partial\theta^\top}\right]

U(θ^)=0U(\hat{\theta}) = 0 を θ∗\theta^* のまわりで展開すると、θ^−θ∗≈J−1U(θ∗)/n\hat{\theta} - \theta^* \approx J^{-1}U(\theta^*)/n です。U(θ∗)U(\theta^*) は独立な和で、平均は0、分散行列は nKnK なので、次の近似が成り立ちます。この種の展開は、次回の記事で詳しく証明します。

n E[(θ^−θ∗)(θ^−θ∗)⊤]≈J−1KJ−1n\,E\big[(\hat{\theta} - \theta^*)(\hat{\theta} - \theta^*)^\top\big] \approx J^{-1} K J^{-1}

いよいよ導出です。差を3つの部分に分けます。

ℓ(θ^)−nη(θ^)=ℓ(θ^)−ℓ(θ∗)⏟A+ℓ(θ∗)−nη(θ∗)⏟B+nη(θ∗)−nη(θ^)⏟C\ell(\hat{\theta}) - n\eta(\hat{\theta}) = \underbrace{\ell(\hat{\theta}) - \ell(\theta^*)}_{A} + \underbrace{\ell(\theta^*) - n\eta(\theta^*)}_{B} + \underbrace{n\eta(\theta^*) - n\eta(\hat{\theta})}_{C}

1つずつ見ていきましょう。

【3つの部分の期待値】

  1. BB: θ∗\theta^* は定数なので、E[ℓ(θ∗)]=nη(θ∗)E[\ell(\theta^*)] = n\eta(\theta^*) で、E[B]=0E[B] = 0 である
  2. AA: ℓ\ell は θ^\hat{\theta} が山頂なので、ℓ(θ∗)≈ℓ(θ^)−n2(θ^−θ∗)⊤J(θ^−θ∗)\ell(\theta^*) \approx \ell(\hat{\theta}) - \frac{n}{2}(\hat{\theta} - \theta^*)^\top J (\hat{\theta} - \theta^*) である。よって E[A]≈12tr(KJ−1)E[A] \approx \frac{1}{2}\text{tr}(KJ^{-1}) である
  3. CC: η\eta は θ∗\theta^* が山頂なので、η(θ^)≈η(θ∗)−12(θ^−θ∗)⊤J(θ^−θ∗)\eta(\hat{\theta}) \approx \eta(\theta^*) - \frac{1}{2}(\hat{\theta} - \theta^*)^\top J (\hat{\theta} - \theta^*) である。よって E[C]≈12tr(KJ−1)E[C] \approx \frac{1}{2}\text{tr}(KJ^{-1}) である

期待値の計算には、E[v⊤Jv]=tr(J E[vv⊤])E[\boldsymbol{v}^\top J\boldsymbol{v}] = \text{tr}(J\,E[\boldsymbol{v}\boldsymbol{v}^\top]) を使いました。合計すると、b≈tr(KJ−1)b \approx \text{tr}(KJ^{-1}) です。

ここで、モデルが正しいとします。すると θ∗\theta^* は真の値で、K=J=I(θ)K = J = I(\theta)(フィッシャー情報行列)です。これは有効スコア関数の等式の多変量版です。したがって、tr(KJ−1)=tr(Ip)=p\text{tr}(KJ^{-1}) = \text{tr}(I_p) = p となります。

n η(θ^)≈ℓ(θ^)−p,AIC=−2ℓ(θ^)+2pn\,\eta(\hat{\theta}) \approx \ell(\hat{\theta}) - p, \qquad \text{AIC} = -2\ell(\hat{\theta}) + 2p

これが罰則の正体です。甘さ pp の内訳は、AA の p/2p/2 と CC の p/2p/2 です。AA は「過去問に合わせすぎて上がった点」で、CC は「合わせすぎたせいで本番で下がった点」です。2つの山がよく似た形をしているので、同じ量になるのです。最後に −2-2 を掛けるのは、尤度比検定の統計量と同じ尺度にそろえるためです。

例題:補正が厳密に1になることを示す

X1,…,XnX_1, \dots, X_n は独立に N(μ,1)N(\mu, 1) に従い、モデルは N(θ,1)N(\theta, 1)(母数は θ\theta の1個)である。E[ℓ(θ^)−nη(θ^)]=1E[\ell(\hat{\theta}) - n\eta(\hat{\theta})] = 1 を示せ。

【答案】

  1. θ^=Xˉ\hat{\theta} = \bar{X} である。ℓ(Xˉ)=−n2log⁡(2π)−12∑i(Xi−Xˉ)2\ell(\bar{X}) = -\frac{n}{2}\log(2\pi) - \frac{1}{2}\sum_i (X_i - \bar{X})^2 で、E[∑i(Xi−Xˉ)2]=n−1E[\sum_i (X_i - \bar{X})^2] = n - 1 より、次の式を得る
E[ℓ(Xˉ)]=−n2log⁡(2π)−n−12E[\ell(\bar{X})] = -\frac{n}{2}\log(2\pi) - \frac{n-1}{2}
  1. Z∼N(μ,1)Z \sim N(\mu, 1) として、η(θ)=−12log⁡(2π)−12{1+(θ−μ)2}\eta(\theta) = -\frac{1}{2}\log(2\pi) - \frac{1}{2}\{1 + (\theta - \mu)^2\} である
  2. E[(Xˉ−μ)2]=1/nE[(\bar{X} - \mu)^2] = 1/n より、E[n η(Xˉ)]=−n2log⁡(2π)−n2−12E[n\,\eta(\bar{X})] = -\frac{n}{2}\log(2\pi) - \frac{n}{2} - \frac{1}{2} である
  3. 差をとると、−n−12+n2+12=1-\frac{n-1}{2} + \frac{n}{2} + \frac{1}{2} = 1 である

答えは、近似なしでちょうど1です。この例では、E[A]=E[C]=1/2E[A] = E[C] = 1/2 も厳密に成り立ちます。

では、モデルが間違っているときはどうでしょうか? 正規モデルで、真の分布の尖度が κ=μ4/σ4\kappa = \mu_4/\sigma^4 だとします。JJ と KK を計算すると、tr(KJ−1)=1+(κ−1)/2\text{tr}(KJ^{-1}) = 1 + (\kappa - 1)/2 です。真が正規分布なら κ=3\kappa = 3 で値は2で、母数の個数に一致します。ところが、真が指数分布なら κ=9\kappa = 9 で値は5です。AICの罰則は、このときの甘さを小さく見積もっています。

補足(ここは読み飛ばしてOK): pp の代わりに tr(K^J^−1)\text{tr}(\hat{K}\hat{J}^{-1}) を使う規準を、TIC(竹内の情報量規準)と呼びます。また、AICの罰則を plog⁡np\log n に重くしたBICは、一致性をもちます。

AICの性質で、もう1つ大切なことがあります。真のモデルに不要な母数を1つ足したモデルと比べると、−2ℓ-2\ell の差は漸近的に χ2(1)\chi^2(1) に従います。AICが大きいほうのモデルを選ぶのは、この値が2を超えるときです。P(χ2(1)>2)≈0.157P(\chi^2(1) > 2) \approx 0.157 は nn が増えても0になりません。つまり、AICは真のモデルを選ぶ一致性をもちません。AICの目的は真のモデルの特定ではなく、本番の成績の見積もりだからです。

試験のポイント:AICの導出は「甘さを3つに分ける」「E[A]E[A] と E[C]E[C] が p/2p/2 ずつ」「モデルが正しいと K=JK = J」の3段階です。この流れを答案に書けるようにしましょう。

以上が、情報量規準AICの導出についての説明です。

クロスバリデーション:1個抜きの公式を導く

最後に、もう1つの方法であるクロスバリデーションの理論を見ていきましょう。

AICは数式で甘さを補正しました。一方、クロスバリデーション は、データの一部を隠して、本番を実際に再現します。過去問の一部を解かずにとっておいて、模試として使うイメージです。

【kk 分割クロスバリデーションの手順】

  1. データを kk 個のグループに分ける
  2. 1グループを検証用に残し、残りでモデルを推定する
  3. 検証用のデータで予測誤差を計算する
  4. 検証用のグループを順に変えて、誤差を平均する

k=nk = n のときが、1個ずつ抜く 1個抜きクロスバリデーション(LOOCV)です。損失を2乗誤差とすると、次の量になります。

CV=1n∑i=1n(yi−y^i(−i))2\text{CV} = \frac{1}{n}\sum_{i=1}^{n}\big(y_i - \hat{y}_i^{(-i)}\big)^2

y^i(−i)\hat{y}_i^{(-i)} は、ii 番目を除いて推定したモデルの、ii 番目の予測値です。素直に計算すると nn 回の推定が必要ですが、線形模型なら1回で済みます。ei=yi−y^ie_i = y_i - \hat{y}_i を通常の残差、hiih_{ii} をハット行列の対角成分とすると、次の式が成り立ちます。

yi−y^i(−i)=ei1−hiiy_i - \hat{y}_i^{(-i)} = \frac{e_i}{1 - h_{ii}}

証明には、データを1つ書き換える技を使います。

【証明】

  1. yiy_i だけを y^i(−i)\hat{y}_i^{(-i)} に置き換えたデータ y~\tilde{\boldsymbol{y}} を考える
  2. y~\tilde{\boldsymbol{y}} の最小二乗推定量は、ii 番目を除いた推定量 β^(−i)\hat{\boldsymbol{\beta}}^{(-i)} に一致する。β^(−i)\hat{\boldsymbol{\beta}}^{(-i)} は ii 番目以外の2乗和を最小にし、ii 番目の残差も0にするからである
  3. よって y^i(−i)\hat{y}_i^{(-i)} は y~\tilde{\boldsymbol{y}} の予測値で、y^i(−i)=∑jhijy~j=y^i−hiiyi+hiiy^i(−i)\hat{y}_i^{(-i)} = \sum_j h_{ij}\tilde{y}_j = \hat{y}_i - h_{ii}y_i + h_{ii}\hat{y}_i^{(-i)} である
  4. 整理すると (1−hii)y^i(−i)=y^i−hiiyi(1 - h_{ii})\hat{y}_i^{(-i)} = \hat{y}_i - h_{ii}y_i である。(1−hii)yi(1 - h_{ii})y_i から引けば、(1−hii)(yi−y^i(−i))=yi−y^i(1 - h_{ii})(y_i - \hat{y}_i^{(-i)}) = y_i - \hat{y}_i を得る

hiih_{ii} は てこ比 とも呼ばれ、その点が自分の予測値をどれだけ引っ張るかを表します。声の大きい人ほど、抜けたときの変化が大きいのです。

例題:公式で1個抜きCVを計算する

データ (x,y)(x, y) が (1, 1), (2, 3), (3, 2), (4, 5), (5, 4) である。単回帰モデルと定数のみのモデルの1個抜きCVを、公式で計算して比べよ。

【答案】

  1. 単回帰の結果は y^=0.6+0.8x\hat{y} = 0.6 + 0.8x で、残差は −0.4,0.8,−1,1.2,−0.6-0.4, 0.8, -1, 1.2, -0.6 である
  2. hii=1/n+(xi−xˉ)2/Sxxh_{ii} = 1/n + (x_i - \bar{x})^2/S_{xx} より、てこ比は順に 0.6,0.3,0.2,0.3,0.60.6, 0.3, 0.2, 0.3, 0.6 である
  3. ei/(1−hii)e_i/(1 - h_{ii}) は −1,1.143,−1.25,1.714,−1.5-1, 1.143, -1.25, 1.714, -1.5 で、2乗の平均は約1.81である
  4. 定数のみのモデルは y^=3\hat{y} = 3、hii=0.2h_{ii} = 0.2 である。ei/(1−hii)e_i/(1 - h_{ii}) は −2.5,0,−1.25,2.5,1.25-2.5, 0, -1.25, 2.5, 1.25 で、2乗の平均は3.125である

答えは単回帰が約1.81、定数のみが3.125で、単回帰を選びます。なお、単回帰の手元の成績 ∑iei2/n\sum_i e_i^2/n は0.72でした。手元の成績より、模試の成績のほうが2.5倍ほど悪いことが分かります。これが甘さの大きさです。

最後に、AICとの関係を見ましょう。hiih_{ii} をその平均 p/np/n で置き換えると、CV≈(RSS/n)/(1−p/n)2\text{CV} \approx (\text{RSS}/n)/(1 - p/n)^2 です。対数をとって nn 倍すると、次のようになります。

nlog⁡CV≈nlog⁡RSSn−2nlog⁡(1−pn)≈nlog⁡RSSn+2pn\log\text{CV} \approx n\log\frac{\text{RSS}}{n} - 2n\log\left(1 - \frac{p}{n}\right) \approx n\log\frac{\text{RSS}}{n} + 2p

右辺は、正規線形回帰のAICから定数を除いたものです。答え合わせのように、計算で補正したAICと、実際に隠してみたCVが一致しました。一般に、対数尤度を損失とする1個抜きCVとAICは、漸近的に同等です。

方法見積もり方強み弱み
AIC甘さを pp で補正計算が1回で済む尤度が必要。モデルのズレに弱い
クロスバリデーションデータを隠して再現尤度が不要で、損失を選べる計算が多い。時系列などでは分け方に注意

以上が、クロスバリデーションについての解説です。

細かい展開を忘れてしまっても問題ありません。ここでは、「手元の成績は平均して母数の個数だけ甘い」という点だけ押さえておきましょう。

今回のまとめ

それでは、今回の内容を振り返ります。

今回は、KL情報量の性質を証明し、AICの罰則を導き、クロスバリデーションの公式を証明しました。AICとクロスバリデーションは、どちらも「過去問の点から本番の点を見積もる」方法です。式で補正するか、問題を隠すかの違いで、行き着く先は同じでした。

今回のポイントを振り返りましょう。

振り返り

  • KL情報量の非負性は log⁡x≤x−1\log x \le x - 1 で証明できる。モデルの比較には平均対数尤度だけを見ればよい。
  • 手元の成績の甘さは tr(KJ−1)\text{tr}(KJ^{-1}) で、モデルが正しければ母数の個数 pp になる。
  • 甘さの内訳は、過去問で上がる p/2p/2 と本番で下がる p/2p/2 である。
  • AICは一致性をもたないが、本番の成績の見積もりとしてよい性質をもつ。
  • 線形模型の1個抜きCVは ei/(1−hii)e_i/(1 - h_{ii}) で計算でき、AICと漸近的に同等である。

キーワードのおさらい

  • カルバック・ライブラー情報量: Eg[log⁡(g/f)]E_g[\log(g/f)] で定める、真の分布から見たモデルのズレ
  • 平均対数尤度: 新しいデータでの対数尤度の期待値 Eg[log⁡f(Z;θ)]E_g[\log f(Z;\theta)]
  • 擬似真値: 平均対数尤度を最大にする母数の値
  • 情報量規準AIC: −2ℓ(θ^)+2p-2\ell(\hat{\theta}) + 2p で定める、本番の成績の見積もり
  • クロスバリデーション: データの一部を隠して予測誤差を測る方法
  • てこ比: ハット行列の対角成分で、その点が自分の予測値に与える影響

さて、AICの導出では、最尤推定量が真の値の近くで正規分布に近づくことを先取りして使いました。この性質は、どう証明すればよいのでしょうか?

次回は、クラーメル・ラオの不等式と最尤推定量の漸近正規性を証明し、デルタ法を導きます。

それでは、次の記事へ進んでいきましょう。

練習問題5問開いて挑戦する
問1. 平均4のポアソン分布を gg、平均2のポアソン分布を ff とする。D(g ∥ f)D(g \,\|\, f) の値に最も近いものはどれか。
問2. 最大対数尤度が −120.5-120.5、母数の個数が4のモデルのAICとして正しいものはどれか。
問3. 真の分布の尖度が κ=5\kappa = 5 で、正規モデル(母数2個)を当てはめる。手元の成績の甘さ tr(KJ−1)=1+(κ−1)/2\text{tr}(KJ^{-1}) = 1 + (\kappa - 1)/2 の値として正しいものはどれか。
問4. 線形模型のある点の残差が ei=1.2e_i = 1.2、てこ比が hii=0.4h_{ii} = 0.4 である。その点を除いて推定したモデルでの予測誤差として正しいものはどれか。
問5. n=50n = 50 の正規線形回帰で、係数2個のモデルAの残差平方和が200、係数3個のモデルBの残差平方和が190である。nlog⁡(RSS/n)+2(係数の個数)n\log(\text{RSS}/n) + 2(\text{係数の個数}) で比べたとき、正しいものはどれか。