統計モデルの比較

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

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

g(y)が真の分布、f(x)がモデルの分布。

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

そのための自然な発想として、対数尤度の平均 (1)L(θ;X)=1N=1Nlogf(x;θ) を最大化するような θ を選べば良さそうである。すなわち (2)θiL(θ^;X)=0 となるようなパラメータの組 θ^ を選べばよさそうだ。

ところが、データは、サンプリングの都度異なるので、最尤推定したパラメータも、確率的に変動する。 fという分布(モデル)を仮定した下での「真の」パラメータは、 (3)L(θ;Y)=EY[logf(y;θ)]=g(y)logf(y;θ)dy を最大化するような値であって、 (4)θiL(θ;Y)=0 となるようなパラメータの組を θ と書くことにする。 そして、N の極限では、両者は一致する(そのような場合を考える)。

では、L(θ^;X) ができるだけ大きくなるようなモデル f(x) が最善のモデルと言えるのだろうか? 統計理論の教えるところによれば、分布関数 g(y)f(x) の隔たり具合(「距離」)は、 カルバック・ライブラー情報量(Kullback-Leibler information) KL(g;f)=g(x)log(g(x)f(x))dx=g(x)logg(x)dxg(x)logf(x)dx を客観的な指標として比較すべきである(KL0で、両者が一致する場合のみ KL=0となる)。 そうすると、最大化しなければならないのは、上式の第2項に対応する (5)L(θ^;Y)=g(θ^;y)logf(θ^;y)dy=EY[logf(y;θ^)] ということになる。

一方、この L(θ^;Y) と、サンプル平均として求めた尤度 L(θ^;X) との間には、 N でも0にならない偏差(バイアス)が存在することが知られている。

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

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

平均対数尤度の微分を、θ の周りで展開すると、以下のようになる。 (6)θiL(θ^;X)=θiL(θ;X)+jk2L(θ;X)θjθi(θ^jθj)+ ここで、定義から上式の左辺は0である。

以下で用いるため、行列 A の各成分を Aij=2logf(x;θ)θjθi とおいて、 2L(θ;X)θjθi=EX[A] と表しておく。また、ベクトル b の成分 bi=θiL(θ;X) とおき、パラメータもベクトル表記 θ=(θ1,θ2,,θk) すると (6)式は (7)EX[A](θ^θ)=b と書ける。

ここで、 (8)bi=1N=1Nlogf(x;θ)θi は、N個の独立な確率変数の和の形をしているので、中心極限定理から正規分布に漸近することが分かる。 その平均 EX[bi] は0で、分散・共分散行列の i,j 成分は (9)EX[bibj]=1NEX[logf(x;θ)θilogf(x;θ)θj]=1NEX[B] となる。 ここで、行列Bの成分を Bij=logf(x;θ)θilogf(x;θ)θj と置いた。

(7)式から θ^θ=EX[A]1b であるから、θ^θ の分散・共分散行列は (10)EX[(θ^θ)(θ^θ)]=EX[A]1EX[bb](EX[A]1)=1NEX[A]1EX[B]EX[A]1 となる。

Nが十分大きければ、Xについての平均はYについてのそれに代えられるので(大数の法則) EX[(θ^θ)(θ^θ)]1NEY[A]1EY[B]EY[A]1

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

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

ここで、(11)式を (12)EX[L(θ^(X);Y)L(θ;Y)]+EX[L(θ;Y)L(θ;X)]+EX[L(θ;X)L(θ^(X);X)] と書き直し、1行目を C1、2行目を C2, 3行目を C3 とする。 言うまでもなく、C=C1+C2+C3 である。

C1の評価

対数尤度をθの周りでテイラー展開すると (13)EY[logf(y;θ^(X))]=EY[logf(y;θ)]+i=1kEY[logf(y;θ)]θi(θ(X)iθi)+12ij2EY[logf(y;θ)]θiθj(θ(X)iθi)(θ(X)jθj)+ を得る。 ここで、θの定義より、右辺2項目(一階微分の項)は0となる。 上式を、対称行列 (14)J=EY[2logf(y;θ)θiθj] を使って整理すると、 (15)L(θ^(X);Y)L(θ;Y)12(θ^θ)J(θ^θ) となる。

