『医学のための因果推論』のための統計学入門：9.一般化線形モデル
『医学のための因果推論Ⅰ 一般化線形モデル』を勉強した時のまとめです。
結構難しい本でしたので自分なりの解釈を入れながらまとめました。
一般化線型モデル（GLM）は、分散分析やポアソン回帰などの回帰モデル全般をまとめたものです。
この章では一般化線型モデルを統一的な形式で説明します。
このモデルはパラメトリック分布の一種であるため、大量のデータがあるときには最尤法（最もデータに合ったパラメータを探す方法）を用いることができます。
一般化線型モデルのパラメータは、スコア方程式を数値的に解くことで推定されます。
また統計的な検定や信頼区間も最尤法によって構成することが可能です。
指数型分布族
定義
確率変数 $${Y}$$ が、単一のパラメータ $${\theta}$$ によって決定される確率分布に従い、その確率密度関数または確率関数が次のような形式で表されるとき、この分布を指数型分布族（exponential family）と呼びます。
$$
p(y;θ)=exp[a(y)b(θ)+c(θ)+d(y)]p(y ; \theta) = \exp[a(y) b(\theta) + c(\theta) + d(y)]p(y;θ)=exp[a(y)b(θ)+c(θ)+d(y)]
$$
ここでの $${a(x), b(x), c(x), d(x)}$$ の選び方によって、正規分布や二項分布などが決まります。
逆に言えば、全ての確率分布が扱いやすいわけではないので、指数型分布族に注目することで一般化線型モデルへの拡張が容易になります。
特に、$${a(y)=y}$$ のとき、この分布は正準（canonical）と呼ばれます。
正準という言葉は、「標準」や「基本的な形」という意味合いをもちます。
正規分布や二項分布、ポアソン分布はすべて正準な指数型分布族です。
$${b(\theta)}$$ は自然パラメータと呼ばれます。
自然パラメータは指数型分布族において分布の形や特性を直接決定する基本的な要素です。
$${b(\theta)}$$ は分布の位置やスケールといった重要な特徴を制御するため、このパラメータが指数部分に現れることで確率密度関数は指数型分布族の形を保つことができます。
指数型分布族には便利な性質がいくつかあります。
関数 $${f(x)}$$ の導関数を $${f'(x)}$$と表すとき、$${a(y)}$$ の期待値と分散は次のように書けます。
$$
\begin{gathered} \mathrm{E}[a(Y)] = -\frac{c'(\theta)}{b'(\theta)} \\ \operatorname{Var}[a(Y)] = \frac{b''(\theta) c'(\theta) - c''(\theta) b'(\theta)}{[b'(\theta)]^3} \end{gathered}
$$
指数型分布族の対数尤度関数は$${l(\theta) = a(y) b(\theta) + c(\theta) + d(y)}$$ なので、スコア関数はこれを微分して $${U(\theta) = a(y) b'(\theta) + c'(\theta)}$$ となります。
スコア関数の期待値は$${\mathrm{E}[U(\theta)] = 0}$$ であり、Fisher情報量はスコア関数の分散であるため $${I(\theta) = \operatorname{Var}[U(\theta)] = [b'(\theta)]^2 \operatorname{Var}[a(Y)]}$$となります。
例えば、正規分布の確率密度関数は以下のように変形できます。
$$
\begin{aligned} p(y ; \mu, \sigma^2) &= \frac{1}{\sqrt{2 \pi \sigma^2}} \exp\left[-\frac{1}{2}\left(\frac{y - \mu}{\sigma}\right)^2\right] \\ &= \exp\left[\frac{y \mu}{\sigma^2} - \frac{y^2}{2 \sigma^2} - \frac{\mu^2}{2 \sigma^2} - \frac{1}{2} \log(2 \pi \sigma^2)\right] \end{aligned}
$$
ここで、 $${b(\mu) = \frac{\mu}{\sigma^2},c(\mu) = -\frac{\mu^2}{2 \sigma^2} - \frac{1}{2} \log(2 \pi \sigma^2),d(y) = -\frac{y^2}{2 \sigma^2}}$$ とおけば、正準形の指数型分布族であることがわかります。
以下のようにして期待値や分散も確認できます。
$$
\mathrm{E}(Y) = -\frac{c'(\theta)}{b'(\theta)} = \mu,
$$
$$
\operatorname{Var}(Y) = \frac{b''(\theta) c'(\theta) - c''(\theta) b'(\theta)}{[b'(\theta)]^3} = \sigma^2
$$
一般化線型モデル（GLM）
モデルの構造
これまでは独立同一な指数型分布族の最尤法について述べてきました。
それでは、「平均は共通」という仮定が満たされない場合にこの方法を適用するにはどうしたらよいでしょうか。
一般化線型モデルを提案したNelderとWedderburn（1972）のアイデアは以下のようなものでした。
まず、$${Y_i}$$ の分布が正準（$${a(y)=y}$$）であり、単一のパラメータ $${\theta_i}$$ で決まると仮定します。
また、$${\theta_i}$$ を除いたその他の関数は共通とします。
$$
p(yi;θi)=exp[yib(θi)+c(θi)+d(yi)]p(y_i ; \theta_i) = \exp\left[y_i b(\theta_i) + c(\theta_i) + d(y_i)\right]p(yi;θi)=exp[yib(θi)+c(θi)+d(yi)]
$$
このとき、$${Y_i(i = 1, \ldots, N)}$$の同時分布は次のように書けます。
$$
\begin{aligned} p(y_1, \ldots, y_N) &= \prod_{i=1}^N \exp\left[y_i b(\theta_i) + c(\theta_i) + d(y_i)\right] \\ &= \exp\left[\sum_{i=1}^N y_i b(\theta_i) + \sum_{i=1}^N c(\theta_i) + \sum_{i=1}^N d(y_i)\right] \end{aligned}
$$
このモデルでは、個々のデータ点 $${Y_i}$$ の平均とデザイン行列 $${\boldsymbol{X}_i}$$の関係を次のように表します。
$$
\left[\mathrm{E}(Y_i \mid \boldsymbol{X}_i)\right] = \boldsymbol{X}_i \boldsymbol{\beta}
$$
このような構造を持つ確率分布を一般化線型モデル（GLM）と呼びます。
このとき、個人レベルのパラメータ $${\theta_i}$$ は p 次元のパラメータ $${\beta}$$ によって表され、パラメータの数は $${N>p}$$ である必要があります。
$${g(x)}$$ はリンク関数（link function）と呼ばれ、1対1の単調な変換です。
実際の解析では、リンク関数はデータに合ったものが選択されます。
指数型分布族の特性から、 $${\theta_i, \mathrm{E}(Y_i \mid \boldsymbol{X}_i),\beta}$$ は、関数 $${b(x),c(x),g(x)}$$ を用いて相互に変換することが可能です。
十分統計量と正準リンク関数
GLMの同時分布は積の形になっており、対数尤度関数が次のように和の形式で表されるため、最尤法の計算上便利です
$$
l(\theta_1, \ldots, \theta_N) = \sum_{i=1}^N y_i b(\theta_i) + \sum_{i=1}^N c(\theta_i) + \sum_{i=1}^N d(y_i)
$$
リンク関数としてさまざまな変換が使われますが、特に重要なのは次のようなリンク関数です。
$$
\sum_{i=1}^N y_i b(\theta_i) = \sum_{i=1}^N y_i \boldsymbol{X}_i \beta
$$
このとき、パラメータ $${\beta}$$ の最尤推定量は $${\sum_{i=1}^N y_i \boldsymbol{X}_i}$$ だけに依存し、特にサンプルサイズが小さいときでも計算が安定します。
例)２項分布におけるロジスティック回帰
この $${\sum_{i=1}^N y_i \boldsymbol{X}_i}$$i をGLMの十分統計量と呼び、上記のような性質を満たすリンク関数を正準リンク関数と呼びます。
例えば、正規分布、2項分布、ポアソン分布の正準リンク関数は次のようになります。
$$
g(\mu_i) = \mu_i \text{正規分布}
$$
import numpy as np
import matplotlib.pyplot as plt
# 説明変数 X と応答変数 Y の関係を正規分布で表現
def normal_glm(x, beta0, beta1, sigma):
mu = beta0 + beta1 * x # 恒等リンクによる平均
y = np.random.normal(mu, sigma, size=len(x)) # 正規分布に基づく応答変数の生成
return y, mu
# パラメータ設定
x = np.linspace(-10, 10, 100)
beta0 = 5
beta1 = 0.8
sigma = 2
# データ生成
y, mu = normal_glm(x, beta0, beta1, sigma)
# 図示
plt.figure(figsize=(10, 6))
plt.plot(x, mu, label="Expected Mean (μ = β0 + β1 * X)", color='red', linestyle='--')
plt.scatter(x, y, alpha=0.5, label="Generated Data (Y ~ Normal(μ, σ^2))", color='blue')
plt.xlabel("X")
plt.ylabel("Y")
plt.title("Normal Distribution in GLM with Identity Link")
plt.legend()
plt.grid(True)
plt.show()$$
g(\mu_i) = \log\left(\frac{\mu_i}{1 - \mu_i}\right) \text{2項分布}
$$
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import bernoulli
# 説明変数 X と応答変数 Y の関係を二項分布で表現
def binomial_glm(x, beta0, beta1):
eta = beta0 + beta1 * x # リニア予測子
p = 1 / (1 + np.exp(-eta)) # ロジットリンク関数で成功確率を計算
y = bernoulli.rvs(p) # 二項分布（ベルヌーイ分布）に基づく応答変数の生成
return y, p
# パラメータ設定
x = np.linspace(-5, 5, 100) # 説明変数の範囲
beta0 = 0.0 # 切片
beta1 = 1.0 # 傾き
# データ生成
y, p = binomial_glm(x, beta0, beta1)
# 図示
plt.figure(figsize=(10, 6))
plt.plot(x, p, label="Success Probability (p = 1 / (1 + exp(-(β0 + β1 * X))))", color='purple', linestyle='--')
plt.scatter(x, y, alpha=0.5, label="Generated Data (Y ~ Bernoulli(p))", color='blue')
plt.xlabel("X")
plt.ylabel("Y (Success Probability)")
plt.title("Binomial Distribution in GLM with Logit Link")
plt.legend()
plt.grid(True)
plt.show()
$$
g(\mu_i) = \log(\mu_i) \text{ポアソン分布}
$$
# 説明変数 X と応答変数 Y の関係をポアソン分布で表現
def poisson_glm(x, beta0, beta1):
eta = beta0 + beta1 * x # リニア予測子
mu = np.exp(eta) # 対数リンク関数で期待値を計算
y = np.random.poisson(mu) # ポアソン分布に基づく応答変数の生成
return y, mu
# パラメータ設定
x = np.linspace(0, 5, 50)
beta0 = 1
beta1 = 0.8
# データ生成
y, mu = poisson_glm(x, beta0, beta1)
# 図示
plt.figure(figsize=(10, 6))
plt.plot(x, mu, label="Expected Mean (μ = exp(β0 + β1 * X))", color='green', linestyle='--')
plt.scatter(x, y, alpha=0.5, label="Generated Data (Y ~ Poisson(μ))", color='orange')
plt.xlabel("X")
plt.ylabel("Y")
plt.title("Poisson Distribution in GLM with Log Link")
plt.legend()
plt.grid(True)
plt.show()正準リンク関数を用いたモデルは、他のリンク関数に比べて小さいサンプルサイズでの最尤推定量の挙動が良いとされています。
パラメータの推定
一般化線型モデル（GLM）では、パラメータベクトル $${\boldsymbol{\beta} = (\beta_1, \ldots, \beta_p)^T}$$ の最尤推定量を求めるために、これまでと同様の方針で対数尤度を偏微分して計算します。
具体的には、各パラメータ $${\beta_j}$$ に対してスコア関数 $${U_j = \frac{\partial l(\theta_1, \ldots, \theta_N)}{\partial \beta_j}}$$ を求め、すべてのパラメータに対してスコア方程式 $${\boldsymbol{U}(\boldsymbol{\beta}) = 0}$$ を解きます。
具体的には、スコア関数の要素は次の式で表されます。
$$
U_j = \sum_{i=1}^N \frac{y_i - \mu_i}{\operatorname{Var}(Y_i \mid \boldsymbol{X}_i)} x_{ij} \frac{\partial \mu_i}{\partial \eta_i}
$$
ここで、期待値 $${\mathrm{E}(Y_i \mid \boldsymbol{X}_i) = \mu_i}$$ であり、リンク関数の関係式 $${g(\mu_i) = \eta_i}$$ を用います。
Fisher情報行列 $${\boldsymbol{I} = \mathrm{E}}$$ はスコア関数を微分して求めることができ、行列の各要素は次のように表されます。
$$
I_{jk} = \sum_{i=1}^N \frac{x_{ij} x_{ik}}{\operatorname{Var}(Y_i \mid \boldsymbol{X}_i)} \left( \frac{\partial \mu_i}{\partial \eta_i} \right)^2
$$
この情報行列を計算する際は、対角行列 $${\boldsymbol{W}}$$ を定義すると便利です。
対角成分 $${w_{ii} = \frac{1}{\operatorname{Var}(Y_i \mid \boldsymbol{X}_i)} \left( \frac{\partial \mu_i}{\partial \eta_i} \right)^2}$$ で構成された $${N \times N}$$ 行列として $${\boldsymbol{W}}$$ を定義すると、Fisher情報行列は次のように表されます。
$$
\boldsymbol{I} = \boldsymbol{X}^T \boldsymbol{W} \boldsymbol{X}
$$
最後に、パラメータ $${\beta_j}$$ のスコア関数がどのように導かれるのかを見ていきます。
一般化線型モデルにおいて、対数尤度関数は各データ点の和として書けるため、$${i}$$ 番目のデータが尤度関数にどのように寄与するかは $${l_i = y_i b(\theta_i) + c(\theta_i) + d(y_i)}$$ で表すことができます。
スコア関数は、微分の連鎖法則を利用して次のように計算されます。
$$
U_j = \sum_{i=1}^N \frac{\partial l_i}{\partial \beta_j} = \sum_{i=1}^N \frac{\partial l_i}{\partial \theta_i} \frac{\partial \theta_i}{\partial \mu_i} \frac{\partial \mu_i}{\partial \beta_j}
$$
これを導出するために、各部分をひとつひとつ微分して計算します。
目標 は $${b\left(\theta_i\right), c\left(\theta_i\right), d(y) }$$を具体的な形に置き換えることです.
$$
\begin{aligned}\frac{\partial l_i}{\partial \theta_i}=y_i b^{\prime}\left(\theta_i\right)+c^{\prime}\left(\theta_i\right)=b^{\prime}\left(\theta_i\right)\left(y_i-\mu_i\right)\end{aligned}
$$
となり、 $${\begin{aligned}\mu_i=-\frac{c^{\prime}\left(\theta_i\right)}{b^{\prime}\left(\theta_i\right)}\end{aligned}}$$ を微分すると
$$
\begin{aligned}\frac{\partial \mu_i}{\partial \theta_i}=\frac{-c^{\prime \prime}\left(\theta_i\right)}{b^{\prime}\left(\theta_i\right)}+\frac{c^{\prime}\left(\theta_i\right) b^{\prime \prime}\left(\theta_i\right)}{b^{\prime}\left(\theta_i\right)^2}=b^{\prime}\left(\theta_i\right) \operatorname{Var}\left(Y_i \mid \boldsymbol{X}_i\right)\end{aligned}
$$
となることを利用すれば、
$$
\begin{aligned}\frac{\partial \theta_i}{\partial \mu_i}=\left(\frac{\partial \mu_i}{\partial \theta_i}\right)^{-1}=\frac{1}{b^{\prime}\left(\theta_i\right) \operatorname{Var}\left(Y_i \mid \boldsymbol{X}_i\right)}\end{aligned}
$$
となります。
最後の部分は $${\eta_i=X_i \beta}$$ の偏微分を考えることで、
$$
\begin{aligned}\frac{\partial \mu_i}{\partial \beta_j}=\frac{\partial \mu_i}{\partial \eta_i} \frac{\partial \eta_i}{\partial \beta_j}=\frac{\partial \mu_i}{\partial \eta_i} x_{i j}\end{aligned}
$$
とります。
個々のデータのスコア関数への貢献は、３つの偏微分の積で決まりますから、総和をとれば
$$
\begin{aligned}U_j=\sum_{i=1}^N \frac{y_i-\mu_i}{\operatorname{Var}\left(Y_i \mid \boldsymbol{X}_i\right)} x_{i j} \frac{\partial \mu_i}{\partial \eta_i}\end{aligned}
$$
が導かれます。
Newton-Raphson法によるスコア方程式の解法
スコア方程式を解くにはどのようにすれば良いでしょうか。
一般化線型モデルのスコア関数は通常、パラメータ$${\beta}$$ が非線型な関数であるため、簡単に解を求めることができません。
そのため、Newton-Raphson法などの計算アルゴリズムを用いることが一般的です。
Newton-Raphson法は、ある関数 $${f(x)}$$ とその導関数 $${f′(x)}$$、初期値 $${x^{(a=0)}}$$ が与えられたときに方程式 $${f(x)=0}$$ の解を求めるアルゴリズムです。
以下の反復計算によって$${x}$$ を更新します。
$$
x^{(a+1)} = x^{(a)} - \frac{f(x^{(a-1)})}{f'(x^{(a-1)})}
$$
$${a}$$ は反復回数であり、左辺の$${x^{(a+1)}}$$ は$${a}$$回目の計算結果です。
Newton-Raphson法では、反復計算を停止するための基準が必要です。
例えば正の実数$${\varepsilon > 0}$$ を指定して、反復ごとに関数の値が十分ゼロに近くなるか（$${|f(x^{(a)})| < \varepsilon}$$）、または$${x}$$ の変化が小さくなるか（$${|x^{(a)} - x^{(a-1)}| < \varepsilon}$$）を確認するなどの手続きが用いられます。
コンピュータでNewton-Raphson法を実行する際、計算精度は浮動小数点の桁数に依存しているため、計算ごとに丸め誤差が生じます。
そのため実用上は関数$${f(x) = 0}$$の厳密な解を求めるのではなく、$${f(x) }$$が十分にゼロに近い$${x}$$ を探すことが実用的な方法とされています。
import numpy as np
import matplotlib.pyplot as plt
# 関数 f(x) とその導関数 f'(x) の定義
def f(x):
return x**3 - 2 * x + 2
def f_prime(x):
return 3 * x**2 - 2
# Newton-Raphson法による解法（ダンピング因子付き）
def newton_raphson_damped(x0, alpha=0.5, tol=1e-6, max_iter=20):
x_values = [x0] # 反復ごとの x の値を保存
for _ in range(max_iter):
fx = f(x0)
fpx = f_prime(x0)
# 導関数がゼロに近い場合、更新を停止して二分法に切り替え
if abs(fpx) < 1ｒe-6:
print("導関数がゼロに近いため、Newton-Raphson法を中止し二分法に切り替える")
return bisection_method(-2, 2, tol=tol) # 解が存在しそうな範囲を指定
# Newton-Raphson法による更新（ダンピング因子 alpha を掛ける）
x_new = x0 - alpha * fx / fpx
x_values.append(x_new)
# 収束判定
if abs(f(x_new)) < tol:
break
x0 = x_new
return x_values
# 二分法による解法
def bisection_method(a, b, tol=1e-6, max_iter=100):
if f(a) * f(b) > 0:
print("指定した範囲に解が含まれていません。範囲を変更してください。")
return []
x_values = []
for _ in range(max_iter):
c = (a + b) / 2
x_values.append(c)
# 収束判定
if abs(f(c)) < tol or abs(b - a) < tol:
break
if f(a) * f(c) < 0:
b = c
else:
a = c
return x_values
# 初期値とパラメータ
x0 = -1.5 # 初期値を調整
alpha = 0.5 # ダンピング因子
x_values = newton_raphson_damped(x0, alpha=alpha)
# プロット設定
x_range = np.linspace(-2, 2, 400)
y_range = f(x_range)
plt.figure(figsize=(10, 6))
plt.plot(x_range, y_range, label="f(x) = x^3 - 2x + 2", color="blue")
plt.axhline(0, color="black", linestyle="--", linewidth=0.7)
# Newton-Raphson法の各反復を図示
for i in range(1, len(x_values)):
x_prev = x_values[i - 1]
x_curr = x_values[i]
plt.plot([x_prev, x_prev], [0, f(x_prev)], color="red", linestyle="--", linewidth=0.5)
plt.plot([x_prev, x_curr], [f(x_prev), 0], color="green", linestyle="--", linewidth=0.5)
plt.scatter(x_curr, 0, color="orange", zorder=5)
# 最終結果表示
plt.scatter(x_values[-1], 0, color="red", label="Approximate Root", zorder=5)
plt.xlabel("x")
plt.ylabel("f(x)")
plt.title("Damped Newton-Raphson Method with Bisection Backup")
plt.legend()
plt.grid(True)
plt.show()図は関数 $${f(x) = x^3 - 2x + 2}$$ とその導関数 $${f'(x) = 3x^2 - 2}$$ の$${x=0}$$の解において、Newton-Raphson法を初期値 $${x_0 = 1.5}$$ から始め、反復的に $${x}$$ の値を更新し、近似解に収束させたものです。
情報量規準とモデルの選択
赤池情報量規準（AIC）
尤度比検定は対数尤度を比較してどちらのモデルがデータに適しているかを判定する方法です。
しかしモデルが真に適合しているかどうかだけでなく、モデルがデータにどの程度合っているかも評価したい場面があります。
特にGLMは異なるデザイン行列を持つ無数のモデル候補があります。
このとき、どのモデルがデータに最もよく適合しているかを評価するための指標として、赤池情報量規準（AIC）が用いられます（Akaike, 1973）。
AICは対数尤度関数に最尤推定量を代入した値に基づきますが、複雑なモデルを選びがちになるため、バイアス補正がなされています。
具体的には、次の式でAICを定義します。
$$
\text{AIC for Model A} = -2 l(\widehat{\boldsymbol{\beta}}_A) + 2q, \quad \text{AIC for Model B} = -2 l(\widehat{\boldsymbol{\beta}}_B) + 2p
$$
ここで、$${\widehat{\boldsymbol{\beta}}_Aβ}$$ と $${\widehat{\boldsymbol{\beta}}_B}$$ はモデルAとBそれぞれの最尤推定量であり、$${q}$$ と $${p}$$ はそれぞれのパラメータ数です。
AICが小さいほど予測性能が良いと判断されます。
AICとKullback-Leibler情報量の関係
モデルがデータにどれほど適合しているかを測る別の指標のひとつとして、Kullback-Leibler情報量（KL情報量）があります。
このKL情報量は、データに適合させたモデル $${p(y ; \boldsymbol{\beta})}$$ と真のモデル $${q(y)}$$ の間の「距離」を測る指標として、次のように定義されます。
$$
KL = \int \log\left(\frac{q(y)}{p(y ; \boldsymbol{\beta})}\right)
$$
AICを小さくすることが真のモデルにより近いモデルを選択することにつながるのは、このKL情報量が小さいモデルを選ぶという意味になります。
KL情報量は、次のように2つの項に分けることができます。
片方の項はデータに適合させたモデル$${p(y ; \boldsymbol{\beta})}$$に依存せず、もう片方の項がモデル$${p(y ; \boldsymbol{\beta})}$$に依存する部分です。
$$
KL = \int \log[q(y)] q(y) dy - KL(\widehat{\boldsymbol{\beta}})
$$
ここで、$${KL(\boldsymbol{\beta}) = \mathrm{E}[\log(p(y ; \boldsymbol{\beta}))]}$$ とおきます。
この定義に基づくと、AICを最小化することは期待値としてのKL情報量を最小化する、つまり真のモデルにより近いモデルを選ぶことに相当します。
右側のプロット： AICの値が候補モデルの平均 𝜇 μ に対してどのように変化するかを示しています。 真のモデルに最も近いと考えられるモデルは、AICが最小のモデルであることを示しています。
図から、AICの最小化がKL情報量の最小化に相当する、すなわち対数尤度が最大のモデルが結果としてAICの値を最小にし、そのモデルが真のモデルに最も近いという概念を示しています。
import numpy as np
import matplotlib.pyplot as plt
# 真のモデル q(y) を標準正規分布と仮定
def true_model(y):
return (1 / np.sqrt(2 * np.pi)) * np.exp(-0.5 * y**2)
# 候補モデル p(y; β) を平均 μ の正規分布と仮定
def candidate_model(y, mu):
return (1 / np.sqrt(2 * np.pi)) * np.exp(-0.5 * (y - mu)**2)
# 対数尤度を計算する関数
def log_likelihood(y, p):
return np.sum(np.log(p)) * (y[1] - y[0])
# yの範囲と真の分布 q(y)
y = np.linspace(-4, 4, 100)
q_y = true_model(y)
# 候補モデルの異なる μ に対する対数尤度と AIC の計算
mu_values = np.linspace(-2, 3, 50)
log_likelihood_values = [log_likelihood(y, candidate_model(y, mu)) for mu in mu_values]
aic_values = [-2 * ll + 2 for ll in log_likelihood_values] # AIC = -2 * log(likelihood) + 2 * p (ここではパラメータ数p=1)
# 図示
plt.figure(figsize=(14, 6))
# 対数尤度と AIC の関係
plt.subplot(1, 2, 1)
plt.plot(mu_values, log_likelihood_values, color='blue', label="Log-Likelihood", linestyle='-', marker='o')
plt.xlabel(r'Candidate model mean $\mu$')
plt.ylabel("Log-Likelihood")
plt.title("Log-Likelihood of Candidate Models")
plt.axvline(0, color='red', linestyle='--', label="True Model ($\mu = 0$)")
plt.legend()
plt.grid(True)
plt.subplot(1, 2, 2)
plt.plot(mu_values, aic_values, color='green', label="AIC", linestyle='-', marker='o')
plt.xlabel(r'Candidate model mean $\mu$')
plt.ylabel("AIC")
plt.title("AIC of Candidate Models")
plt.axvline(0, color='red', linestyle='--', label="True Model ($\mu = 0$)")
plt.legend()
plt.grid(True)
plt.suptitle("Relationship Between Log-Likelihood, AIC, and True Model")
plt.tight_layout(rect=[0, 0, 1, 0.95])
plt.show()AICとKL情報量の関係をわかりやすくするために、次に-2倍を掛けたものを考えます。
以下のように期待値の観点からAICを考えると、モデルのKL情報量が小さい（すなわち真のモデルに近い）ものが選ばれることがわかります。
$$
\begin{aligned} -2N \times KL(\widehat{\boldsymbol{\beta}}) & \approx -2N \times KL(\boldsymbol{\beta}) + (\widehat{\boldsymbol{\beta}} - \boldsymbol{\beta})^T \boldsymbol{I}(\boldsymbol{\beta})(\widehat{\boldsymbol{\beta}} - \boldsymbol{\beta}) \\ & \approx -2 \mathrm{E}[l(\boldsymbol{\beta})] + (\widehat{\boldsymbol{\beta}} - \boldsymbol{\beta})^T \boldsymbol{I}(\boldsymbol{\beta})(\widehat{\boldsymbol{\beta}} - \boldsymbol{\beta}) \\ & \approx -2 \mathrm{E}[l(\widehat{\boldsymbol{\beta}})] + 2(\widehat{\boldsymbol{\beta}} - \boldsymbol{\beta})^T \boldsymbol{I}(\boldsymbol{\beta})(\widehat{\boldsymbol{\beta}} - \boldsymbol{\beta}) \\ & \approx -2 \mathrm{E}[l(\widehat{\boldsymbol{\beta}})] + 2 \chi_p^2 \end{aligned}
$$
この式の両辺の期待値を取ると、$${\mathrm{E}(\chi_p^2) = p}$$ より次のように整理できます：
$$
-2N \times \mathrm{E}[KL(\widetilde{\boldsymbol{\beta}})] \approx -2 \mathrm{E}[l(\widetilde{\boldsymbol{\beta}})] + 2p
$$
AICの期待値は次のように表され、モデルの選択基準として用いることができます。
$$
\mathrm{E}(AIC) = -2 \mathrm{E}[l(\widehat{\boldsymbol{\beta}})] + 2p
$$
これにより、AICとKL情報量の期待値には近似的な関係があることがわかります。
あるモデルのAICが小さい場合、そのモデルがKL情報量の観点から真のモデルに近いことを意味します。
ただしこの議論はあくまで期待値に基づいているため、AICのばらつきが考慮されていない点に注意が必要です。
参考文献
『医学のための因果推論Ⅰ 一般化線形モデル』