Zenn
🟦

情報量規準の数理を理解したい

2024/04/08に公開

はじめに

この記事は情報量規準の定義にまつわる数理を理解するために試行錯誤した記録である.

もともと統計モデル・機械学習モデルやその数理に興味があり,統計検定準1級の勉強をする中で出てきたAIC赤池情報量規準)や,より一般の情報量規準について深掘りしたくなった.以下の書籍を手に取ったところ「まさにこれだ」と感銘を受けたので,重要と思われる主に3.4節の内容を中心に行間を自分なりに噛み砕き,一部の数値例についてPythonによる簡単なシミュレーションを行った上で,記事の形にアウトプットしておく.

https://www.asakura.co.jp/detail.php?book_code=12782

(間違いがあるかもしれませんがその際はご指摘いただけますと幸いです.)

要約

  • 最尤法によりモデルを推定するケースを考える.
  • 推定された最尤モデルが真のモデルを予測の意味でよく近似していることを,Kullback–Leibler情報量ではかる.
  • Kullback–Leibler情報量の第2項に出てくる平均対数尤度を,モデル推定に使った最大対数尤度を使って推定する.
  • しかし,最大対数尤度だけではモデル選択基準とするには不十分であり,補正項を考慮する必要がある.
  • 情報量規準はそうした補正項を精密に評価した上で定義される.(次回の記事を書く機会があれば触れたい)

モデルの最尤推定

真の確率分布をGG,その確率密度関数をggとする.この現象の構造を捉えるために,データ分析者はパラメータθΘ(Rp)\bm{\theta} \in \Theta (\subset \mathbb{R}^p)によりf(θ)f(\cdot | \bm{\theta})という確率密度関数で表されるモデルを想定する.そしてggを何らかの意味でよく再現するようなモデルf(θ)f(\cdot | \bm{\theta})を探したい状況を考える.その際によく用いられるのが最尤法である.

GGからデータxn={x1,,xn}\bm{x}_n = \{ x_1, \cdots, x_n \}が得られているとし,各θ\bm{\theta}に対して定まる対数尤度

(θ)=logi=1nf(xiθ)=i=1nlogf(xiθ) \ell(\bm{\theta}) = \log \prod_{i=1}^n f(x_i | \bm{\theta}) = \sum_{i=1}^n \log f(x_i | \bm{\theta})

を最大化するパラメータ

θ^=θ^(xn)=arg maxθΘ(θ) \hat{\bm{\theta}} = \hat{\bm{\theta}}(\bm{x}_n) = \argmax_{\theta \in \Theta}{\ell(\bm{\theta})}

を持つモデルf(θ^)f(\cdot | \hat{\bm{\theta}})を採用する.このようなパラメータ推定を最尤推定という.

モデルの"良さ"とKullback–Leibler情報量

さて,こうしてデータxn\bm{x}_nからできたモデルf(θ^)f(\cdot | \hat{\bm{\theta}})はどのくらい"良い"のだろうか?GGから新しい未知のデータZ=zZ = zが得られた時,f(zθ^)f(z | \hat{\bm{\theta}})g(z)g(z)をうまく予測できるようなモデルが望ましいだろう.それをはかる際の指標として,二つの確率分布間の"距離"[1]をはかるKullback–Leibler情報量を考える.

さて,g(z)g(z)f(zθ^)f(z | \hat{\bm{\theta}})Kullback–Leibler情報量

I(g(z);f(zθ^)):=EG[logg(Z)f(Zθ^)]=EG[logg(Z)]EG[logf(Zθ^)] \begin{align*} I(g(z); f(z | \hat{\bm{\theta}})) &:= E_G \left[ \log \frac{g(Z)}{f(Z | \hat{\bm{\theta}})} \right] \\ &= E_G[\log g(Z)] - E_G[\log f(Z | \hat{\bm{\theta}})] \end{align*}

と定義される.値が小さいほど二つの分布が近いことを意味する.第一項は未知の分布ggにのみ依存しており分析者が制御できないため,注目するのは第二項EG[logf(Zθ^)]E_G[\log f(Z | \hat{\bm{\theta}})]である.これを平均対数尤度と呼ぶ.平均対数尤度が大きいほど,構築したモデルが真の分布をよく近似していると言えるだろう.