trace()は行列の対角要素の和(跡,トレース)。

xをベクトル、Aを対称行列とすると、次の関係が成り立つ: xAx=trace(Axx)

右辺の平均を求めるために、まず、以下の計算をおこなう。 (16)EX[(θ^θ)J(θ^θ)]=EX[trace(J(θ^θ)(θ^θ))]=trace(JEX[(θ^θ)(θ^θ)])1Ntrace(JEY[A]1EY[B]EY[A]1)

ここで、EY[A]=J、すなわち EY[A]1=J1 に他ならないこと、そして、 (17)I=EY[B]=EY[logf(y;θ)θilogf(y;θ)θj] と改めて置き直すことで、(16)式は (18)=1Ntrace(IJ1) となる。

さらに、IJはFisher情報行列の異なる表現に他ならない(I=J)から、結局 (19)=1Ntrace(E)=kN を得る(Eは単位行列)。

以上をまとめると、N での C1 の見積もりとして (20)C1=EX[L(θ^(X);Y)L(θ;Y)]k2N が得られた。

C2の評価

すでに、これまでの計算でも仮定されていたとおり、XYからのサンプルであるから、N では (21)C2=EX[L(θ;Y)L(θ;X)]0 となる。

C3の評価

C1と同様に評価にすると、 (22)C3=EX[L(θ;X)L(θ^(X);X)]k2N を示すことができる(確認せよ)。

以上から、全体的なバイアスの量 (23)C=C1+C2+C3=kN が得られた。

赤池情報量基準(AIC)

前節までの結果から、 十分サンプル数が多く、各サンプルが独立であるような場合、サンプルから推定したパラメータから計算した対数尤度にはバイアスが含まれており、 その母平均は (24)L(θ^(X);Y)=L(θ^(X);X)kN=1N{=1Nlogf(x;θ^)k} と見積もられる。ここで、kはパラメータの個数である。

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

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

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

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

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

パラメータの最尤推定はfit()関数を使って行うことができる。 変量Xは (x - loc)/scale のように、パラメータlocscaleによって変換されるが、 フィッティングの際、それらをある値に固定したい場合は、それぞれ、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)

icon-pc 練習:文長のモデルの比較

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


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

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

データ点の数をNとすると、モデルからの残差の二乗和 iN(yih(xi;θ))2 を最小化するようなパラメータをθ^を与えると、対数尤度を最大化できる(最小二乗法)。 そのときの対数尤度は L(θ^,σ^;{x,y})=N2log(2πσ^2)12σ^2i=1N(yih(xi;θ^))2 となる。

分散の最尤推定値は、不偏分散とはならない点に注意。 実際に L(σ2)=0 を計算して確かめてみよ。

ここで、分散の最尤推定値は σ^2=1Ni=1N(yih(xi;θ))2 になるから、最大対数尤度は、 L(θ^,σ^)=12[log(2π)+log(σ^2)+1] と書ける。

モデル関数のパラメータθの個数をmとすれば、ゆらぎの分散 σ2 の分を加え、モデル全体のパラメータは m+1 ある勘定になるから、 AICは AIC=N[log(2π)+log(σ^2)+1]+2×(m+1) となる。

用いたデータは、Rのサンプルデータ(cars)に基づいている。

以下は、自動車の車速 x と制動距離 y=h(x) の関係を3通りの多項式 h1(x)=a1x+a0h2(x)=a2x2+a1x+a0h3(x)=a3x3+a2x2+a1x+a0 でモデル化し、(古い時代の)実測データからそれぞれのAICを計算するコードの例である。

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

コードの補足:

numpy.polyfit()関数は、データの組から最小二乗法でパラメータの組(上記のai)を計算する。

numpy.poly1d()は、与えられたパラメータの組を元に、多項式関数を生成する。

# 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()

コードを少し変更して、データ点に加え、二次関数でのフィッティング結果もプロットした例。

icon-pc 練習:燃費と車重の関係

データファイル:mtcars-mpg-wt.txt
(Rのサンプルデータ mtcars から作成)

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

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