こんにちは！データサイエンティストの青木和也（https://twitter.com/kaizen_oni）です！
今回の記事では、統計学の青本「自然科学の統計学」の第7章-演習問題3「最尤推定量」の導出について、Pythonで解いていきたいと思います。
正規分布のような扱いやすい分布でない場合でも、今回の方法を使用すれば近似的な最尤推定量を求めることができます！
扱いづらい分布における母数を推定する際にぜひ利用してみてください！
問題文
＜最尤推定量＞
したのデータに分布(7.8)を当てはめた時の$\theta$の最尤推定量を、本文に述べた方法で求めよ。
中央値を最初の近似解$\theta^{(0)}$をして、逐次近似解$\theta^{(k)}$と係数$W_1^{(k)}$の収束の様子を観察せよ。
$$8.2~~~7.5~~~8.7~~~8.4~~~9.6$$
※(7.8)
$$f_t(x) = \frac2\pi\cdot\frac1{(1+x^2)^2}$$
東京大学教養学部統計学教室『自然科学の統計学』(東京大学出版社/2001) 第7章 演習問題 P228~229
本書に記載の最尤推定量を求める方法
最尤推定量$\hat{\theta}_{ML}$の満たすべき方程式（尤度方程式）は
$$\frac\partial{\partial \theta}\{\sum_i \log f_{\theta}(x_i)\} = \sum_i \frac{x_i – \theta}{1 + (x_i – \theta)^2} = 0~~~(7.16)$$
上式の解$\theta$は直接得られないが、適当な近似解から出発して、逐次近似によって求めることができる。今、
$$W_i = \frac1{1 + (x_i – \theta)^2}$$
とおくと、(7.16)は形式的に$\sum W_i(x_i – \theta) = 0$、したがって、
$$\theta = \frac{\sum W_i x_i}{\sum W_i}$$
と表される。そこで、最初の近似解$\theta^{(0)}$を用いて、
$$W_i^{(0)} = \frac1{1 + (x_i – \theta^{(0)})^2}, ~~~ i = 1,2,\cdots, n$$
とおき、
$$\theta^{(1)} = \frac{\sum W_i^{(0)}x_i}{\sum W_i^{(0)}$$
を新しい近似解とし、収束するまで次々と新しい近似解を計算すれば良い。
東京大学教養学部統計学教室『自然科学の統計学』(東京大学出版社/2001) 第7章 演習問題 P209~210
実装コード
本問題を解くための実装コードは以下のとおりです。
import numpy as np
# データ
arr = np.array([8.2, 7.5, 8.7, 8.4, 9.6])
# 最初の近似解
med = np.median(arr)
theta0 = med
# 近似解の分母
def W(arr, theta):
sum_ = (1/(1+ (arr - theta)**2)).sum()
return sum_
# 近似解の分子
def Wx(arr, theta):
sum_ = (arr/(1+ (arr - theta)**2)).sum()
return sum_
theta = theta0
print(f'theta0: {theta}')
while True:
old_theta = theta
new_theta = Wx(arr, theta) / W(arr, theta)
print(f'W: {W(arr, theta)}')
print(f'theta: {new_theta}')
theta = new_theta
# 前回との差分を比較、差分0.00001未満で収束と判定
diff = np.abs(new_theta - old_theta)
if diff < 0.00001:
print('Finish')
break解説
本実装コードでは、以下の3ステップを実行することにより近似的な最尤推定量を導出しています。
- データから最初の近似解$\theta^{(0)}$を求める
- 近似解の分母・分子を導出する関数を定義する
- 近似解を随時更新し、近似解の差分が0.00001未満になったら収束と判定し、処理を終了する
順を追って解説していきます。
データから最初の近似解$\theta^{(0)}$を求める
本問題においては、データの中央値を近似解の初期値と設定します。
# データ
arr = np.array([8.2, 7.5, 8.7, 8.4, 9.6])
# 最初の近似解
med = np.median(arr)
theta0 = med近似解の分母・分子を導出する関数を定義する
逐次近似解$\theta^{(k)}$を求める式
$$\theta^{(k)} = \frac{\sum W_i x_i}{\sum W_i}$$
の分母と分子を求める関数を定義します。
# 近似解の分母
def W(arr, theta):
sum_ = (1/(1+ (arr - theta)**2)).sum()
return sum_
# 近似解の分子
def Wx(arr, theta):
sum_ = (arr/(1+ (arr - theta)**2)).sum()
return sum_近似解を随時更新し、近似解の差分が0.00001未満になったら収束と判定し、処理を終了する
近似解の初期値$\theta^{(0)}$を使用して逐次近似解の導出をスタートし、前回求めた近似解と今回求めた近似解の差分が0.00001未満になった段階で処理を終了し、その時点での近似解を最尤推定量とします。
theta = theta0
print(f'theta0: {theta}')
while True:
old_theta = theta
new_theta = Wx(arr, theta) / W(arr, theta)
print(f'W: {W(arr, theta)}')
print(f'theta: {new_theta}')
theta = new_theta
# 前回との差分を比較、差分0.00001未満で収束と判定
diff = np.abs(new_theta - old_theta)
if diff < 0.00001:
print('Finish')
break
まとめ
今回の記事では、統計学の青本「自然科学の統計学」の第7章-演習問題3「最尤推定量」の導出について、Pythonで解く実装コードを紹介し、その手順を簡単に説明しました。
このような近似的な手法を使えば、正規分布のような扱いやすい分布が仮定できない場合についても、対数微分をすることができれば、近似的な最尤推定量を求めることができます。
皆さんもぜひ、分布の母数を求める際にPythonコードを利用した近似的手法を用いてみてください！
コメント