しかし,平均対数尤度も依然として真の確率分布GGに依存している.そこで,GGから得られたデータxn={x1,,xn}\bm{x}_n = \{ x_1, \cdots, x_n \}から平均対数尤度を推定することを考えよう.それには経験分布を考える.nn個のデータの各点で等確率1/n1 / nをもつ確率関数g^(xi)=1/n(i=1,,n)\hat{g}(x_i) = 1 / n (i = 1, \cdots, n)を用いて未知の確率分布GGを経験分布G^\hat{G}で置き換えた推定量を考えると

EG^[logf(Zθ^)]=logf(zθ^)dG^(z)=zxnlogf(zθ^)g^(z)=1ni=1nlogf(xiθ^)=1n(θ^) \begin{align*} E_{\hat{G}}[\log f(Z | \hat{\bm{\theta}})] &= \int \log f(z | \hat{\bm{\theta}}) d\hat{G}(z) \\ &= \sum_{z \in \bm{x}_n} \log f(z | \hat{\bm{\theta}}) \hat{g}(z) \\ &= \frac{1}{n} \sum_{i=1}^n \log{f(x_i | \hat{\bm{\theta}})} = \frac{1}{n} \ell(\hat{\bm{\theta}}) \end{align*}

となり,両辺をnn倍することで

(θ^)=nEG^[logf(Zθ^)] \ell(\hat{\bm{\theta}}) = n E_{\hat{G}}[\log f(Z | \hat{\bm{\theta}})]

を得る.

以上により,「平均対数尤度」の推定量として「最大対数尤度」を採用してよさそうである.

最大対数尤度の補正

複数,例えば2つのモデルf1(θ1),f2(θ2)f_1(\cdot | \bm{\theta}_1), f_2(\cdot | \bm{\theta}_2)と各々のパラメータの最尤推定値θ^1,θ^2\hat{\bm{\theta}}_1, \hat{\bm{\theta}}_2が得られている状況を考える.上記の議論を踏まえると,最大対数尤度1(θ^1)\ell_1(\hat{\bm{\theta}}_1)2(θ^2)\ell_2(\hat{\bm{\theta}}_2)の大小比較を行うことでどちらのモデルが良いかを判定できるのではないか?

しかし,そのまま単純に最大対数尤度を考えるだけではモデル比較に使えない.以下の等式

nEG(z)[logf(Zθ^(xn))]平均対数尤度=(θ^(xn))最大対数尤度{(θ^(xn))nEG(z)[logf(Zθ^(xn))]}補正項 n \underbrace{ E_{G(z)}[\log f(Z | \hat{\bm{\theta}}(\bm{x}_n))] }_{平均対数尤度} = \underbrace{ \ell(\hat{\bm{\theta}}(\bm{x}_n)) }_{最大対数尤度} - \underbrace{ \left\{ \ell(\hat{\bm{\theta}}(\bm{x}_n)) - n E_{G(z)}[\log f(Z | \hat{\bm{\theta}}(\bm{x}_n))] \right\}}_{補正項}

が成り立つことから,最大対数尤度に対して補正項の考慮を行うことが,モデル選択基準としての情報量規準の定義につながっていく.

なお,一般には真の分布ggが未知であるため,補正項をデータから計算することはできない.そこで,理論的には補正項(θ^(Xn))nEG(z)[logf(Zθ^(Xn))]\displaystyle \ell(\hat{\bm{\theta}}(\bm{X}_n)) - n E_{G(z)}[\log f(Z | \hat{\bm{\theta}}(\bm{X}_n))]に対し,期待値をとった量

b(G)=EG(xn)[(θ^(Xn))nEG(z)[logf(Zθ^(Xn))]] b(G) = E_{G(\bm{x}_n)} \left[ \ell(\hat{\bm{\theta}}(\bm{X}_n)) - n E_{G(z)}[\log f(Z | \hat{\bm{\theta}}(\bm{X}_n))] \right]

を介する.

【具体例で考察】正規分布の場合

参考書籍にg,fg, fどちらも正規分布の場合が載っているので,数式の行間を埋め,シミュレーションを行ってみる.

分析者が現象を正規分布モデルで説明しようとした状況を考えよう.モデルの確率密度関数はf(zμ,σ2)=12πσ2e(zμ)22σ2\displaystyle f(z | \mu, \sigma^2) = \frac{1}{\sqrt{2 \pi \sigma^2}} e^{- \frac{(z-\mu)^2}{2 \sigma^2}}である[2].データxn={x1,,xn}\bm{x}_n = \{ x_1, \cdots, x_n \}に基づくモデルの対数尤度は

