統計モデルの比較

このページでは、モデルの「良さ」の比較について考えてみる。

パラメーター最尤値と尤度のバイアス

ある確率分布 $g(x)$ に従っている独立なサンプル $x_1,x_2,\cdots ,x_N$ が得られたとしよう。 そのデータを、$k$個のパラメータの組 ${\theta }_1,{\theta }_2,\cdots ,{\theta }_k$ を持つような分布関数 $f(x;\{\theta \})$ で表現したい。 以下では、混乱が生じない限り、パラメータの組を単に $\theta$ と表記することにする。

そのための自然な発想として、対数尤度の平均 $$\begin{array}{}\text{(1)}& L(\theta ;X)=\frac{1}{N}\sum _{\ell =1}^N\log f(x_{\ell };\theta )\end{array}$$ を最大化するような $\theta$ を選べば良さそうである。すなわち $$\begin{array}{}\text{(2)}& \frac{\partial }{\partial {\theta }_i}L(\hat{\theta };X)=0\end{array}$$ となるようなパラメータの組 $\hat{\theta }$ を選べばよさそうだ。

ところが、データは、サンプリングの都度異なるので、最尤推定したパラメータも、確率的に変動する。 $f$という分布（モデル）を仮定した下での「真の」パラメータは、 $$\begin{array}{}\text{(3)}& L(\theta ;Y)=E_Y[\log f(y;\theta )]=\int g(y)\log f(y;\theta )dy\end{array}$$ を最大化するような値であって、 $$\begin{array}{}\text{(4)}& \frac{\partial }{\partial {\theta }_i}L({\theta }^{\ast };Y)=0\end{array}$$ となるようなパラメータの組を ${\theta }^{\ast }$ と書くことにする。 そして、$N\to \mathrm{\infty }$ の極限では、両者は一致する（そのような場合を考える）。

では、$L(\hat{\theta };X)$ ができるだけ大きくなるようなモデル $f(x)$ が最善のモデルと言えるのだろうか？ 統計理論の教えるところによれば、分布関数 $g(y)$ と $f(x)$ の隔たり具合（「距離」）は、 カルバック・ライブラー情報量(Kullback-Leibler information) $$\begin{array}{r}KL(g;f)=\int g(x)\log \left(\frac{g(x)}{f(x)}\right)dx\\ =\int g(x)\log g(x)dx-\int g(x)\log f(x)dx\end{array}$$ を客観的な指標として比較すべきである（$\text{KL}\ge 0$で、両者が一致する場合のみ $\text{KL}=0$となる）。 そうすると、最大化しなければならないのは、上式の第２項に対応する $$\begin{array}{}\text{(5)}& L(\hat{\theta };Y)=\int g(\hat{\theta };y)\log f(\hat{\theta };y)dy=E_Y[\log f(y;\hat{\theta })]\end{array}$$ ということになる。

一方、この $L(\hat{\theta };Y)$ と、サンプル平均として求めた尤度 $L(\hat{\theta };X)$ との間には、 $N\to \mathrm{\infty }$ でも0にならない偏差（バイアス）が存在することが知られている。

パラメーターのばらつき具合

サンプルから求めたパラメータ$\hat{\theta }$は、真の分布に対して尤度を最大とするパラメータ値${\theta }^{\ast }$の周りでどのようにばらついているのだろうか。

平均対数尤度の微分を、${\theta }^{\ast }$ の周りで展開すると、以下のようになる。 $$\begin{array}{}\text{(6)}& \frac{\partial }{\partial {\theta }_i}L(\hat{\theta };X)=\frac{\partial }{\partial {\theta }_i}L({\theta }^{\ast };X)+\sum _j^k\frac{{\partial }^{2}L({\theta }^{\ast };X)}{\partial {\theta }_j\partial {\theta }_i}({\hat{\theta }}_j-{\theta }_j^{\ast })+\cdots \end{array}$$ ここで、定義から上式の左辺は0である。

