[統計学]最尤法が苦手?最尤推定法をPythonでやってみる。
深層学習にしても、ベイズ統計にしても、最尤法の考え方が利用されます。
最尤法はどこかわすれがち。
統計学では、学び始めのころは、推定や検定にウエイトが置いて、学習します。
最尤法は独学で、という人も多いでしょう。
わたしもそうでした。
最尤法って、でてくる数式は、二項分布や正規分布など、学び始めに知ったよく使われる確率分布が用いられます。
でも、最尤法はロジックが異なる。
最尤推定法では、$${\hat{\theta}=\text{arg}\medspace\underset{\theta}{\max}\medspace L(\theta)}$$なる、いわば逆関数を求めます。尤度関数を微分などの演算して、最尤推定値$${\hat{\theta}}$$を探します。
$$
\begin{align*}
\hat{\theta}&=\text{arg}\medspace\underset{\theta}{\max}\medspace L(\theta)
=\text{arg}\medspace\underset{\theta}{\max} \displaystyle\prod_{i=1}^n f(x_i:\theta) \\
&=\text{arg}\medspace\underset{\theta}{\max} \left(f(x_1:\theta) \times f(x_2:\theta) \times \cdots \times f(x_n:\theta)\right)
\end{align*}
$$
から想像できるように、ほとんどの場合で、尤度関数の形は複雑になります。最尤推定値は$${\hat{\theta}}$$を具体的に数値計算しないと算出できません。
そこで、
この本、いいです。
最尤法に限定して深く説明した内容ですので、類書はありません。
けっして★★☆☆☆なんて評価の書籍ではありません!
本での勉強では、最尤推定値を具体的な数値を求めるケースはありません。最尤推定が、演習問題にも取り上げられるケースも稀です。
この本では、最尤法を具体的に数値計算して、算出しています。
Excelを使っていますけど、Excel特有の高度な手法ではありません。
Google Spreadsheetでも可能ですし、PythonやRでも可能です。
※ 下でPythonコードを掲載しています。
扱われている最尤法は次の通り。(目次はコチラ)
正規分布の最尤推定値の導出、数値での算出
二項分布分布の最尤推定値の導出、数値での算出
ポアソン分布の最尤推定値の導出、数値での算出
正規分布の最尤推定値の導出、数値での算出
ガンマ分布の最尤推定値の導出、数値での算出
一様分布の最尤推定値の導出、数値での算出
多変量正規分布(2次元)の最尤推定値の導出、数値での算出
回帰(線形/ロジスティクス/ポアソン)の最尤法
そのほかの推定法として
最小二乗法
モーメント法
このように、統計学の範囲内で、最尤推定のトピックは網羅されています。回帰(線形/ロジスティクス/ポアソン)の最尤法を、勉強したことのある人は少ないのでは?
最尤法が苦手かも?と学習がストップしてしまったら、手にとって見てください。
ユーズドで安価で手に入ります。
最尤推定法をPython化してみた。
母集団の総数$${N}$$を、超幾何分布$${H(H,15,26)}$$と二項分布$${Bin(26,15/N)}$$を用いて、最尤推定しています。
最尤推定法での典型問題です。
本ではExcelで計算しています。マクロ・VBA不要で出来ます。
Excelでは、36から86まで尤度関数の算出を、計算式のセルを選択してダブルクリックし、フィルダウンさせて計算させます。
# Google Colaboratoryで実行するときは、あらかじめ
# !pip install japanize-matplotlib
# をしてください。
#%matplotlib inline
import numpy as np
import japanize_matplotlib
import matplotlib.pyplot as plt
from scipy.stats import hypergeom
x = 7
N_values = np.arange(36, 86)
prob_values = []
for N in N_values:
prob = hypergeom.pmf(x, N, 15, 26)
prob_values.append(prob)
japanize_matplotlib.japanize()
plt.plot(N_values, prob_values, 'o-')
plt.xlabel('N')
plt.ylabel('尤度関数')
plt.title('超幾何分布 H(N, 15, 26) for x = 7')
plt.grid(True)
max_index = np.argmax(prob_values)
max_N = N_values[max_index]
max_prob = prob_values[max_index]
plt.plot(max_N, max_prob, 'ro')
plt.text(max_N, max_prob, f'({max_N}, {max_prob:.4f})', ha='left', va='bottom')
plt.show()
print(f"最大の尤度: {max_prob:.4f}")
print(f"それに対応するNの値: {max_N}")
# Google Colaboratoryで実行するときは、あらかじめ
# !pip install japanize-matplotlib
# をしてください。
#%matplotlib inline
import numpy as np
import japanize_matplotlib
import matplotlib.pyplot as plt
from scipy.stats import binom
x = 7
N_values = np.arange(36, 86)
prob_values = []
for N in N_values:
p = 15 / N
prob = binom.pmf(x, 26, p)
prob_values.append(prob)
japanize_matplotlib.japanize()
plt.plot(N_values, prob_values, 'o-')
plt.xlabel('N')
plt.ylabel('尤度関数')
plt.title('二項分布 Bin(26, 15/N) for x = 7')
plt.grid(True)
max_index = np.argmax(prob_values)
max_N = N_values[max_index]
max_prob = prob_values[max_index]
plt.plot(max_N, max_prob, 'ro')
plt.text(max_N, max_prob, f'({max_N}, {max_prob:.4f})', ha='left', va='bottom')
plt.show()
print(f"最大の尤度: {max_prob:.4f}")
print(f"それに対応するNの値: {max_N}")
この記事が気に入ったらサポートをしてみませんか？