#標準正規分布から標本を発生させ、標本から分散の最尤推定量を求める

qiita.rb

import numpy as np
import matplotlib.pyplot as plt
import math

size数分の標本をランダムに抽出し、分散の最尤推定量を求める

qiita.rb

size = 10**7
def max_lilelifood_estimation (size):
#標準正規分布から標本を発生
  a = np.random.randn(size)
  expectedvalue = sum(a)/size
#最尤推定料の方程式
  dispersion = ( sum( [(a[i] - expectedvalue)**2 for i in range (size)] ))/ size
  return dispersion
print(max_lilelifood_estimation (size))

#最尤推定量の分布
$f_N(x) = \frac{1}{\sqrt{2\pi}\sigma} \exp(-\frac{(x-\mu)^2}{2\sigma^2}) $
ここでは、$ x = \sigma$

qiita.rb

number = 10000
dispersion_distribution = [max_lilelifood_estimation(size) for i in range(number)]
dispersion_distribution.sort()
print(dispersion_distribution)
expectation_dispersion = sum (dispersion_distribution) / number
dispersion_dispersion = ( sum( [(dispersion_distribution[i] - expectation_dispersion)**2 for i in range (number)] ))/ number
f = lambda x: math.exp(-((x-expectation_dispersion)**2)/(2*(dispersion_dispersion)**2))/(math.sqrt(2*math.pi)*dispersion_dispersion)
ff = [ f(sel) for sel in dispersion_distribution ]
plt.plot(dispersion_distribution, ff)

#最尤推定法の漸近性
理論上では、フィッシャー情報行列を用いて、証明される
今回は、理論は省略し、実装することで、漸近性を確かめる。

qiita.rb

T = list()
for i in range(4,8):
  a = [((10**i)*(j+1))/100 for j in range(99)]
  T.extend(a)
logT = [math.log10(t) for t in T]
plt.scatter(logT,[max_lilelifood_estimation(int(t)) for t in T])

横軸が、標準正規分布からランダムに抽出する数の対数
縦軸が分散の最尤推定量である

大体収束することがわかる。
これ以上は、分散の最尤推定量を求めるのに時間を消費するので、割愛

#参考資料
統計機械学習の授業