以下で用いるため、行列 $A$ の各成分を $$A_{ij}=-\frac{{\partial }^2\log f(x;{\theta }^{\ast })}{\partial {\theta }_j\partial {\theta }_i}$$ とおいて、 $$\frac{{\partial }^{2}L({\theta }^{\ast };X)}{\partial {\theta }_j\partial {\theta }_i}=-E_X\left[A\right]$$ と表しておく。また、ベクトル $\mathit{b}$ の成分 $$b_i=\frac{\partial }{\partial {\theta }_i}L({\theta }^{\ast };X)$$ とおき、パラメータもベクトル表記 $$\mathit{\theta }=({\theta }_1,{\theta }_2,\cdots ,{\theta }_k{)}^{\mathrm{\top }}$$ すると (6)式は $$\begin{array}{}\text{(7)}& E_X[A](\hat{\mathit{\theta }}-{\mathit{\theta }}^{\ast })=\mathit{b}\end{array}$$ と書ける。

ここで、 $$\begin{array}{}\text{(8)}& b_i=\frac{1}{N}\sum _{\ell =1}^N\frac{\partial \log f(x_{\ell };{\mathit{\theta }}^{\ast })}{\partial {\theta }_i}\end{array}$$ は、$N$個の独立な確率変数の和の形をしているので、中心極限定理から正規分布に漸近することが分かる。 その平均 $E_X[b_i]$ は0で、分散・共分散行列の $i,j$ 成分は $$\begin{array}{}\text{(9)}& E_X[b_i{b}_j]=\frac{1}{N}{E}_X\left[\frac{\partial \log f(x;{\mathit{\theta }}^{\ast })}{\partial {\theta }_i}\frac{\partial \log f(x;{\mathit{\theta }}^{\ast })}{\partial {\theta }_j}\right]=\frac{1}{N}{E}_X[B]\end{array}$$ となる。 ここで、行列$B$の成分を $$B_{ij}=\frac{\partial \log f(x;{\mathit{\theta }}^{\ast })}{\partial {\theta }_i}\frac{\partial \log f(x;{\mathit{\theta }}^{\ast })}{\partial {\theta }_j}$$ と置いた。

(7)式から $$\hat{\mathit{\theta }}-{\mathit{\theta }}^{\ast }=E_X[A{]}^{-1}\mathit{b}$$ であるから、$\hat{\mathit{\theta }}-{\mathit{\theta }}^{\ast }$ の分散・共分散行列は $$\begin{array}{}\text{(10)}& E_X[(\hat{\mathit{\theta }}-{\mathit{\theta }}^{\ast })(\hat{\mathit{\theta }}-{\mathit{\theta }}^{\ast }{)}^{\mathrm{\top }}]=E_X[A{]}^{-1}{E}_X[\mathit{b}{\mathit{b}}^{\mathrm{\top }}]{(E_X[A{]}^{-1})}^{\mathrm{\top }}=\frac{1}{N}{E}_X[A{]}^{-1}{E}_X[B]E_X[A{]}^{-1}\end{array}$$ となる。

$N$が十分大きければ、$X$についての平均は$Y$についてのそれに代えられるので（大数の法則） $$E_X[(\hat{\mathit{\theta }}-{\mathit{\theta }}^{\ast })(\hat{\mathit{\theta }}-{\mathit{\theta }}^{\ast }{)}^{\mathrm{\top }}]\approx \frac{1}{N}{E}_Y[A{]}^{-1}{E}_Y[B]E_Y[A{]}^{-1}$$

対数尤度のバイアス量の見積もり

データから推定したパラメータの組 $\{\hat{\theta }\}$ から導出される尤度の差の平均量 $$\begin{array}{}\text{(11)}& C=E_X[L(\hat{\theta }(X);Y)-L(\hat{\theta }(X);X)]\end{array}$$ を考える。ここで、$\hat{\theta }$はサンプル$X$から計算される量であることめ明示するため $\hat{\theta }(X)$ と表記した。

