最尤推定値の漸近分布
簡単なシミュレーション例
パラメータの推定値が尤度関数を最大にする値として与えられる最尤推定値には、データ数を増やすと正規分布に近付くという漸近定理がある。すなわち、
である。ここで、
はパラメータの真値、
は最尤推定値、
はデータ数、
はフィッシャー情報行列である。
N試行の2項分布において、成功がｋ回であったとき、成功確率の最尤推定値は次式で与えられる。
また、
である（岡本、2009、p. 158参照）。
いま。コイン投げにおいて、表と裏のデザインの違いにより表の出る確率が0.6であるとする。このコインの表の出る確率を、コインをN=1000投げて最尤推定値で推定するものとする。この1000回投げるコイン投げによる最尤推定値の分布を2000回のシミュレーションによって調べるスクリプトを以下のように作成した。スクリプトファイルは圧縮ファイルfilesMLMaxEst.zipとしてまとめたので、ダウンロード解凍すれば、自由に利用できる。
import numpy
as np
import scipy.stats
as ss
import matplotlib.pyplot
as plt
p0 = 0.6
# 真値
NTrials = 2000 # シミュレーションの試行数
N = 1000
# 2項分布の試行数
p_hat_est = [] # 推定値のリスト
for t in range(NTrials):
if (t % 100) == 0:
print('t = ', t)
R = np.empty((N),
dtype = int)
for i
in range(N):
if ss.uniform.rvs() < p0:
R[i] = 1
else:
R[i] = 0
p_hat
= R.mean()
# 最尤推定値
p_hat_est.append(p_hat)
#
#
推定値の変換
#
n_p_est = (np.array(p_hat_est) - p0) * (N**0.5)
n_p_est_s = sorted(n_p_est) # 累積分布を描くための並べ替え
#
# 累積度数の比率
#
cum_v = (np.arange(0,
len(n_p_est_s)) + 1) / len(n_p_est_s)
plt.plot(n_p_est_s, cum_v)
#
# 累積正規分布の描画
#
xvalues_n = np.arange(min(n_p_est_s), max(n_p_est_s), 0.01)
norm_v = ss.norm.cdf(xvalues_n, scale = (p0 * (1.0 - p0)) ** 0.5)
plt.plot(xvalues_n, norm_v)
plt.show()
上のスクリプトを実行すると図１のグラフが表示される。
図１
細かい階段状の折れ線が最尤推定値の累積度数の比率であり、滑らかな赤の曲線が累積正規分布である。最尤推定値の分布が正規分布でよく近似されていることが分かる
データ数Nを５に減らしてみると、図２のグラフを得る。
図２
最尤推定値の累積度数の比率を表す折れ線は、滑らかな累積正規分布とはかなり異なる。
参考文献
岡本安晴（2009）「データ分析のための統計学入門」おうふう．