統計モデルの比較
このページでは、モデルの「良さ」の比較について考えてみる。
パラメーター最尤値と尤度のバイアス
ある確率分布
そのための自然な発想として、対数尤度の平均
ところが、データは、サンプリングの都度異なるので、最尤推定したパラメータも、確率的に変動する。
では、
一方、この
パラメーターのばらつき具合
サンプルから求めたパラメータ
平均対数尤度の微分を、
以下で用いるため、行列
ここで、
(7)式から
対数尤度のバイアス量の見積もり
データから推定したパラメータの組
ここで、(11)式を
の評価
対数尤度を
trace()は行列の対角要素の和(跡,トレース)。
右辺の平均を求めるために、まず、以下の計算をおこなう。
ここで、
さらに、
以上をまとめると、
の評価
すでに、これまでの計算でも仮定されていたとおり、
の評価
以上から、全体的なバイアスの量
赤池情報量基準(AIC)
前節までの結果から、
十分サンプル数が多く、各サンプルが独立であるような場合、サンプルから推定したパラメータから計算した対数尤度にはバイアスが含まれており、
その母平均は
尤度の最大値
上式のカッコ内の量を-2倍した量は、赤池情報量基準(AIC)と呼ばれる:
現実のデータに適用するに際しては、いくつかの改良版の基準が提案されているので、ここから先は専門書などに当たること。
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)
練習:文長のモデルの比較
夏目漱石の「こころ」の文の長さのデータをこちらに用意した。 文の長さが確率的に生成されていると仮定(単純化)した場合、どのようなどのような確率分布で表現するのがよいか、検討してみなさい。
最小二乗法による回帰曲線
データ点の組
データ点の数を
分散の最尤推定値は、不偏分散とはならない点に注意。
実際に
ここで、分散の最尤推定値は
モデル関数のパラメータ
用いたデータは、Rのサンプルデータ(cars)に基づいている。
以下は、自動車の車速
実際にコードを実行すると、二次関数の場合にAICが最も小さくなることがわかる。
自動車の運動エネルギーは速度の2乗に比例する一方で、制動力(摩擦力)が一定とすると、制動で失われるエネルギーは(力)
コードの補足:
numpy.polyfit()
関数は、データの組から最小二乗法でパラメータの組(上記の
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()
コードを少し変更して、データ点に加え、二次関数でのフィッティング結果もプロットした例。
練習:燃費と車重の関係
データファイル:mtcars-mpg-wt.txt
(Rのサンプルデータ mtcars から作成)
複数の自家用車についての燃費(マイル/ガロン)と車重(単位:1000ポンド)のデータをこちらに 用意した。
AICを用いて、両者の関係を表すモデルについて考察してみなさい。