ここで、(11)式を $$\begin{array}{}\text{(12)}& \begin{array}{r}{E}_X[L(\hat{\theta }(X);Y)-L({\theta }^{\ast };Y)]\\ +E_X[L({\theta }^{\ast };Y)-L({\theta }^{\ast };X)]\\ +E_X[L({\theta }^{\ast };X)-L(\hat{\theta }(X);X)]\end{array}\end{array}$$ と書き直し、１行目を $C_1$、2行目を $C_2$, 3行目を $C_3$ とする。 言うまでもなく、$C=C_1+C_2+C_3$ である。

$C_1$の評価

対数尤度を${\theta }^{\ast }$の周りでテイラー展開すると $$\begin{array}{}\text{(13)}& \begin{array}{r}{E}_Y[\log f(y;\hat{\theta }(X))]=E_Y[\log f(y;{\theta }^{\ast })]+\sum _{i=1}^k\frac{\partial E_Y[\log f(y;{\theta }^{\ast })]}{\partial {\theta }_i}(\theta (X{)}_i-{\theta }_i^{\ast })\\ +\frac{1}{2}\sum _i\sum _j\frac{{\partial }^2{E}_Y[\log f(y;{\theta }^{\ast })]}{\partial {\theta }_i\partial {\theta }_j}(\theta (X{)}_i-{\theta }_i^{\ast })(\theta (X{)}_j-{\theta }_j^{\ast })+\cdots \end{array}\end{array}$$ を得る。 ここで、${\theta }^{\ast }$の定義より、右辺２項目（一階微分の項）は0となる。 上式を、対称行列 $$\begin{array}{}\text{(14)}& J=-E_Y\left[\frac{{\partial }^2\log f(y;{\theta }^{\ast })}{\partial {\theta }_i\partial {\theta }_j}\right]\end{array}$$ を使って整理すると、 $$\begin{array}{}\text{(15)}& L(\hat{\theta }(X);Y)-L({\theta }^{\ast };Y)\approx -\frac{1}{2}(\hat{\mathit{\theta }}-{\mathit{\theta }}^{\ast }{)}^{\mathrm{\top }}J(\hat{\mathit{\theta }}-{\mathit{\theta }}^{\ast })\end{array}$$ となる。

右辺の平均を求めるために、まず、以下の計算をおこなう。 $$\begin{array}{}\text{(16)}& \begin{array}{r}{E}_X[(\hat{\mathit{\theta }}-{\mathit{\theta }}^{\ast }{)}^{\mathrm{\top }}J(\hat{\mathit{\theta }}-{\mathit{\theta }}^{\ast })]=E_X[\text{trace}(J(\hat{\mathit{\theta }}-{\mathit{\theta }}^{\ast })(\hat{\mathit{\theta }}-{\mathit{\theta }}^{\ast }{)}^{\mathrm{\top }})]\\ =\text{trace}\left(J{E}_X[(\hat{\mathit{\theta }}-{\mathit{\theta }}^{\ast })(\hat{\mathit{\theta }}-{\mathit{\theta }}^{\ast }{)}^{\mathrm{\top }}]\right)\\ \approx \frac{1}{N}\text{trace}(J{E}_Y[A{]}^{-1}{E}_Y[B]E_Y[A{]}^{-1})\end{array}\end{array}$$

ここで、$E_Y[A]=J$、すなわち $E_Y[A{]}^{-1}=J^{-1}$ に他ならないこと、そして、 $$\begin{array}{}\text{(17)}& I=E_Y[B]=E_Y\left[\frac{\partial \log f(y;{\mathit{\theta }}^{\ast })}{\partial {\theta }_i}\frac{\partial \log f(y;{\mathit{\theta }}^{\ast })}{\partial {\theta }_j}\right]\end{array}$$ と改めて置き直すことで、(16)式は $$\begin{array}{}\text{(18)}& =\frac{1}{N}\text{trace}\left(I{J}^{-1}\right)\end{array}$$ となる。