(μ,σ2)=logi=1nf(xiμ,σ2)=n2log(2πσ2)12σ2i=1n(xiμ)2 \ell(\mu, \sigma^2) = \log \prod_{i=1}^n f(x_i | \mu, \sigma^2) = -\frac{n}{2} \log(2 \pi \sigma^2) - \frac{1}{2 \sigma^2} \sum_{i=1}^n (x_i-\mu)^2

である.最尤法によって最尤推定値

μ^=1ni=1nxi,σ^2=1ni=1n(xiμ^)2 \hat{\mu} = \frac{1}{n} \sum_{i=1}^n x_i, \quad \hat{\sigma}^2 = \frac{1}{n} \sum_{i=1}^n ( x_i - \hat{\mu} )^2

が求まり,最大対数尤度

(μ^,σ^2)=n2log(2πσ^2)12σ^2i=1n(xiμ^)2=n2log(2πσ^2)n2 \begin{align*} \ell(\hat{\mu}, \hat{\sigma}^2) &= -\frac{n}{2} \log(2 \pi \hat{\sigma}^2) - \frac{1}{2 \hat{\sigma}^2} \sum_{i=1}^n ( x_i - \hat{\mu} )^2 \\ &= -\frac{n}{2} \log(2 \pi \hat{\sigma}^2) - \frac{n}{2} \end{align*}

と最尤モデルf(zμ^,σ^2)=12πσ^2e(zμ^)22σ^2\displaystyle f(z | \hat{\mu}, \hat{\sigma}^2) = \frac{1}{\sqrt{2 \pi \hat{\sigma}^2}} e^{- \frac{(z-\hat{\mu})^2}{2 \hat{\sigma}^2}}が求まる.これにより

logf(zμ^,σ^2)=12log(2πσ^2)12σ^2(zμ^)2 \log f(z | \hat{\mu}, \hat{\sigma}^2) = -\frac{1}{2} \log(2 \pi \hat{\sigma}^2) - \frac{1}{2 \hat{\sigma}^2} (z-\hat{\mu})^2

である.

ここまでは真の分布ggが何であっても成り立つ議論であった.以降,真の分布ggが正規分布N(μg,σg2)N(\mu_g, \sigma_g^2)であるとしよう[3].平均対数尤度は,xn\bm{x}_nを固定したままzN(μg,σg2)z \sim N(\mu_g, \sigma_g^2)に関して期待値をとることで

EG(z)[logf(zμ^,σ^2)]=12log(2πσ^2)12σ^2EG(z)[(zμ^)2]=12log(2πσ^2)12σ^2{σg2+(μgμ^)2} \begin{align*} E_{G(z)}[\log f(z | \hat{\mu}, \hat{\sigma}^2)] &= -\frac{1}{2} \log(2 \pi \hat{\sigma}^2) - \frac{1}{2 \hat{\sigma}^2} E_{G(z)}[(z-\hat{\mu})^2] \\ &= -\frac{1}{2} \log(2 \pi \hat{\sigma}^2) - \frac{1}{2 \hat{\sigma}^2} \{ \sigma_g^2 + (\mu_g-\hat{\mu})^2 \} \end{align*}

と計算できる.

補正項,つまり最大対数尤度と平均対数尤度(のnn倍)の差は

(μ^,σ^2)nEG(z)[logf(zμ^,σ^2)]=n2σ^2{σg2+(μgμ^)2}n2 \ell(\hat{\mu}, \hat{\sigma}^2) - n E_{G(z)}[\log f(z | \hat{\mu}, \hat{\sigma}^2)] = \frac{n}{2 \hat{\sigma}^2} \{ \sigma_g^2 + (\mu_g-\hat{\mu})^2 \} - \frac{n}{2}

となる.期待値をとることで

