Maximum Likelihood Estimator(최대가능도추정량)
Maximum likelihood estimator는 frequentist가 중요하게 여기는 또 다른 유형의 estimator로 likelihood가능도를 최대화하여 얻어지는 estimator이다. 이 estimator는 우리가 model로 가정한 probability distribution이 어느 parameter $\theta$ 값이 주어진 data를 가장 잘 설명할 수 있는지를 찾는 것과 같다. 이를 찾는 방법으로써 likelihood를 최대화하는 parameter의 값을 찾는 방식을 사용한다. Likelihood는 이전 글에서도 몇 번 언급했지만, data의 joint distribution을 parameter의 함수로 다시금 바라본 함수이다. Likelihood가 결국 설명하는 것은, 주어진 parameter $\theta$에 대하여 우리가 가정한 probability distribution $P_{\theta}$ 아래에서 이미 수집된 data가 그와같이 나올 가능성, 경향성 정도를 의미한다. 결국 Likelihood를 최대화하는 $\theta$를 찾겠다는 것은 data와 가장 잘 어울리는 probability distribution의 $\theta$값을 찾겠다는 것을 의미한다. 해당하는 $\theta$ 값은 true value의 $\theta$를 estimate하기 위한 performance 좋은 estimator로써 자리매김하는데, 이 estimator의 이름이 바로 maximum liklihood estimator최대가능도추정량 이다. 상당히 직관적인 철학을 바탕으로 정의된 듯한 maximum likelihood estimator가 실질적으로 효율적이고 성능이 좋은 estimator인지 알아보기 위해서는 조금 더 많은 과정을 거쳐야 한다. 하지만, 적어도 이렇게 정의한 estimator는 우리의 믿음과 직관으로 보았을 땐 상당히 타당한 과정에 입각하여 정의되었고, 좋은 성질을 가질 것이라 기대할만도 하다.
먼저 data $X$의 joint density를 $p_{\theta}$라고 하자. 만약 data가 i.i.d.였다면 density $p_{\theta}$는 각 data $X_i$ joint density의 전체 곱으로 나타내어질 것이다. 이 density함수는 아직 data $x$의 함수이다. 하지만 다른 관점으로 $p_{\theta}$를 $\theta$의 함수로 나타낼 수 있고, 이를 likelihood라고 한다.
\[L(\theta) = p\_{\theta}(X)\]그리고 $L(\cdot)$를 최대화하는 $\hat{\theta} = \hat{\theta}(X)$를 $\theta$의 maximum likelihood estimator(MLE)라고 한다. 그리고 $g(\theta)$의 MLE는 $g(\hat{\theta})$로 둔다(이렇게 두어도 원래 MLE 정의와 충돌이 딱히 생기지 않는다). 하지만 likelihood $L(\cdot)$는 주로 곱으로 이루어진 항이 많기 때문에 미분이 어렵다. MLE를 구하기 위해서 logarithm를 씌운 log-likelihood로그가능도함수를 활용하면 미분하기 매우 수월하다.
\[l(\theta) = \log L(\theta)\]이제 exponential family를 density로 가지는 data의 MLE 형태를 살펴보자.
Example 20.1 Data $X$의 density가 canonical one-parameter exponential family로부터 비롯되었을 때, 우리는 density를 다음과 같이 작성할 수 있다.
\[p_{\eta}(x) = \exp\{\eta T(x) -A(\eta)\}h(x)\]MLE를 구하기 위해서 log-likelihood를 살펴보자.
\[l(\eta) = \log p_{\eta}(X) = \eta T - A(\eta) +\log h(X)\]우리는 다음을 얻는다.
\[\begin{align*} l'(\eta) &= T - A'(\eta) \\ l''(\eta) &= -A''(\eta) - Var_{\eta}(T)\lt 0 \\ \end{align*}\]이계도함수가 항상 음수이므로 local minimum이 unique하여 global minimum이 된다. 그러므로 $A’$의 inverse를 $\psi$라고 두면 $\eta$의 MLE는
\[\hat{\eta} = \psi(T)\]이다.
이제 data가 $n$개의 i.i.d. random sample $X_1, \cdots, X_n$으로 이루어져 있고, density $p_{\eta}$를 가지고 있을 때 전체 data $X = (X_1, \cdots, X_n)$의 joint density는 $\prod_{i=1}^n p_{\eta}(x_i)$이고 log-likelihood는
\[l(\eta) = \eta\sum_{i=1}^n T(X_i) - nA(\eta) + \sum_{i=1}^n \log h(X_i)\]이다. 위와 마찬가지 방법으로 MLE를 구하면
\[\hat{\eta} = \psi(\overline{T}) \quad\text{where}\quad \overline{T} = \frac{1}{n}\sum_{i=1}^n T(X_i)\]이다. 이것이 다른 한편으로 보여주는 것은 $E_{\eta} T(X_i) = A’(\eta)$의 MLE가
\[A'(\hat{eta}) = \overline{T}\]즉 그저 $T$의 sample mean이라는 것이다.
한 가지 더 관찰할 수 있는 점은 $\eta$의 MLE가 complete sufficient statistic $T$의 함수라는 것이다. 그러므로 $\eta$의 MLE 역시 UMVUE임을 알 수 있다. 하지만 일반적인 경우에서는 항상 UMVUE는 아니고 심지어 bias를 가지고 있을 수도 있다.
이 상황에서 MLE가 i.i.d. sample들의 mean의 함수이기 떄문에 CLT와 delta method를 이용하여 limiting distribution를 알 수 있다. 먼저 $Var_{\eta}(T) = A’’(\eta)$이기 때문에 CLT에 의하여
\[\sqrt{n}(\overline{T} - A'(\eta)) \leadsto N(0, A''(\eta))\]이다. 이제 $\psi(\cdot)$로 delta method를 적용해야 한다. $\psi$의 도함수는 inverse function derivative를 구하는 방법을 참고하면
\[\psi' = \frac{1}{A''(\psi(\cdot))}\]이다. 그러므로 $\psi’(A’(\eta)) = 1/A’’(\eta)$이고, delta method에 의해
\[\sqrt{n}(\hat{\eta} - \eta) \leadsto N(0, 1/A''(\eta))\]이다. 한편, 각 data $X_i$가 전달하는 Fisher information은
\[I(\eta) = Var_{\eta}\left(\frac{\partial}{\partial\eta}\log p_{\eta}(X_i)\right) = Var_{\eta}(T(X_i)) = A''(\eta)\]이므로, 전체 data에 대한 Fisher information은 $J(\theta) = nI(\theta) = nA’’(\eta)$이며, Cramer-Rao bound에 의하여 $\eta$의 unbiased estimator $\tilde{\eta}$에 대하여
\[Var_{\eta}(\tilde{\eta}) \ge \frac{1}{nA''(\eta)}\]임에 따라
\[Var_{\eta}(\sqrt{n}(\tilde{\eta} - \eta)) = nVar_{\eta}(\tilde{\eta})\ge \frac{1}{A''(\eta)}\]가 성립한다. 그러므로 MLE $\hat{\eta}$의 variance는 $n\to\infty$임에 따라 가장 최소화됨을 알 수 있다. 즉, asymptotic점근적으로하게 Cramer-Rao bound에 도달한다. 이러한 형태의 estimator들을 asymptotically efficient하다고 한다.
다른 model에 대해서도 MLE를 한 번 구해보자.
Example 20.2 Data $Y_1, \cdots, Y_n$이 서로 independent하게 $Y_i \sim N(\alpha+\beta x_i, \sigma^2),\, i = 1, \cdots, n$을 따른다고 하고, $x_1, \cdots, x_n$은 알려진 constant라고 가정하자. 이는 사실 가장 단순한 regression회귀 model임을 알 수 있다. 이때, covariate 혹은 explanatory variable, independent variable로 불리는 known parameter $x_i$가 함께 model을 구성하고 있음을 알 수 있다. 만약 data와 연관이 깊은 관측가능한 parameter가 발굴되었을 때, 그 parameter를 covariate으로 활용하면 data에 대해 더 많은 정보를 전달할 수 있을 것이다. Regression은 바로 이 점을 이용하여 만들어진 model인데, 특히 각 data가 따르는 distribution의 mean을 covariate을 통해 linear하게 조절할 수 있게 함으로써 data에 더 적합한 modeling을 할 수 있도록 도와주었다. 그러므로 고작 두 개의 unknown parameter로 distribution의 mean이 설명된다고 가정할 수 있게 되었고, 충분히 정확도 높은 modeling을 구사할 수 있게 되었다.
이 상황에서 MLE를 한 번 구해보자. 어렵지 않게 $Y_i$를 data라고 생각하면서 과정을 밟아나가면 된다. 이때 unknown parameter는 $\theta = (\alpha, \beta, \sigma^2)$이다. 각 data point가 따르는 distribution의 density는
\[p_{\theta}(y_i) = (2\pi\sigma^2)^{-1/2}\exp\left[-\frac{1}{2\sigma^2}(y_i - \alpha -\beta x_i)^2\right]\]각 data는 서로 independent하므로 전체 data의 joint density는 각 data point density의 곱이다. 그러므로 likelihood는
\[L(\theta) = \prod_i p_{\theta}(y_i) = (2\pi\sigma^2)^{-n/2}\exp\left[-\frac{1}{2\sigma^2}\sum_i (y_i -\alpha-\beta x_i)^2\right]\]그리고 log-likelihood는
\[l(\theta) = \log L(\theta) = -\frac{n}{2}\log{(2\pi)} -\frac{n}{2}\log{\sigma^2}-\frac{1}{2\sigma^2}\sum_i(y_i-\alpha-\beta x_i)^2\]이다. 이를 통해 partial derivative를 살펴보면
\[\begin{align*} \frac{\partial l}{\partial\sigma^2}\bigg\vert_{\hat{\theta}} &= -\frac{n}{2\sigma^2}+\frac{1}{2\sigma^4}\sum_i (y_i -\alpha-\beta x_i)^2 =0 \\ \frac{\partial l}{\partial\alpha}\bigg\vert_{\hat{\theta}} &= -\frac{1}{\sigma^2}\sum_i(y_i -\alpha-\beta x_i) = 0 \\ \frac{\partial l}{\partial\beta}\bigg\vert_{\hat{\theta}} &= -\frac{1}{\sigma^2}\sum_i x_i(y_i-\alpha-\beta x_i) = 0 \\ \end{align*}\]각 derivative에 대하여 endpoint끝점으로 limit을 취했을 때 lower endpoint으로의 limit은 $-\infty$, upper endpoint으로의 limit은 $+\infty$인 것으로 보아 $l(\cdot)$은 global maximum을 가지는 것을 알 수 있다. 그리고 그 maximum point는 위의 연립된 식을 풀면 도출된다. 그러므로 MLE는
\[\begin{align*} \hat{\beta} &= \frac{\sum x_i\sum y_i - n\sum_i {x_iy_i}}{(\sum x_i)^2 - n\sum x_i^2} \\ \hat{\alpha} &= \frac{1}{n}(\sum y_i - \hat{\beta}\sum x_i) \\ \hat{\sigma^2} &= \frac{1}{n}\sum_i (y_i - \hat{\alpha} - \hat{\beta}x_i)^2 \end{align*}\]이다. 식이 길어지는 탓에 모두 대입하여 표현하진 않았다. $\sigma^2$은 각 data가 mean으로부터 어느 정도의 퍼짐으로 값을 가지는지, 혹은 mean으로부터 어느 정도의 error오차을 가지는지를 표현하는 항이다. 이를 조금 더 명확히 하면 error term오차항이라고 하는 $\epsilon_i := Y_i - \alpha - \beta x_i \sim N(0, \sigma^2)$는 data가 가질 error를 random variable로 나타낸 것인데, 이의 variance가 $\sigma^2$로 세워졌음을 확인할 수 있다. $\sigma^2$의 MLE는 $y_i - \hat{\alpha} - \hat{\beta}x_i$의 제곱의 합 형태로 도출이 된 것을 확인할 수 있는데, $y_i - \hat{\alpha} - \hat{\beta}x_i$를 주로 $e_i$라고 표기하며 residual잔차라고 일컫는다. 그러므로 error의 variance의 MLE는 residual들의 sample variance(그런데 나누는 수가 $n-1$이 아닌 $n$)로 꽤나 자연스럽게 도출되었음을 알 수 있다.
Example 20.3 Data $X_1, \cdots, X_n$이 jointly distributed되어 있어, 첫 번째 data는 $X_1\sim N(0, 1)$이었고, $j = 1, \cdots, n-1$에 대하여
\[X_{j+1}\vert X_1 = x_1, \cdots, X_j = x_j \sim N(\rho x_j, 1)\]이라고 하자. $\rho$는 unknown parameter이다. 이는 autoregression자기회귀 model 중 가장 간단한 경우로, 바로 직전(즉, lag가 1인 경우)의 data에 일정한 rate로 mean이 변화하는 구조를 띠고 있다. Autoregression은 data point가 independent한 것이 아니라 시간 등에 따라 어느 정도의 correlation을 가져서 앞의 data의 값에 따라 영향을 받는 구조에 적합한 model이다. 여기서 Likelihood를 구해보면
\[\begin{align*} L(\rho) &= p_{\rho}(x_1)\prod_{j=1}^{n-1} p_{\rho}(x_{j+1}\vert x_1, \cdots, x_j) \\ &= (2\pi)^{-n/2}\exp\left[\frac{-1}{2}\left(x_1^2 + \sum_{j=1}^{n-1}(x_{j+1}-\rho x_j)^2\right)\right] \end{align*}\]이고, log-likelihood는
\[l(\rho) = -\frac{n}{2}\log{(2\pi)}-\frac{1}{2}\left(x_1^2 +\sum_{j=1}^{n-1}(x_{j+1} - \rho x_j)^2\right)\]이다.
\[\begin{align*} \frac{\partial l}{\partial \rho}\bigg\vert_{\hat{\rho}} &= \sum_{j=1}^{n-1} x_j(x_{j+1} - \rho x_j) = 0 \\ \frac{\partial^2 l}{\partial\rho^2}\bigg\vert_{\hat{\rho}} &= -\sum_{j=1}^{n-1} x_j^2 \lt 0 \end{align*}\]이므로 우리는 다음과 같은 MLE를 얻는다.
\[\hat{\rho} = \frac{\sum_{j=1}^{n-1} x_jx_{j+1}}{\sum_{j=1}^{n-1} x_j^2}\]
MLE의 limiting distribution은 굉장한 성질이 있기에 나중에 여러 편에 걸쳐서 소개하게 될 것이다. 또한 asymptotically efficient한 estimator들의 property를 다루기도 쉽지 않은 편이라 이 시리즈의 뒷편에서 다루도록 한다.
다음 편에서는 median이 estimator로써의 역할을 하게 되면 어떠한 성질들이 있는지 살펴본다. 특히 median의 limiting distribution과 performance를 정량적으로 살펴보게 될 것이다.