さらに、$I$と$J$はFisher情報行列の異なる表現に他ならない（$I=J$）から、結局 $$\begin{array}{}\text{(19)}& =\frac{1}{N}\text{trace}(E)=\frac{k}{N}\end{array}$$ を得る($E$は単位行列）。

以上をまとめると、$N\to \mathrm{\infty }$ での $C_1$ の見積もりとして $$\begin{array}{}\text{(20)}& C_1=E_X[L(\hat{\theta }(X);Y)-L({\theta }^{\ast };Y)]\approx -\frac{k}{2N}\end{array}$$ が得られた。

$C_2$の評価

すでに、これまでの計算でも仮定されていたとおり、$X$は$Y$からのサンプルであるから、$N\to \mathrm{\infty }$ では $$\begin{array}{}\text{(21)}& C_2=E_X[L({\theta }^{\ast };Y)-L({\theta }^{\ast };X)]\to 0\end{array}$$ となる。

$C_3$の評価

$C_1$と同様に評価にすると、 $$\begin{array}{}\text{(22)}& C_3=E_X[L({\theta }^{\ast };X)-L(\hat{\theta }(X);X)]\approx -\frac{k}{2N}\end{array}$$ を示すことができる（確認せよ）。

以上から、全体的なバイアスの量 $$\begin{array}{}\text{(23)}& C=C_1+C_2+C_3=-\frac{k}{N}\end{array}$$ が得られた。

赤池情報量基準（AIC)

前節までの結果から、 十分サンプル数が多く、各サンプルが独立であるような場合、サンプルから推定したパラメータから計算した対数尤度にはバイアスが含まれており、 その母平均は $$\begin{array}{}\text{(24)}& L(\hat{\theta }(X);Y)=L(\hat{\theta }(X);X)-\frac{k}{N}=\frac{1}{N}\{\sum _{\ell =1}^N\log f(x_{\ell };\hat{\theta })-k\}\end{array}$$ と見積もられる。ここで、$k$はパラメータの個数である。

尤度の最大値 $\hat{\theta }$ を求め、それを用いて$X$から尤度を計算する際に$k$個の拘束条件を課していることになる。 それら$k$個の自由度が減った分だけの補正が生じ、その値を（いくつかの仮定の下に）評価すると、上記のような結果になる、というわけである。

上式のカッコ内の量を-2倍した量は、赤池情報量基準(AIC)と呼ばれる： $$\begin{array}{}\text{(25)}& \text{AIC}=-2\sum _{\ell =1}^N\log f(x_{\ell };\hat{\theta })+2k\end{array}$$ AICが小さいほど、「適切な」モデルであると考えられる。 パラメータの数を増やせば、それだけデータをよくフィッティングできる一方で、AICでは$2k$だけのペナルティが生じるので、 尤度の最大値とのバランスの中で、適切なモデル ($f(\cdot ;\theta )$) を決めるための判断材料となる。

現実のデータに適用するに際しては、いくつかの改良版の基準が提案されているので、ここから先は専門書などに当たること。

SciPyを使った連続分布のAICの計算

SciPyを使うと、代表的な連続確率分布について、比較的簡単に統計量を計算することができる。 以下のPythonコードは、COVID-19のserial lengthの分析で使われているデータを用い、いくつかの確率分布についてAICを比較した例である。

パラメータの最尤推定はfit()関数を使って行うことができる。 変量Xは (x - loc)/scale のように、パラメータlocとscaleによって変換されるが、 フィッティングの際、それらをある値に固定したい場合は、それぞれ、floc=値、fscale=値のように指定すればよい。 各分布についてのパラメータの個数や意味については scipy.statsのマニュアル を参照のこと。