b(G)=EG(xn)[n2σ^2{σg2+(μgμ^)2}n2]=n2EG(xn)[1σ^2/σg2+1σ^2(μ^μg)2]n2=n2EG(xn)[nχn12+1σg2nχn12(μ^μg)2]n2=n2{E[nχn12]+1σg2E[nχn12]E[(μ^μg)2]}n2=n2{nn3+1σg2nn3σg2n}n2=2nn3 \begin{align*} b(G) &= E_{G(\bm{x}_n)}\left[ \frac{n}{2 \hat{\sigma}^2} \{ \sigma_g^2 + (\mu_g-\hat{\mu})^2 \} - \frac{n}{2} \right] \\ &= \frac{n}{2} E_{G(\bm{x}_n)}\left[ \frac{1}{\hat{\sigma}^2 / \sigma_g^2} + \frac{1}{\hat{\sigma}^2} \cdot (\hat{\mu} - \mu_g)^2 \right] - \frac{n}{2} \\ &= \frac{n}{2} E_{G(\bm{x}_n)}\left[ \frac{n}{\chi^2_{n-1}} + \frac{1}{\sigma_g^2} \frac{n}{\chi^2_{n-1}} \cdot (\hat{\mu} - \mu_g)^2 \right] - \frac{n}{2} \\ &= \frac{n}{2} \left\{ E\left[ \frac{n}{\chi^2_{n-1}} \right] + \frac{1}{\sigma_g^2} E\left[ \frac{n}{\chi^2_{n-1}} \right] \cdot E\left[ (\hat{\mu} - \mu_g)^2 \right] \right\} - \frac{n}{2} \\ &= \frac{n}{2} \left\{ \frac{n}{n-3} + \frac{1}{\sigma_g^2} \cdot \frac{n}{n-3} \cdot \frac{\sigma_g^2}{n} \right\} - \frac{n}{2} \\ &= \frac{2 n}{n - 3} \end{align*}

となる[4].上記の式変形では,

σ^2σg2=1ni=1n(xiμ^σg)2=1nχn12 \frac{\hat{\sigma}^2}{\sigma_g^2} = \frac{1}{n} \sum_{i=1}^n \left( \frac{x_i - \hat{\mu}}{\sigma_g} \right)^2 = \frac{1}{n} \chi^2_{n-1}

であり,E[1/χr2]=1/(r2)E[1 / \chi^2_r] = 1 / (r-2)が成り立つこと(χr2\chi^2_rは自由度rrのカイ二乗統計量)と,標本平均と標本分散が独立であることを証明無しに認めて用いた.

数式の議論まとめ

以上をまとめる.正規分布N(μg,σg2)N(\mu_g, \sigma_g^2)に従うデータxn={x1,,xn}\bm{x}_n = \{ x_1, \cdots, x_n \}があるとする.このデータを正規分布モデルで表現することを想定する.最尤推定により推定された最尤モデルf(μ^,σ^2)f(\cdot | \hat{\mu}, \hat{\sigma}^2)の良さを「最大対数尤度 - 補正項」で評価する.ここで

補正項=n2σ^2{σg2+(μgμ^)2}n2 補正項 = \frac{n}{2 \hat{\sigma}^2} \{ \sigma_g^2 + (\mu_g-\hat{\mu})^2 \} - \frac{n}{2}

であり,その期待値は

b(G)=2nn3 b(G) = \frac{2 n}{n - 3}

である.ちなみにnn \to \inftyにおいてb(G)2b(G) \approx 2であり,これは正規分布モデルのパラメータ数に等しい.

Pythonでのシミュレーション

補正項とその期待値をシミュレーションしてみよう.環境は以下の通り.

Python 3.8.13
numpy 1.23.5
import numpy as np

# 真の分布gの母数
mu = 3
sigma = 0.1

# 補正項を計算する関数
def correction_calc(mu=0, sigma=0.1, n=100, seed=0):
    
    # ランダムシードの固定(0または未入力だったらスキップ)
    if seed != 0:
        np.random.seed(seed)
    
    # 真の分布から正規乱数の生成
    x = np.random.normal(loc=mu, scale=sigma, size=n)
    
    # 標本平均・標本分散(ともに最尤推定値)の計算
    x_mean = np.mean(x)
    x_var = np.var(x)
    
    # 補正項の計算
    correction = (n / 2) * ((sigma**2 + (mu - x_mean)**2) / x_var - 1)
    
    return correction

この準備のもと,サンプルサイズnnを変えながら補正項を5つずつ計算していく.(ランダムシードを固定している点はご容赦を.)

n=10n=10

# サンプルサイズを指定
n = 10

# 補正項を5つ計算
for i in range(1, 6):
     print(correction_calc(mu, sigma, n, i))
-1.4412341801890833
1.289367281291841
1.7921243249806118
4.009671376323563
0.5121271185925691

n=1000n=1000

# サンプルサイズを指定
n = 1000

