最尤法
岡本安晴
2019.6、2025.3
データ分析において設定された確率モデルのパラメータ値の推定には、いくつかの方法がある。簡単なものとしては、データから算出される平均値や分散（これらはモーメントと呼ばれる）と確率モデルに基づいて算出される平均値や分散を等しいと置いてパラメータ値を推定する方法で、モーメント法と呼ばれている。ここで説明する最尤法では、データの確率が最大となるパラメータの値として求められ、最尤推定値と呼ばれている。プログラムによるシミュレーションを行いながら具体的に説明していく。
本文中のスクリプトファイルは、圧縮ファイルmlfiles.zipにまとめた。
いま、表と裏の模様の凹凸がかなり異なるコインの表の出る確率
を調べるために、コインを１0回投げたところ表が７回出たとする。この結果に対して、二項分布のモデルを適用して確率を推定する。
二項分布は次式で与えられる。
ここで、
は試行数、は確率、
は表の出た回数を表す。
試行において
回となる確率は
と表されるが、これはの値は不明なのでの関数の形である。データが与えられたとき、そのデータが得られる確率を、確率分布のパラメータの関数として見たものは尤度関数と呼ばれている。上の尤度関数のグラフを描くスクリプトを以下のように用意した（リスト１）。
import matplotlib.pyplot
as plt
import numpy
as np
import scipy.special
def L_func(p):
return scipy.special.binom(10,
7) * (p ** 7) * ((1 - p) ** (10 -7))
p_values = np.arange(0.0,
1.005, 0.01)
L_values = L_func(p_values)
plt.plot(p_values, L_values)
plt.show()
実行すると図１に示すグラフが表示される。
図１のグラフは、
のときに最大値をとる。すなわち、であると考えると、データの得られる確率が最大となるので、がデータに対して最も尤もらしい値と考えられる。尤度関数を最大にするパラメータの値を最尤推定値と呼ぶので、が最尤推定値となる。
試行数、事象の生起確率、生起回数の二項分布の場合の尤度関数（likelihood function）は、の関数であることを明示して
と書くが、に関して定数項は最大値を考えるとき関係ないので省いて、
と書くことが多い。
最尤推定値は
を最大にするの値であるが、は
のとき最大となる。の推定値であることを示すため上にハット＾を付けて表すので、
である。
最尤推定値の性質を調べるために、確率の真値が
であるとしてシミュレーションを行う。
回の試行を行い事象の生起回数を数えるという実験を10000回行うことは次のスクリプトで表せる（リスト２）。
N = 10
p = 0.7
p_hat = []
for t in range(10000):
k = scipy.stats.binom.rvs(N,
p)
p_hat.append(k
/ N)
リストp_hatに格納されている
のサンプルのヒストグラムは、次のスクリプトで描ける。
plt.hist(p_hat, bins
= [-0.025 + i / 20 for i in
range(21)])
plt.show()
実行すると、図２のグラフが表示される。
真値を中心とした分布になっているが、最尤推定値
が常に真値と一致しているわけではない。
図2に描かれている分布の累積比率を描いてみる。スクリプトは以下のように用意した。
s_p_hat = np.sort(p_hat)
n_len = s_p_hat.size
cum_s_p_hat = []
for i in
range(n_len):
cum_s_p_hat.append((i + 1) / n_len)
plt.plot(s_p_hat, cum_s_p_hat)
plt.show()
実行すると図３のグラフを得る。
より詳しく見るために、最尤推定値の分布を標準化したグラフを描いてみる。スクリプトを以下に示す。
z_p_hat = (s_p_hat
- s_p_hat.mean()) / s_p_hat.std()
min_x = z_p_hat.min()
max_x = z_p_hat.max()
x_range = np.arange(min_x, max_x, 0.01)
cum_norm = scipy.stats.norm.cdf(x_range)
plt.plot(z_p_hat, cum_s_p_hat)
plt.plot(x_range, cum_norm, 'r--')
plt.show()
実行すると図４のグラフが描かれる。赤の曲線は累積標準正規分布関数である。標準化最尤推定値は
間隔の離散量であるので、連続量の累積分布関数である累積標準正規分布関数のグラフとの違いが目立つ。
試行数の値を
と大きくとって最尤推定値
の分布を調べてみる。
スクリプトを以下のように用意した（リスト３）。
N = 1000
p = 0.7
p_hat = []
for t in
range(10000):
k = scipy.stats.binom.rvs(N,
p)
p_hat.append(k
/ N)
plt.hist(p_hat, bins
= [0.0005 + i/1000 for i in
range(600, 800, 2)])
plt.show()
スクリプトを実行すると図５のヒストグラムを得る。最尤推定値はかなり連続量に近い分布になっている。
累積比率分布を描くスクリプトを次のように用意した。
s_p_hat = np.sort(p_hat)
n_len = s_p_hat.size
cum_s_p_hat = []
for i in
range(n_len):
cum_s_p_hat.append((i + 1) / n_len)
plt.plot(s_p_hat, cum_s_p_hat)
plt.show()
このスクリプトを実行すると図６のグラフを得る。
図５のヒストグラムでは、ランダムなサンプリングによる分布の不規則な凸凹の変動が見られるが、累積比率曲線の場合はランダムな変動が相殺されるので不規則な凸凹は見られない。滑らかな階段状である。
最尤推定値を標準化して累積比率関数を描いてみる。スクリプトを以下に示す。
z_p_hat = (s_p_hat
- s_p_hat.mean()) / s_p_hat.std()
min_x = z_p_hat.min()
max_x = z_p_hat.max()
x_range = np.arange(min_x, max_x, 0.01)
cum_norm = scipy.stats.norm.cdf(x_range)
plt.plot(z_p_hat, cum_s_p_hat)
plt.plot(x_range, cum_norm, 'r--')
plt.show()
実行すると図７のグラフを得る。
赤い曲線は累積標準正規分布関数である。かなり滑らかな最尤推定値の累積比率関数と累積標準正規分布関数が極めてよく一致していることが分かる。
一般に最尤推定値の分布は、データ数（二項分布の場合は総試行数）が多くなると真値を平均値とする正規分布に近付き、その正規分布の分散はフィッシャーの情報量の逆数で与えられるという性質がある（岡本、2014、pp. 225-227）。
import matplotlib.pyplot
as plt
import numpy
as np
import scipy.special
def L_func(p):
return scipy.special.binom(10,
7) * (p ** 7) * ((1 - p) ** (10 -7))
p_values = np.arange(0.0,
1.005, 0.01)
L_values = L_func(p_values)
plt.plot(p_values, L_values)
plt.show()
import matplotlib.pyplot
as plt
import numpy
as np
import scipy.stats
N = 10
p = 0.7
p_hat = []
for t in range(10000):
k = scipy.stats.binom.rvs(N,
p)
p_hat.append(k
/ N)
plt.hist(p_hat, bins
= [-0.025 + i / 20 for i in
range(21)])
plt.show()
s_p_hat = np.sort(p_hat)
n_len = s_p_hat.size
cum_s_p_hat = []
for i in
range(n_len):
cum_s_p_hat.append((i + 1) / n_len)
plt.plot(s_p_hat, cum_s_p_hat)
plt.show()
z_p_hat = (s_p_hat
- s_p_hat.mean()) / s_p_hat.std()
min_x = z_p_hat.min()
max_x = z_p_hat.max()
x_range = np.arange(min_x, max_x, 0.01)
cum_norm = scipy.stats.norm.cdf(x_range)
plt.plot(z_p_hat, cum_s_p_hat)
plt.plot(x_range, cum_norm, 'r--')
plt.show()
import matplotlib.pyplot
as plt
import numpy
as np
import scipy.stats
N = 1000
p = 0.7
p_hat = []
for t in range(10000):
k = scipy.stats.binom.rvs(N,
p)
p_hat.append(k
/ N)
plt.hist(p_hat, bins
= [0.0005 + i/1000 for i in
range(600, 800, 2)])
plt.show()
s_p_hat = np.sort(p_hat)
n_len = s_p_hat.size
cum_s_p_hat = []
for i in
range(n_len):
cum_s_p_hat.append((i + 1) / n_len)
plt.plot(s_p_hat, cum_s_p_hat)
plt.show()
z_p_hat = (s_p_hat
- s_p_hat.mean()) / s_p_hat.std()
min_x = z_p_hat.min()
max_x = z_p_hat.max()
x_range = np.arange(min_x, max_x, 0.01)
cum_norm = scipy.stats.norm.cdf(x_range)
plt.plot(z_p_hat, cum_s_p_hat)
plt.plot(x_range, cum_norm, 'r--')
plt.show()
参考文献
岡本安晴（2019）．いまさら聞けないPythonでデータ分析．丸善出版