サンプルに対する確率の対数を求める stats.poisson.logpmf() 関数や stats.norm.logpdf() は、 配列データを与えると配列を返すので、numpy.sum()でその総和を取って対数尤度を求めている。

from scipy import stats
import numpy as np

length=[3, 2, 4, 2, 2, 1, 4, 2, 10, 4, 3, 4, 4, 3, 3, 9, 9, 3, 4, 3, 9, 4, 5, 9, 7, 7, 3, 7]

mu = np.mean(length)
logL_poisson = np.sum(stats.poisson.logpmf(length, mu=mu))
aic_poisson = -2*logL_poisson + 2*1
print('AIC(POISSON)=',aic_poisson)

params_norm = stats.norm.fit(length)
logL_norm = np.sum( stats.norm.logpdf(length, loc=params_norm[0], scale=params_norm[1]) )
aic_norm = -2*logL_norm + 2*2
print('AIC(NORMAL)=',aic_norm)

params_gamma = stats.gamma.fit(length,floc=0)
logL_gamma = np.sum( stats.gamma.logpdf(length, params_gamma[0], loc=params_gamma[1], scale=params_gamma[2]) )
aic_gamma = -2*logL_gamma + 2*2
print('AIC(GAMMA)=',aic_gamma)

params_lognorm = stats.lognorm.fit(length,floc=0)
logL_lognorm = np.sum( stats.lognorm.logpdf(length, params_lognorm[0], loc=params_lognorm[1], scale=params_lognorm[2]) )
aic_lognorm = -2*logL_lognorm + 2*2
print('AIC(LOGNORM)=',aic_lognorm)

params_weibull = stats.weibull_min.fit(length,floc=0)
logL_weibull = np.sum( stats.weibull_min.logpdf(length, params_weibull[0], loc=params_weibull[1], scale=params_weibull[2]) )
aic_weibull = -2*logL_weibull + 2*2
print('AIC(WEIBULL)=',aic_weibull)

　練習：文長のモデルの比較

夏目漱石の「こころ」の文の長さのデータをこちらに用意した。 文の長さが確率的に生成されていると仮定（単純化）した場合、どのようなどのような確率分布で表現するのがよいか、検討してみなさい。

最小二乗法による回帰曲線

データ点の組 $(x_i,y_i)$ を表すモデルが、パラメータを含む関数 $h(x;\theta )$ と、平均0、分散 ${\sigma }^2$ の正規分布に従うゆらぎ $\xi$ を使って、 $$y_i=h(x_i;\theta )+\xi$$ のように表されるとしよう。

データ点の数を$N$とすると、モデルからの残差の二乗和 $$\sum _i^N(y_i-h(x_i;\theta ){)}^2$$ を最小化するようなパラメータを$\hat{\theta }$を与えると、対数尤度を最大化できる（最小二乗法）。 そのときの対数尤度は $$L(\hat{\theta },\hat{\sigma };\{x,y\})=-\frac{N}{2}\log (2\pi {\hat{\sigma }}^2)-\frac{1}{2{\hat{\sigma }}^2}\sum _{i=1}^N(y_i-h(x_i;\hat{\theta }){)}^2$$ となる。

ここで、分散の最尤推定値は $${\hat{\sigma }}^2=\frac{1}{N}\sum _{i=1}^N(y_i-h(x_i;\theta ){)}^2$$ になるから、最大対数尤度は、 $$L(\hat{\theta },\hat{\sigma })=-\frac{1}{2}[\log (2\pi )+\log ({\hat{\sigma }}^2)+1]$$ と書ける。

モデル関数のパラメータ$\theta$の個数を$m$とすれば、ゆらぎの分散 ${\sigma }^2$ の分を加え、モデル全体のパラメータは $m+1$ ある勘定になるから、 AICは $$\text{AIC}=N[\log (2\pi )+\log ({\hat{\sigma }}^2)+1]+2\times (m+1)$$ となる。

