# MLE estimator

앞선 포스팅에서 언급한 바이지만, 빈도론자는 주어진 데이터를 표상하는 고정된 모수(parameter)가 이미 존재한다고 생각합니다. 그러한 데이터를 발생시킬 가능성(likelihood)이 가장 높은 파라메터 값을 찾아낼 수 있는 방법론중의 하나가 바로 MLE(Maximum Likelihood Estimation; 최대가능도추정)입니다.

수식으로 표현하면 아래와 같이 먼저 likelihood function을 정의하고,

$\mathcal{L}(\theta;x)$

이를 최대화하는 파라메터 $\theta$를 다음과 같이 추정하는 것입니다.

$\hat{\theta}=\arg\max_\theta\mathcal{L}(\theta;x)$

$\beta$를 maximum likelihood estimator를 이용해서 추정 할 수도 있습니다.

모델과 가정은 OLS estimator의 경우와 동일하게 출발합니다.

$y_i=x_i\beta+\epsilon_i$

$y=X\beta+\epsilon$

$\epsilon\sim{N(0,\sigma^2I)}$

• Linear relationship between input and output
• Full rank design matrix
• Multivariate normal distribution of $\epsilon$
• Diagonal covariance matrix of $\epsilon$
• Equal diagonal entries of $\epsilon$

단, 여기서는 $N$개의 샘플이 아래의 IID를 만족한다고 가정합니다. 왜냐하면 해당 가정이 있어야만 likelihood function을 product($\prod$)로 세울 수 있고, 여기에 log를 취해 sum($\Sigma$)으로 바꾸어 계산할 수 있기 때문입니다.

• Independent and Identically Distributed

이제 정규 확률 변수의 선형변환도 정규분포를 따른다는 성질에 따라 종속변수 $y_i$ 역시 평균 $x_i\beta_0$, 분산 $\sigma_0^2$의 조건부 정규분포를 따르게 되며, 해당 pdf는 아래의 식과 같이 표현할 수 있습니다.

$f_Y(y_i|X)=(2\pi\sigma_0^2)^{-1/2}exp(-\frac{1}{2}\frac{(y_i-x_i\beta_0)^2}{\sigma_0^2})$

위 가정에 따라 샘플들은 모두 독립적이므로, likelihood는 아래와 같이 개별 샘플의 likelihood의 product로 구합니다.

$\begin{aligned} \mathcal{L}(\beta,\sigma^2;y,X)&=\prod_{i=1}^Nf_Y(y_i|X;\beta,\sigma^2) \\&=\prod_{i=1}^N(2\pi\sigma^2)^{-1/2}exp(-\frac{1}{2}\frac{(y_i-x_i\beta)^2}{\sigma^2}) \\&=(2\pi\sigma^2)^{-N/2}exp(-\frac{1}{2\sigma^2}\sum_{i=1}^N(y_i-x_i\beta)^2) \end{aligned}$

여기에 자연로그를 취하면,

