はじめに
この記事は情報量規準の定義にまつわる数理を理解するために試行錯誤した記録である．
もともと統計モデル・機械学習モデルやその数理に興味があり，統計検定準1級の勉強をする中で出てきたAIC（赤池情報量規準）や，より一般の情報量規準について深掘りしたくなった．以下の書籍を手に取ったところ「まさにこれだ」と感銘を受けたので，重要と思われる主に3.4節の内容を中心に行間を自分なりに噛み砕き，一部の数値例についてPythonによる簡単なシミュレーションを行った上で，記事の形にアウトプットしておく．
https://www.asakura.co.jp/detail.php?book_code=12782
（間違いがあるかもしれませんがその際はご指摘いただけますと幸いです．）
要約
- 最尤法によりモデルを推定するケースを考える．
- 推定された最尤モデルが真のモデルを予測の意味でよく近似していることを，Kullback–Leibler情報量ではかる．
- Kullback–Leibler情報量の第2項に出てくる平均対数尤度を，モデル推定に使った最大対数尤度を使って推定する．
- しかし，最大対数尤度だけではモデル選択基準とするには不十分であり，補正項を考慮する必要がある．
- 情報量規準はそうした補正項を精密に評価した上で定義される．（次回の記事を書く機会があれば触れたい）
モデルの最尤推定
真の確率分布をG，その確率密度関数をgとする．この現象の構造を捉えるために，データ分析者はパラメータ\bm{\theta} \in \Theta (\subset \mathbb{R}^p)によりf(\cdot | \bm{\theta})という確率密度関数で表されるモデルを想定する．そしてgを何らかの意味でよく再現するようなモデルf(\cdot | \bm{\theta})を探したい状況を考える．その際によく用いられるのが最尤法である．
Gからデータ\bm{x}_n = \{ x_1, \cdots, x_n \}が得られているとし，各\bm{\theta}に対して定まる対数尤度
\ell(\bm{\theta}) = \log \prod_{i=1}^n f(x_i | \bm{\theta}) = \sum_{i=1}^n \log f(x_i | \bm{\theta})
を最大化するパラメータ
\hat{\bm{\theta}} = \hat{\bm{\theta}}(\bm{x}_n) = \argmax_{\theta \in \Theta}{\ell(\bm{\theta})}
を持つモデルf(\cdot | \hat{\bm{\theta}})を採用する．このようなパラメータ推定を最尤推定という．
モデルの"良さ"とKullback–Leibler情報量
さて，こうしてデータ\bm{x}_nからできたモデルf(\cdot | \hat{\bm{\theta}})はどのくらい"良い"のだろうか？Gから新しい未知のデータZ = zが得られた時，f(z | \hat{\bm{\theta}})でg(z)をうまく予測できるようなモデルが望ましいだろう．それをはかる際の指標として，二つの確率分布間の"距離"をはかるKullback–Leibler情報量を考える．
さて，g(z)とf(z | \hat{\bm{\theta}})のKullback–Leibler情報量は
\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*}
と定義される．値が小さいほど二つの分布が近いことを意味する．第一項は未知の分布gにのみ依存しており分析者が制御できないため，注目するのは第二項E_G[\log f(Z | \hat{\bm{\theta}})]である．これを平均対数尤度と呼ぶ．平均対数尤度が大きいほど，構築したモデルが真の分布をよく近似していると言えるだろう．
しかし，平均対数尤度も依然として真の確率分布Gに依存している．そこで，Gから得られたデータ\bm{x}_n = \{ x_1, \cdots, x_n \}から平均対数尤度を推定することを考えよう．それには経験分布を考える．n個のデータの各点で等確率1 / nをもつ確率関数\hat{g}(x_i) = 1 / n (i = 1, \cdots, n)を用いて未知の確率分布Gを経験分布\hat{G}で置き換えた推定量を考えると
\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*}
となり，両辺をn倍することで
\ell(\hat{\bm{\theta}}) = n E_{\hat{G}}[\log f(Z | \hat{\bm{\theta}})]
を得る．
以上により，「平均対数尤度」の推定量として「最大対数尤度」を採用してよさそうである．
最大対数尤度の補正
複数，例えば2つのモデルf_1(\cdot | \bm{\theta}_1), f_2(\cdot | \bm{\theta}_2)と各々のパラメータの最尤推定値\hat{\bm{\theta}}_1, \hat{\bm{\theta}}_2が得られている状況を考える．上記の議論を踏まえると，最大対数尤度\ell_1(\hat{\bm{\theta}}_1)と\ell_2(\hat{\bm{\theta}}_2)の大小比較を行うことでどちらのモデルが良いかを判定できるのではないか？
しかし，そのまま単純に最大対数尤度を考えるだけではモデル比較に使えない．以下の等式
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\}}_{補正項}
が成り立つことから，最大対数尤度に対して補正項の考慮を行うことが，モデル選択基準としての情報量規準の定義につながっていく．
なお，一般には真の分布gが未知であるため，補正項をデータから計算することはできない．そこで，理論的には補正項\displaystyle \ell(\hat{\bm{\theta}}(\bm{X}_n)) - n E_{G(z)}[\log f(Z | \hat{\bm{\theta}}(\bm{X}_n))]に対し，期待値をとった量
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, fどちらも正規分布の場合が載っているので，数式の行間を埋め，シミュレーションを行ってみる．
分析者が現象を正規分布モデルで説明しようとした状況を考えよう．モデルの確率密度関数は\displaystyle f(z | \mu, \sigma^2) = \frac{1}{\sqrt{2 \pi \sigma^2}} e^{- \frac{(z-\mu)^2}{2 \sigma^2}}である．データ\bm{x}_n = \{ x_1, \cdots, x_n \}に基づくモデルの対数尤度は
\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
である．最尤法によって最尤推定値
\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
が求まり，最大対数尤度
\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*}
と最尤モデル\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}}が求まる．これにより
\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
である．
ここまでは真の分布gが何であっても成り立つ議論であった．以降，真の分布gが正規分布N(\mu_g, \sigma_g^2)であるとしよう．平均対数尤度は，\bm{x}_nを固定したままz \sim N(\mu_g, \sigma_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*}
と計算できる．
補正項，つまり最大対数尤度と平均対数尤度（のn倍）の差は
\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}
となる．期待値をとることで
\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*}
となる．上記の式変形では，
\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 / \chi^2_r] = 1 / (r-2)が成り立つこと（\chi^2_rは自由度rのカイ二乗統計量）と，標本平均と標本分散が独立であることを証明無しに認めて用いた．
数式の議論まとめ
以上をまとめる．正規分布N(\mu_g, \sigma_g^2)に従うデータ\bm{x}_n = \{ x_1, \cdots, x_n \}があるとする．このデータを正規分布モデルで表現することを想定する．最尤推定により推定された最尤モデルf(\cdot | \hat{\mu}, \hat{\sigma}^2)の良さを「最大対数尤度 - 補正項」で評価する．ここで
補正項 = \frac{n}{2 \hat{\sigma}^2} \{ \sigma_g^2 + (\mu_g-\hat{\mu})^2 \} - \frac{n}{2}
であり，その期待値は
である．ちなみにn \to \inftyにおいてb(G) \approx 2であり，これは正規分布モデルのパラメータ数に等しい．
Pythonでのシミュレーション
補正項とその期待値をシミュレーションしてみよう．環境は以下の通り．
Python 3.8.13
numpy 1.23.5
import numpy as np
mu = 3
sigma = 0.1
def correction_calc(mu=0, sigma=0.1, n=100, seed=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
この準備のもと，サンプルサイズnを変えながら補正項を5つずつ計算していく．（ランダムシードを固定している点はご容赦を．）
n=10
n = 10
for i in range(1, 6):
print(correction_calc(mu, sigma, n, i))
-1.4412341801890833
1.289367281291841
1.7921243249806118
4.009671376323563
0.5121271185925691
n=1000
n = 1000
for i in range(1, 6):
print(correction_calc(mu, sigma, n, i))
20.333828816127774
-3.2085759986053564
-8.150289753807993
53.105746883741475
10.831668304333908
この2つの場合だけでも，サンプルサイズnが大きいほど補正項の分散が大きそうなことがわかる．
ここまでは補正項を5つだけ計算していたが，より多くの補正項を計算し，期待値としてその標本平均を求めてみる．
n=10
n = 10
list = []
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 = []
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=1000
n = 1000
list = []
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 = []
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
先ほど見たようにサンプルサイズnが大きいほど標本分散が大きくなっているが，標本平均は理論値2の周りをうろうろしているように見えなくもない．
なお，理論的には期待値E_{G(\bm{x}_n)}[\cdot]をとるが，このシミュレーションでは標本平均を考えているので，もしかしたら見逃している事項があるかもしれない．
今後について
情報量規準の定義に至るには，
\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*}
という分解を行い，D_2 = 0であること，またn \to \inftyにおいて漸近的にD_1 = D_3 = \frac{1}{2} \mathrm{tr} \{I(\bm{\theta}_0) J(\bm{\theta}_0)^{-1}\}（諸々の定義は次回の記事で行う）が成り立つことを示すのが重要となる．次回の記事を書く機会があれば，この辺りをざっと眺めたい（厳密な証明はできないと思うので）．
参考文献
Discussion