以下は、自動車の車速 $x$ と制動距離 $y=h(x)$ の関係を３通りの多項式 $$\begin{array}{r}{h}_1(x)=a_{1}x+a_0\\ h_2(x)=a_2{x}^2+a_{1}x+a_0\\ h_3(x)=a_3{x}^3+a_2{x}^2+a_{1}x+a_0\end{array}$$ でモデル化し、（古い時代の）実測データからそれぞれのAICを計算するコードの例である。

実際にコードを実行すると、二次関数の場合にAICが最も小さくなることがわかる。 自動車の運動エネルギーは速度の２乗に比例する一方で、制動力（摩擦力）が一定とすると、制動で失われるエネルギーは（力）$\times$ (距離)に等しいから、 制動距離は速度の二乗に比例するはずである。 その点で、AICによるモデル比較の結果は妥当と解釈できる。

# coding: utf-8

import scipy
import numpy as np
import matplotlib.pyplot as plt

data = np.array([
    [6.4,0.6],[6.4,3.0],[11.3,1.2],[11.3,6.7],[12.9,4.9],[14.5,3.0],[16.1,5.5],[16.1,7.9],
    [16.1,10.4],[17.7,5.2],[17.7,8.5],[19.3,4.3],[19.3,6.1],[19.3,7.3],[19.3,8.5],
    [20.9,7.9],[20.9,10.4],[20.9,10.4],[20.9,14.0],[22.5,7.9],[22.5,11.0],[22.5,18.3],
    [22.5,24.4],[24.2,6.1],[24.2,7.9],[24.2,16.5],[25.8,9.8],[25.8,12.2],[27.4,9.8],
    [27.4,12.2],[27.4,15.2],[29.0,12.8],[29.0,17.1],[29.0,23.2],[29.0,25.6],[30.6,11.0],
    [30.6,14.0],[30.6,20.7],[32.2,9.8],[32.2,14.6],[32.2,15.8],[32.2,17.1],[32.2,19.5],
    [35.4,20.1],[37.0,16.5],[38.6,21.3],[38.6,28.0],[38.6,28.3],[38.6,36.6],[40.2,25.9]])

n=data.shape[0]

veloc = data[:,0]
dist = data[:,1]

param1 = np.polyfit(veloc, dist, 1)
param2 = np.polyfit(veloc, dist, 2)
param3 = np.polyfit(veloc, dist, 3)

func1 = np.poly1d(param1)
func2 = np.poly1d(param2)
func3 = np.poly1d(param3)

rss_1 = np.sum( np.square(func1(veloc)-dist) ) / n
rss_2 = np.sum( np.square(func2(veloc)-dist) ) / n
rss_3 = np.sum( np.square(func3(veloc)-dist) ) / n

aic_1 = n * (np.log(2*np.pi) + np.log(rss_1) + 1) + 2*(2+1)
aic_2 = n * (np.log(2*np.pi) + np.log(rss_2) + 1) + 2*(3+1)
aic_3 = n * (np.log(2*np.pi) + np.log(rss_3) + 1) + 2*(4+1)

print('AIC(LINEAR FUNCTION):',aic_1)
print('AIC(QUADRATIC FUNCTION):',aic_2)
print('AIC(CUBIC FUNCTION):',aic_3)

plt.plot(veloc,dist, '.', color=(1.0,0,0.0), linewidth=1.0)
plt.xlabel('VELOC[km/h]')
plt.xlim(0,50)
plt.ylabel('BRAKING DISTANCE[m]')
plt.ylim(0,40)
plt.grid(True)
plt.show()

　練習：燃費と車重の関係

複数の自家用車についての燃費（マイル/ガロン）と車重（単位：1000ポンド）のデータをこちらに 用意した。

AICを用いて、両者の関係を表すモデルについて考察してみなさい。