$\begin{aligned} \mathcal{l}(\beta,\sigma^2;y,X)&=\ln(\mathcal{L}(\beta,\sigma^2;y,X)) \\&=\ln\bigg((2\pi\sigma^2)^{-N/2}exp(-\frac{1}{2\sigma^2}\sum_{i=1}^N(y_i-x_i\beta)^2\bigg) \\&=\ln\Big((2\pi\sigma^2)^{-N/2}\Big)+\ln\bigg(exp(-\frac{1}{2\sigma^2}\sum_{i=1}^N(y_i-x_i\beta)^2\bigg) \\&=-\frac{N}{2}\ln(2\pi\sigma^2)-\frac{1}{2\sigma^2}\sum_{i=1}^N(y_i-x_i\beta)^2 \\&=-\frac{N}{2}\ln(2\pi)-\frac{N}{2}\ln(\sigma^2)-\frac{1}{2\sigma^2}\sum_{i=1}^N(y_i-x_i\beta)^2 \end{aligned}$

이제 log-likelihood를 maximize하는 $\beta$와 $\sigma^2$을 구하면 우리가 원하는 목적을 달성할 수 있습니다. 다시 말해 log-likelihood의 각각의 파라메터에 대한 그래디언트 값이 0이 되는 $\beta$와 $\sigma^2$값을 찾는 것입니다.

$\max_{\beta,\sigma^2}\mathcal{l}(\beta,\sigma^2;y,X)$

$\nabla_{\beta}\mathcal{l}(\beta,\sigma^2;y,X)=0$

$\frac{\partial}{\partial\sigma^2}\mathcal{l}(\beta,\sigma^2;y,X)=0$

먼저 $\beta$를 구해보면,

$\begin{aligned} \nabla_{\beta}\mathcal{l}(\beta,\sigma^2;y,X) &=\nabla_{\beta}\bigg(-\frac{N}{2}\ln(2\pi)-\frac{N}{2}\ln(\sigma^2)-\frac{1}{2\sigma^2}\sum_{i=1}^N(y_i-x_i\beta)^2\bigg) \\&=\frac{1}{\sigma^2}\sum_{i=1}^Nx^{\intercal}(y_i-x_i\beta) \\&=\frac{1}{\sigma^2}\bigg(\sum_{i=1}^Nx^{\intercal}y_i-\sum_{i=1}^Nx^{\intercal}xi\beta\bigg) \end{aligned}$

위 식이 0이 되도록 하려면 아래의 조건을 만족하면 됩니다.

$\sum_{i=1}^Nx^{\intercal}y_i-\sum_{i=1}^Nx^{\intercal}xi\beta=0$

OLS estimator의 경우와 마찬가지로 Design Matrix가 full rank임을 가정했으므로, $X^{\intercal}X$가 역행렬을 가지게 되어 아래와 같은 결론을 도출할 수 있게 됩니다.

$\beta=\bigg(\sum_{i=1}^Nx^{\intercal}xi\bigg)^{-1}\sum_{i=1}^Nx^{\intercal}y=(X^{\intercal}X)^{-1}X^{\intercal}y$

마지막으로 $\sigma^2$을 구해보면,

$\begin{aligned} \frac{\partial}{\partial\sigma^2}\mathcal{l}(\beta,\sigma^2;y,X)&=\frac{\partial}{\partial\sigma^2}\bigg(-\frac{N}{2}\ln(2\pi)-\frac{N}{2}\ln(\sigma^2)-\frac{1}{2\sigma^2}\sum_{i=1}^N(y_i-x_i\beta)^2\bigg) \\&=-\frac{N}{2\sigma^2}-\bigg[\frac{1}{2}\sum_{i=1}^N(y_i-x_i\beta)^2\bigg]\frac{d}{d\sigma^2}\bigg(\frac{1}{\sigma^2}\bigg) \\&=-\frac{N}{2\sigma^2}-\bigg[\frac{1}{2}\sum_{i=1}^N(y_i-x_i\beta)^2\bigg]\bigg(-\frac{1}{(\sigma^2)^2}\bigg) \\&=-\frac{N}{2\sigma^2}+\bigg[\frac{1}{2}\sum_{i=1}^N(y_i-x_i\beta)^2\bigg]\frac{1}{(\sigma^2)^2} \\&=\frac{1}{2\sigma^2}+\bigg[\frac{1}{\sigma^2}\sum_{i=1}^N(y_i-x_i\beta)^2-N\bigg] \end{aligned}$

$\sigma^2$이 0이 아니라는 가정 하에, 위 유도식이 0이 되는 경우는 오직 아래를 만족하는 경우입니다.

$\sigma^2=\frac{1}{N}\sum_{i=1}^N(y_i-x_i\beta)$

결국 이 경우에도 OLS estimator와 동일함을 알 수 있습니다. 두 파라메터 모두 OLS estimator와 MLE estimator가 동일함은 실제로 아래의 코드에서 확인이 가능합니다.

# Python implementation

위에서 정리한 내용을 파이썬 코드로 구현하면 아래와 같습니다. likelihood function의 경우엔 최대한 수식을 그대로 나타내려고 했습니다.

여담으로 Optimizer 부분을 보면 느끼시겠지만 실제로 scipy에는 다양한 최적화 알고리즘을 선택할 수 있습니다. 제 경우엔 BFGS(Broyden–Fletcher–Goldfarb–Shanno)를 사용해보았는데요, 이 부분을 구현하다가 우연히 제가 존경하는 다크프로그래머님의 최적화 기법의 직관적 이해라는 글을 읽게 되었는데, 다양한 최적화 기법들의 원리와 특징들을 저 같은 초심자를 위해 아주 잘 정리 해 주셨습니다.

import numpy as np
import statsmodels.api as sm
from scipy.optimize import minimize
from scipy.optimize import differential_evolution

df = sm.datasets.get_rdataset("Duncan", "carData").data

y = df['prestige'].values.reshape(45, 1)
X = df[['income', 'education']].values.reshape(45, 2)
params0 = [1., 1., 1]

N = len(X)

# Likelihood function
def loglikelihood(params):
    beta = np.array([params[0], params[1]]).reshape(2, 1)
    sigma_square = params[2]
    
    error = y-np.matmul(X, beta)
    
    return -(-(N/2)*np.log(2*np.pi)-(N/2)*np.log(sigma_square)-(1/(2*sigma_square))*np.sum(error**2))
    
# Optimizer
res = minimize(loglikelihood, params0, method='BFGS', options={'gtol': 1e-4, 'disp': True})
# res = differential_evolution(loglikelihood, [(0, 1), (0, 1), (0, 1000)], tol=1e-4, mutation=(1e-4, 1e-1), disp=True, maxiter=1000)

# Parameters
print(res.x)

"""
[  0.5482716    0.4957514  174.81864984]
"""

# Analytical Beta & Sigma^2

print(np.matmul(np.matmul(np.linalg.inv(np.matmul(X.T, X)), X.T), y))

"""
[[0.54827173]
 [0.49575132]]
"""

print((1/N)*np.sum((y-np.matmul(X, np.array([res.x[0], res.x[1]]).reshape(2, 1)))**2))

"""
174.82031640531054
"""

MLE Estimator로 구한 베타와 시그마 제곱은 OLS Estimator의 경우와 완전히 동일함을 확인할 수 있습니다.

# Reference

• Linear regression models
• Linear regression - Maximum Likelihood Estimation
• Maximum Likelihood Estimation for Linear Regression
• 다크 프로그래머 :: 최적화 기법의 직관적 이해