ガウス分布の母平均、母分散の最尤推定-Maximum likelihood estimator

「これならわかる最適化数学」p135~138
最尤推定について書いた。
平均μ、分散σ^2の正規分布の確率密度は
$p(x) = \frac{1}{\sqrt{2 \pi \sigma^2}}\mathrm{e}^{\frac{-(x-\mu)^2}{2\sigma^2}}$
となる。

ここで、正規分布に従った
N個のデータ $x_1,x_2,\dots,x_N$ を観測したとき、
このデータが従う正規分布の平均μ、分散 $\sigma^2$ を推定する。

ここで、 $p(x_1),p(x_2),\dots,p(x_N)$ を観測データの尤度(likelihood)と呼ぶ。これを最大化する操作を行って母平均、母分散を推定する。
これが最尤推定であり、得られる値を最尤推定量と呼ぶ。
（続く）

#-*- coding:utf-8 -*-
#ガウス分布の最尤推定
#Maximum likelihood esstimator
from scipy import stats
from numpy import *
#ガウス分布に従う乱数を生成する。
#これが、最尤推定量がこの値に近づけることを目的としている。

#データ数を設定する
N=100000
P=[]

#平均値を推定する。
def f(x):
	a=sum(x)/len(x)
	return a

#分散値を推定する
def g(x):
	a=f(x)
	D=[]
	for i in range(len(x)):
		D.append((x[i]-a)*(x[i]-a))
 	return sum(D)/len(x)

if __name__=="__main__":
	#母集団の平均値を入力する
	print "母集団の平均値を入力してください。"
	u=input()

	print "母集団の分散を入力してください。"
	b=sqrt(input())

	#正規分布に従うデータx_1~x_Nまで生成
	for i in range(N):
		P.append(stats.norm.rvs(u,b))

	print "母平均,母分散: ",u," ",b*b
	print "最尤推定量の平均、分散: ",f(P)," ",g(P)

実行結果

母集団の平均値を入力してください。
0
母集団の分散を入力してください。
1
母平均,母分散:  0   1.0
最尤推定量の平均、分散:  -0.000390309656977   0.997563352688

あとで、PRMLの議論も盛り込む予定。