# 補正項を5つ計算
for i in range(1, 6):
     print(correction_calc(mu, sigma, n, i))
20.333828816127774
-3.2085759986053564
-8.150289753807993
53.105746883741475
10.831668304333908

この2つの場合だけでも,サンプルサイズnnが大きいほど補正項の分散が大きそうなことがわかる.

ここまでは補正項を5つだけ計算していたが,より多くの補正項を計算し,期待値としてその標本平均を求めてみる.

n=10n=10

# サンプルサイズを指定
n = 10

list = []

# 補正項を100回計算
for i in range(1, 100):
     list.append(correction_calc(mu, sigma, n, i))

print(np.mean(list))
print(np.var(list))
2.4730100700241002
22.93279328630055
# サンプルサイズを指定
n = 10

list = []

# 補正項を10000回計算
for i in range(1, 10000):
     list.append(correction_calc(mu, sigma, n, i))

print(np.mean(list))
print(np.var(list))
2.9215396200060386
27.070986686011192

n=1000n=1000

# サンプルサイズを指定
n = 1000

list = []

# 補正項を100回計算
for i in range(1, 100):
     list.append(correction_calc(mu, sigma, n, i))

print(np.mean(list))
print(np.var(list))
2.436518807833855
492.5831353030407
# サンプルサイズを指定
n = 1000

list = []

# 補正項を10000回計算
for i in range(1, 10000):
     list.append(correction_calc(mu, sigma, n, i))

print(np.mean(list))
print(np.var(list))
1.7674007255551742
501.1739318269931

先ほど見たようにサンプルサイズnnが大きいほど標本分散が大きくなっているが,標本平均は理論値22の周りをうろうろしているように見えなくもない.

なお,理論的には期待値EG(xn)[]E_{G(\bm{x}_n)}[\cdot]をとるが,このシミュレーションでは標本平均を考えているので,もしかしたら見逃している事項があるかもしれない.

今後について

情報量規準の定義に至るには,

b(G)=EG(xn)[(θ^(Xn))nEG(z)[logf(Zθ^(Xn))]]=EG(xn)[(θ^(Xn))(θ0)]+EG(xn)[(θ0)nEG(z)[logf(Zθ0)]]+EG(xn)[nEG(z)[logf(Zθ0)]nEG(z)[logf(Zθ^(Xn))]]=D1+D2+D3 \begin{align*} b(G) &= E_{G(\bm{x}_n)} \left[ \ell(\hat{\bm{\theta}}(\bm{X}_n)) - n E_{G(z)}[\log f(Z | \hat{\bm{\theta}}(\bm{X}_n))] \right] \\ &= E_{G(\bm{x}_n)} \left[ \ell(\hat{\bm{\theta}}(\bm{X}_n)) - \ell(\bm{\theta}_0) \right] \\ &\quad + E_{G(\bm{x}_n)} \left[ \ell(\bm{\theta}_0) - n E_{G(z)}[\log f(Z | \bm{\theta}_0)] \right] \\ &\quad + E_{G(\bm{x}_n)} \left[ n E_{G(z)}[\log f(Z | \bm{\theta}_0)] - n E_{G(z)}[\log f(Z | \hat{\bm{\theta}}(\bm{X}_n))] \right] \\ &= D_1 + D_2 + D_3 \end{align*}

という分解を行い,D2=0D_2 = 0であること,またnn \to \inftyにおいて漸近的にD1=D3=12tr{I(θ0)J(θ0)1}D_1 = D_3 = \frac{1}{2} \mathrm{tr} \{I(\bm{\theta}_0) J(\bm{\theta}_0)^{-1}\}(諸々の定義は次回の記事で行う)が成り立つことを示すのが重要となる.次回の記事を書く機会があれば,この辺りをざっと眺めたい(厳密な証明はできないと思うので).

参考文献

脚注
  1. 数学における距離の公理は満たさない. ↩︎

  2. ここでは正規分布モデルffを作ると決めたことがポイントである.平均と分散という2つの母数があることを認識できていれば,文字としてのμ,σ2\mu, \sigma^2は重要でない.ただし,これらの最尤推定値と,得られる最尤モデルは重要である. ↩︎

  3. ffの分布と明確に分けるためμg,σg2\mu_g, \sigma_g^2という記号を使った. ↩︎

  4. 真の分布がわかっているという,シミュレーション上の仮定があるから計算できる. ↩︎

Discussion

ログインするとコメントできます