다시 해보자 최우 추정과 확률 분포의 맞추기 ②연속형 확률 분포

소개

이전 내용 : 이산형 확률 분포의 최대 우도 추정과 적용
이번은 연속형 확률 분포의 최대 우도 추정과 적합합니다.

연속형 확률 분포: 정규 분포

정규 분포를 나타내는 확률 함수는 다음과 같습니다.
$$P(y)=\frac{1}{\sqrt{2πσ^2}}e^{(-\frac{(y_i-\mu)^2}{2\sigma^2})}　(\mu ,\sigma = 매개 변수)$$

추정과 분포를 맞추는 테스트 데이터는 대답을 맞출 수 있도록 미리 대상의 확률 분포로부터 작성해 둡니다.

파이썬

"""正規分布からの乱数"""
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씩으로 구분하여 카운트하고 있습니다.

그런데, 이러한 데이터에 대해서, 정규 분포를 맞추어 그 파라미터를 최대 우도 추정해 보겠습니다.

파이썬

""" 対数尤度関数を偏微分してパラメータ推定 """
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를 요구하면 되므로

파이썬

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의 설정대로 되어 있네요.

파라미터 추정이 복잡한 경우

자세한 것은 전회의 ① . 비선형 방정식을 푸는 방법을 사용합니다.

파이썬

"""確率密度関数"""
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

손쉽게 매개 변수를 알고 싶습니다 ...

자세한 것은 전회의 ① . 비선형 방정식을 푸는 방법을 사용합니다.

파이썬

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

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

데이터에 적용하여 그림을 보자.

추정한 파라미터 μ=-0.01455655, σ=0.93795781로서 확률 분포(정규 분포)를 데이터에 적용해 보겠습니다.

파이썬

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)

좋은 느낌이군요.