はじめに

前回の内容：離散型確率分布の最尤推定と当てはめ
今回は連続型確率分布の最尤推定と当てはめです。

連続型確率分布：正規分布

正規分布を表す確率関数は以下の通りです。
$$P(y)=\frac{1}{\sqrt{2πσ^2}}e^{(-\frac{(y_i-\mu)^2}{2\sigma^2})}　（\mu,\sigma =パラメータ)$$

推定と分布の当てはめをするテストデータは答え合わせができるよう、あらかじめ対象の確率分布から作成しておきます。

Python

"""正規分布からの乱数"""
import numpy as np
import matplotlib.pyplot as plt

np.random.seed(seed=10)
norm_values = np.random.normal(0, 1, size=1000) # 平均0,標準偏差1の正規分布から1000個抽出 
# 描画
plt.hist(norm_values, bins=50, range=[-5, 5], normed=True)
bin_edges

連続変数なので離散型のようにカウントできません。そのため0.2ずつに区切ってカウントしています。

さて、このようなデータに対し、正規分布を当てはめそのパラメータを最尤推定してみます。

Python

""" 対数尤度関数を偏微分してパラメータ推定 """
import sympy
# 変数を定義（v=σ**2としておく）
sympy.var('μ v y')
# 尤度p(パラメータ|x)を定義
fe=(1/sympy.sqrt(2*sympy.pi*v))*sympy.exp(-(y-μ)**2/(2*v))
# 対数化
logf=sympy.log(fe)
# fを偏微分して、式を展開
pdff1 = sympy.expand(sympy.diff(logf, μ)) #μについて偏微分
pdff2 = sympy.expand(sympy.diff(logf, v)) #vについて偏微分

対数尤度をμについて偏微分した式pdff1が$ \frac{y}{v}-\frac{\mu}{v}$ 、vについて偏微分した式pdff2が$\frac{-1}{2v}+\frac{y^2}{2v^2}-\frac{y\mu}{v^2}+\frac{\mu^2}{2v^2}$ですね。$\sum pdff1=\sum pdff2=0$になるμ,vを求めてやればよいので

Python

def L_sympy(fmu,fs,var,values):       
    likelihood_mu = 0 #尤度の初期値
    likelihood_s = 0 #尤度の初期値
    for i in np.arange(len(values)):
        # likelihood
        likelihood_mu += fmu.subs(var,values[i]) #xに値を代入
        likelihood_s += fs.subs(var,values[i]) #xに値を代入
    param = sympy.solve([likelihood_mu,likelihood_s]) #方程式を解く
    return param

parameters = L_sympy(pdff1,pdff2,"y",norm_values)
parameters[0]["σ"]=sympy.sqrt(parameters[0][v])
parameters

[{v: 0.879764754284410, μ: -0.0145566356154705, 'σ': 0.937957757196138}]

大体平均μ=0、σ=1の設定どおりになってますね。

パラメータ推定が複雑なケース

詳細は前回の①。非線形方程式を解く方法を使います。

Python

"""確率密度関数"""
def probability_function(y,param):
    from scipy.special import factorial
    # param[0]s,param[1]mu: パラメータ
    # y: データy
    # return: データyの確率密度P(y|パラメータ) 
    return (1/np.sqrt(2*np.pi*param[1]**2))*np.exp(-0.5*(y-param[0])**2/param[1]**2)

"""対数尤度関数（連続）"""
def L_func_c(param,y):       
    likelihood = 0 #尤度の初期値
    for i in np.arange(len(y)):
        # model output 
        p = probability_function(y[i], param)
        # likelihood 尤度
        likelihood += -np.log(p) #尤度
    return likelihood

"""パラメータ推定"""
"""
準Newton法（BFGS，L-BFGS法：複雑な場合にメモリの節約などが可能）
"""
from scipy import optimize
x0 = [0,0.1] #パラメータの初期値
bound = [(-100, 100),(0, None)]#,(0,None)] #最適パラメータ検索の範囲(min,max)

params_MLE = optimize.minimize(L_func_c,x0,args=(norm_values),method='l-bfgs-b',
                        jac=None, bounds=bound, tol=None, callback=None,
                        options={'disp': None, 'maxls': 20, 'iprint': -1,
                                 'gtol': 1e-05, 'eps': 1e-08, 'maxiter': 15000,
                                 'ftol': 2.220446049250313e-09, 'maxcor': 10,
                                 'maxfun': 15000})
# 最尤推定パラメータ
print('パラメータμ,σ：',params_MLE.x)
# パラメータ数
k=3
# AIC（一番小さい値がはてまりの良い確率分布）
print('AIC：',params_MLE.fun*(2)+2*k)

パラメータμ,σ： [-0.01455655  0.93795781]
AIC： 2715.7763344850605

手っ取り早くパラメータを知りたい・・・

詳細は前回の①。非線形方程式を解く方法を使います。

Python

"""scipy.stats.fitでパラメータ推定するのもあり"""
from scipy.stats import norm
fit_parameter = norm.fit(norm_values)
fit_parameter

# パラメータμ,σ
(-0.014556635615470447, 0.9379577571961389)

データに当てはめて図示してみる

推定したパラメータμ=-0.01455655, σ=0.93795781として、確率分布（正規分布）をデータに当てはめてみます。

Python

acc_mle = probability_function(np.sort(norm_values), params_MLE.x)
plt.hist(norm_values, bins=50, range=[-5, 5], normed=True)
plt.plot(np.sort(norm_values), acc_mle)

良い感じですね。