情報量基準AICとBICのシミュレーション
モデルのパラメータ数を増やしていくと真のモデルとの誤差が増加することがBrunton and Kutz (2019)に説明されている。オーバーフィッティングのために適用モデルが真のモデルから離れていくことをこのウェブサイトにおけるシミュレーションで示した。適用モデルの選択に用いる指標としてBrunton and Kutz (2019)では情報量基準AICとBICについて説明されている。AICとBICについて、Brunton and Kutz (2019)のシミュレーション例であるCode 4.14の設定でシミュレーションを行ってみた。
Code 4.14では、真のモデルは次式で与えられる２次の多項式である。
データは真のモデルに標準偏差0.1の正規分布に従うノイズの加えられたものである。
データの観測点は、区間
において一様に分布するように設定される。
分析に用いられるモデルは、
次多項式である。
情報量基準AICおよびBICは、モデルの推定が最尤法で行われたとき、次式で与えられる(Konishi and Kitagawa, 2008)。
最尤法を適用する確率モデルは、観測点
における観測値が
であるとき、次式で設定する。
このとき、多項式の最尤推定値は次式
の最小2乗法による解で与えられる。すなわち、
である（cf. 岡本、2009、p. 26）。このときの誤差を
とおくとき、
の最尤推定値は次式で与えられる。
シミュレーションプログラムを２本書いてみた。
（１）サンプルサイズn_pnts=100のデータが与えられたときのAICおよびBICの振る舞いを調べた。
（２）データのサンプルサイズをn_pnts=20, 100, 1000と大きくしたときのAICおよびBICの振る舞いを調べた。
（１）サンプルサイズn_pnts=100のときのAICおよびBIC
サンプルスクリプトをリスト１のように用意した。
リスト１のスクリプトをPython 3.11.8/Anacondaで実行すると、まずサンプルサイズn_pnts=100の条件において、n=1000回のシミュレーションの結果が図示される。各シミュレーションにおいて、０次から5次までの多項式モデルが適用され、AICとBICが算出された。０次から５次までの多項式モデルのAICとBICの各1000回のシミュレーションにおける値のボックスプロットと平均値を図１の左側のグラフに示す。２次多項式から5次多項式モデルのグラフは右側のグラフに示した。
図１
0次多項式から２次多項式までAICおよびBICは急激に減少しているが、2次以上の多項式では変化はわずかな増加が認められる。２次以上の多項式モデルでは、ボックスプロットで示される分布の重なりがある。AICあるいはBICで常に正しいモデルが選択できるようではなさそうである。
情報量基準AICおよびBICで正しい選択がどの程度できるか、その頻度を調べた。1000回のシミュレーションにおいて、AICおよびBICにより正しいモデルとして選択された度数を各次数の多項式モデルごとに数え、ヒストグラムで表した（図２）。
図２
正しいモデルである２次多項式モデルの選択頻度がAICおよびBICとも最も高いが、BICの方が正しいモデルである２次多項式モデルの選択がAICより多い。
次に、サンプルサイズを変えたときのAICとBICによるモデル選択の頻度を調べてみる。
（２）サンプルサイズn_pnts=20, 100, 1000のときのAICおよびBIC
サンプルサイズn_pnts=100の条件の他に、それよりサンプルサイズが小さいn_pnts=20の条件と、多いn_pnts=1000の条件についても調べる。スクリプトはリスト２に示すものである。各サンプルサイズの条件において1000回のシミュレーションを行っている。スクリプトを実行すると図３のグラフが表示される。
図３
グラフは、AIC-サンプルサイズ、およびBIC-サンプルサイズごとに正しいモデルとして選択されたモデルの頻度が描かれているが、条件ごとの区別が分かり易いように折れ線グラフで描かれている。AICおよびBICともに、サンプルサイズが大きくなると正しい２次多項式モデルの選択される頻度が高くなることが示されている。また、AICおよびBICともに、サンプルサイズが20の場合に比べて、100の場合は正しいモデルの選択頻度が高くなっているが、サンプルサイズが100の場合と1000の場合とでは差は認め難い。頻度の値は、端末に出力されている。
(py311) PS D:\XXXXXXXX\aic_bic_simultn\samplefiles>
python .\samplesize.py
sample size = 20
Degree
0
1
2
3
4
5
6 7
8
9 10
Freq.(AIC)
0
0
425
81
67
54
55
58
51
80
129
Freq.(BIC)
0
0
700
86
51
34
26 19
23 22 39
sample size = 100
Degree
0
1
2
3
4
5
6
7
8
9
10
Freq.(AIC)
0
0
681
107
68
57
24
13 25
13 12
Freq.(BIC)
0
0
959
36
2
3
0
0
0
0
0
sample size = 1000
Degree
0
1
2
3 4
5
6
7
8
9 10
Freq.(AIC)
0
0
723
105
59
41
22
19
16
7
8
Freq.(BIC)
0
0
994
6 0
0
0
0
0
0
0
(py311) PS D:\XXXXXXXX\aic_bic_simultn\samplefiles>
リスト１ サンプルサイズn_pnts=100のときのAICとBIC
import math
import numpy as np
import matplotlib.pyplot as plt
M = 5
n = 1000
L = 4
n_pnts = 100
x = np.linspace(0, L, n_pnts)
f = x**2
# True value
AICs = np.empty((n, M+1))
BICs = np.empty((n, M+1))
aic_choice_p = []
bic_choice_p = []
for i_trial in range(n):
# True value + Noise
fn = f + 0.1 *
np.random.normal(size=n_pnts)
for p in range(M+1):
A =
np.empty((n_pnts, p+1))
for
j in range(p+1):
A[:,j] = x ** j
b =
np.linalg.pinv(A) @ fn
f_pred = A @ b
sse
= np.sum((fn - f_pred)**2)
sgm
= (sse / n_pnts) ** 0.5
max_log_L = -0.5 * n_pnts * math.log(2*math.pi) \
- n_pnts * math.log(sgm) \
-sse /(2 * (sgm**2))
AICs[i_trial][p] = -2 * max_log_L + 2 * (p + 1 + 1)
BICs[i_trial][p] = -2 * max_log_L + (p +
1 + 1) * math.log(n_pnts)
c_aic =
np.argmin(AICs[i_trial])
aic_choice_p.append(c_aic)
c_bic =
np.argmin(BICs[i_trial])
bic_choice_p.append(c_bic)
plt.figure(figsize=(12,5))
plt.subplot(1,2,1)
Mean_AICs = np.mean(AICs, axis=0)
plt.plot(range(M+1), Mean_AICs, '^-',
c='r', lw=1, label='AIC(mean)')
Mean_BICs = np.mean(BICs, axis=0)
plt.plot(range(M+1), Mean_BICs, 'v-',
c='g', lw=1, label='BIC(mean)')
plt.boxplot(AICs,
positions=np.arange(0,M+1),boxprops={'color':'r'},
flierprops={'markerfacecolor':'r'})
plt.boxplot(BICs,
positions=np.arange(0,M+1),boxprops={'color':'g'},
flierprops={'markerfacecolor':'g'})
plt.xticks(np.arange(0, M+1))
plt.xlabel('Polynomial Degree',fontsize=14)
plt.ylabel('Information Criterion',
fontsize=14)
plt.title('Degree >= 0',
fontsize=16)
plt.legend(fontsize=14)
plt.subplot(1,2,2)
plt.plot(np.arange(M+1)[2:],
Mean_AICs[2:], '^-', c='r', lw=1,
label='AIC(mean)')
plt.plot(np.arange(M+1)[2:],
Mean_BICs[2:], 'v-', c = 'g', lw=1,
label='BIC(mean)')
plt.xticks(np.arange(2, M+1))
plt.xlabel('Polynomial Degree',
fontsize=14)
plt.ylabel('Information Criterion',
fontsize=14)
plt.boxplot(AICs[:,2:],
positions=np.arange(2, M+1), boxprops={'color':'r'},
flierprops={'markerfacecolor':'r'})
plt.boxplot(BICs[:,2:],
positions=np.arange(2, M+1), boxprops={'color':'g'},
flierprops={'markerfacecolor':'g'})
plt.title('Degree >= 2',
fontsize=16)
plt.legend()
plt.tight_layout()
plt.show()
plt.title('Model Selection',
fontsize=16)
plt.hist([aic_choice_p, bic_choice_p],
bins=np.arange(-0.5, M+1, 1.0),
color=['r', 'g'], label=['AIC', 'BIC'])
plt.ylabel('Frequency', fontsize=14)
plt.xlabel('Degree', fontsize=14)
plt.legend(fontsize=14)
plt.tight_layout()
plt.show()
import math
import numpy as np
import matplotlib.pyplot as plt
M = 10
n = 1000
L = 4
choice_aic = np.zeros((3, M+1),
dtype=int)
choice_bic = np.zeros((3, M+1),
dtype=int)
n_pnts_L = [20, 100, 1000]
for i, n_pnts in enumerate(n_pnts_L):
x = np.linspace(0, L,
n_pnts)
f = x**2
# True value
AICs = np.empty((n, M+1))
BICs = np.empty((n, M+1))
for i_trial in range(n):
# True value + Noise
fn =
f + 0.1 * np.random.normal(size=n_pnts)
for
p in range(M+1):
A = np.empty((n_pnts, p+1))
for j in range(p+1):
A[:,j] = x ** j
b = np.linalg.pinv(A) @ fn
f_pred = A @ b
sse = np.sum((fn - f_pred)**2)
sgm = (sse / n_pnts) ** 0.5
max_log_L = -0.5 * n_pnts * math.log(2*math.pi) \
- n_pnts * math.log(sgm) \
-sse /(2 * (sgm**2))
AICs[i_trial][p] = -2
* max_log_L + 2 * (p + 1 + 1)
BICs[i_trial][p] = -2 * max_log_L + (p + 1 + 1) * math.log(n_pnts)
c_aic = np.argmin(AICs[i_trial])
choice_aic[i, c_aic] += 1
c_bic = np.argmin(BICs[i_trial])
choice_bic[i, c_bic] += 1
for i in range(3):
print()
print(f'sample size =
{n_pnts_L[i]}')
print(f'{"Degree":
<15s}', end='')
for j in range(M+1):
print(f'{j:>10d}', end='')
print()
print(f'{"Freq.(AIC)":<15s}', end='')
for j in range(M+1):
print(f'{choice_aic[i,j]:>10d}', end='')
print()
print(f'{"Freq.(BIC)":<15s}', end='')
for j in range(M+1):
print(f'{choice_bic[i,j]:>10d}', end='')
print()
lines = ['--', ':', '-']
for i in range(3):
plt.plot(range(M+1),
choice_aic[i], ls=lines[i], c='r',
label=f'AIC-{n_pnts_L[i]}')
plt.plot(range(M+1),
choice_bic[i], ls=lines[i], c='g',
label=f'BIC-{n_pnts_L[i]}')
plt.title('Model Selection',
fontsize=16)
plt.ylabel('Frequency', fontsize=14)
plt.xlabel('Degree', fontsize=14)
plt.legend(fontsize=16)
plt.tight_layout()
plt.show()
Brunton, S. L., and Kutz, J. N. (2019). Data-driven science and engineering: Machine learning, dynamic systems, and control. Cambridge University Press.
Konishi, S., and Kitagawa, G. (2008). Information criteria and statistical modeling. Springer.
岡本安晴（2009）統計学を学ぶための数学入門［下］．培風館
岡本安晴（2019）いまさら聞けないPythonでデータ分析．丸善出版