다시 해보자 최대 우도 추정과 확률 분포의 적용 ①이산형 확률 분포

소개

데이터 분석의 화라고 하면, 관측 데이터에 숨어 있는 관계성을 설명할 수 있는 통계 모델의 구축이군요. 통계 모델의 구축은 확률 분포를 이해할 수 없습니다.
데이터 분석의 관점에서, 특정 변수의 데이터는 확률적인 현상에서 발생합니다. 반대로 관측된 데이터의 패턴은 어떤 확률 분포로 표현할 수 있을까?

최대 우도 추정법

그것에는 확률 분포를 적용해 보는 것이 하나의 손입니다. 그렇다고해도 확률 분포는 파라미터별로 형태가 변화하기 때문에, 맞추기 위해서는 관측 데이터에 근거한 적절한 파라미터치를 결정해 주어야 합니다. 매개 변수 값을 계산하는 주요 방법으로 최대 우도 추정 방법이 있습니다.

최대 우도 추정법의 방법을 대략 쓰면, 얻어진 데이터가 있는 파라미터 θ를 가지는 확률 분포로부터 추출된 N개의 데이터 Y라고 가정해, 값 하나하나의 추출 확률 $p(y_i|θ )$를 N개씩 합친 것 우도 함수 L로 합니다.
$$L(θ|{y_i}) = p(y_1|θ) × p(y_2|θ) × ··· × p(y_N|θ)=\prod_{i=1}^{N}p(y_N |θ)$$
우도 함수 L이 최대가 되는 파라미터 값을 계산합니다. 계산하기 쉬운 형태로 변환(대수화)해 파라미터의 수만큼 편미분해 해석적으로 풀 뿐...이라고 하는 흐름입니다.

실제로 해보자

확률 분포는 데이터에 따라 이산형과 연속형 2가지가 있지만, 이번에는 이산형의 확률 분포의 예로서 포아송 분포의 최대 우도 추정과 확률 분포의 적합을 하고 있다.

이산형 확률 분포 : 포아송 분포 (이번)

연속형 확률 분포 : 정규 분포 (다음 번)

이산형 확률 분포: 포아송 분포

포아송 분포는 카운트 데이터에 주로 사용됩니다. 분포를 나타내는 확률 함수는 다음과 같습니다.
$$
P(y)=\frac{λ^ye^{-λ}}{y!}(λ=y의 평균을 나타내는 파라미터, y=0,1,2,3...)
$$

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

파이썬

"""ポアソン分布からの乱数 """
import numpy as np
import matplotlib.pyplot as plt

np.random.seed(seed=10)
# パラメータλ=2.4
poisson_values = np.random.poisson(lam=2.4, size=1000)
p_y, bin_edges, patches = plt.hist(poisson_values , bins=11, range=[-0.5, 10.5], normed=True)
y_k = 0.5*(bin_edges[1:] + bin_edges[:-1])
y_k #y
p_y #p(y)

얻은 1000개의 의사 데이터의 값과 그 수가 나오고 있네요.
그런데 이러한 데이터에 대해 포아송 분포를 맞추고 그 파라미터를 최대 우도 추정해 보겠습니다.

파이썬

""" 対数尤度関数を偏微分してパラメータ推定 """
import sympy #代数計算するライブラリ（微分や方程式を解く）
# 変数を定義
sympy.var('λ y ')
# 尤度p(パラメータ|y)を定義
f = (λ**y)*sympy.exp(-λ)/sympy.factorial(y)
# 対数化
logf=sympy.log(f)
# fを偏微分して、式を展開
pdff = sympy.expand(sympy.diff(logf, λ))

로그 우도를 파라미터 λ에 대해 편미분한 식 pdff가 $\frac{y_i}{λ}-1$가 되고 있네요. 이것을 데이터분 합한 $\sum_{i=1}^{1000}(\frac{y_i}{λ}-1)=0$가 되는 λ를 구해 주면 되므로

파이썬

def L_sympy(f,var,values):       
    likelihood = 0 #尤度の初期値
    for i in np.arange(len(values)): #データの個数分
        # model output 
        # print(values[i])
        likelihood += f.subs(var,values[i]) #yに値を代入
        # print(likelihood)
    param = sympy.solve(likelihood, λ) #λについての方程式を解く
    # print(param)
    return param

L_sympy(pdff,"y",poisson_values)

パラメータ推定量：[289/125]=2.312

파라미터 λ는 2.312로 추정되었다. 이 의사 데이터의 λ의 설정은 2.4이므로, 좋은 느낌으로 추정할 수 있네요.

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

교과서적인 방법은 위와 같이 편미분한 식(도함수)을 풀면 파라미터를 계산할 수 있는 식을 직접 얻을 수 있습니다. 그러나 단순한 1변수 데이터에 확률 분포의 적용과 달리 일반적인 통계 모델링에서는 링크 함수로 변환한 선형 모델을 맞추기 위해 추정하는 파라미터가 복수가 되어, 식도 보다 복잡하기 때문에 간단하게는 풀 수 없습니다.

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

따라서 비선형 방정식을 푸는 방법을 사용합니다. (이른바 $\frac{y}{λ}-1$와 같은 선형 공식이 아니라 $λ^2+θ^3+...$와 같은 비선형).

파이썬

"""離散型確率分布の関数（確率質量関数）"""
def probability_poisson(y,λ):
    from scipy.special import factorial
    # λ: パラメータ
    # y: データy：y回発生する
    # return: データYの確率P(y|パラメータ) 
    return (λ**y)*np.exp(-λ)/factorial(y)

"""対数尤度関数：離散型"""
def L_func(param,y):       
    likelihood = 0 #尤度の初期値
    for i in np.arange(len(y)):
        # model output 
        p = probability_poisson(y[i], param)
        # likelihood
        likelihood += -np.log(p) #尤度
    return likelihood 

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

params_MLE = optimize.minimize(L_func,x0,args=(poisson_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=1
# AIC（一番小さい値がはてまりの良い確率分布）
print('AIC：',params_MLE.fun*(2)+2*k)

パラメータ： [2.31200054]
AIC： [3607.0081302]

마찬가지로 파라미터 λ는 2.312로 추정되었다.

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

파이썬

"""curve_fitでとりあえずパラメータ計算するのもあり"""
from scipy.optimize import curve_fit
parameters, cov_matrix = curve_fit(f=probability_poisson, xdata=y_k, ydata=p_y) 
print("パラメータ:",parameters, "共分散:",cov_matrix)

パラメータ: [2.31643586] 共分散: [[0.00043928]]

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

추정한 파라미터 λ=2.312로서 확률 분포(포아송 분포)를 데이터에 적용해 보겠습니다.

파이썬

"""当てはまっているかの図示"""
acc_mle = probability_poisson(y_k, params_MLE.x)
plt.hist(poisson_values , bins=11, range=[-0.5, 10.5], normed=True)
plt.plot(y_k,acc_mle)

좋은 느낌이군요.