In [42]:

    

# 例5.1 コーシー分布の対数尤度最大化問題を数値計算で解く

import numpy as np
from scipy.stats import cauchy
from scipy.optimize import minimize_scalar

# f(x;0,1)にしたがうコーシー乱数を生成
# このデータのみからパラメータを最尤推定する
# 正解はtheta=0とすでにわかっているが・・・
N = 500
X = cauchy.rvs(loc=0, scale=1, size=N)

def L(theta):
    """コーシー分布の対数尤度関数
    Pythonでは最小化する関数しかない
    そのためターゲット関数を-()で囲むことで最大化問題を最小化問題に変換
    theta   : 推定したいコーシー分布のパラメータ"""
    theta = np.array([theta] * len(X))
    return - (-np.sum(np.log(1 + (X - theta) ** 2)))

# グラフを描画
xlist = np.linspace(-10, 10, 10000)
ylist = np.array([L(x) for x in xlist])
plot(xlist, ylist)
plt.show()

res = minimize_scalar(L, bounds=(-10, 10), method="Brent")
print res.x
print res.fun


    

    













    




0.0168411499014
708.111660358