본문으로 바로가기



감마(Gamma) 분포 MLE 추정 with R - 3. 최대 가능도 추정량(Maximum Likelihood Estimator) 구하기

category Statistics/Inference 2017.06.30 00:21
감마(Gamma) 분포 MLE 추정 with R - 3. 최대 가능도 추정량(Maximum Likelihood Estimator) 구하기

이번 포스팅에서 사용될 R 패키지는 다음과 같다.

library(tidyverse)
library(magrittr)
library(RColorBrewer)

감마 분포(Gamma distribution)의 MLE

우리가 이전 포스팅에서 우리는 MOM(Method of moments)을 이용한 모수 추정에 대하여 알아보았다. 이번 시간에는 최대 가능도 추정법을 사용하여 감마 분포의 모수를 추정하는 방법에 대하여 알아보도록 하겠다.

MOM방법을 이용하여 모수에 대하여 추정하는 논리를 이해하기 시작하면, 통계가 재미있어 진다. 그리고 MOM방법의 업그레이드 버젼이자 통계학 추정이론의 꽃이라 할 수 있는 최대 우도 추정량(Maximum Likelihood Estimator)을 구하는 방법을 이해한다면, 통계가 가지는 진짜 재미에 빠져들게 될 것이라 생각한다.

감마분포 데이터, 모수는??

오늘도 역시 데이터가 감마분포에서 나왔다는 전제하에, 과연 이 데이터는 감마분포 어떤 모수에서 나온 것인지에 대한 퀴즈를 최대 가능도 추정법을 통하여 풀어보는 방법에 대하여 알아보자. 모수의 정확한 값은 역시 필자만 알고있으며, 포스팅의 후반부에 공개하도록 하겠다.

필자는 감마분포에서 다음의 30개의 샘플을 뽑았다. 단, shape 모수인 k와 scale 모수인 theta는 필자만 알고 있도록 하겠다.(이 포스팅의 끝에 공개하도록 하겠다.)

set.seed(1234)
my_sample <- rgamma(30, shape = k, scale = theta)
my_sample
##  [1]  2.409058  6.169255  6.306661  3.400211  6.766151  3.761012  8.372326
##  [8]  3.648530  3.996274  1.486417  5.432445  8.287809  4.945073  5.301527
## [15]  7.415598  3.591295  5.633027  3.963306  5.264060  5.144931  1.972472
## [22]  5.282026  5.277706  7.325914  8.536704 19.977603 10.843007  1.925971
## [29] 10.105529  2.679378

최대 가능도 함수(likelihood function) 유도

이전 포스팅에서 알아봤던 감마분포의 확률밀도함수(pdf)의 형태는 다음과 같았다.

\[ f(x;k,\theta) = {\frac {1}{\Gamma (k)\theta ^{k}}}x^{k\,-\,1}e^{-{\frac {x}{\theta }}} \]

\(n\)개의 관측값이 있다는 가정하(우리는 현재 30의 관측치가 있다.)에서 \(k\)\(\theta\)를 모수로 갖는 감마 분포의 가능도 함수는 다음과 같이 각 데이터에서의 pdf값을 곱한 형식을 취하고 있다.

\[ \begin{align*}L(k,\theta;\underline{x}) & =\prod_{i=1}^{n}f(x_{i};k,\theta)\\ & =\frac{1}{\Gamma(k)^{n}\theta^{nk}}\left(\prod_{i=1}^{n}x_{i}\right)^{k-1}e^{-\sum_{i=1}^{n}\frac{x_{i}}{\theta}} \end{align*} \]

왜 가능도 함수 인가?

그런데 왜 가능도 함수라고 부를까? 만약 분포가 연속형 분포가 아닌 이산형(discrete) 분포라고 한다면 가능도보다 더 피부에 와닿는 말은 최대 확률분포 함수가 아닐까 싶다. 하지만 통계시간에 배웠다시피 연속형 분포의 pdf가 특정 포인트에서 갖는 확률은 0이다. 따라서 ’확률’이라는 용어를 사용하여 함수를 표현하기에는 무리가 있다. 하지만 데이터 값을 중심으로 일정한 \(dx\)값(예를 들어 0.001같은)을 곱해준다고 가정해보면 사실 가능도 함수는 실제 데이터가 나올 확률를 훌륭하게 근사하고 있다고 볼 수 있다.

파란선은 확률이 0지만, 빨간 박스는 확률이 존재한다.

파란선은 확률이 0지만, 빨간 박스는 확률이 존재한다.

가능도 함수는 x대한 함수가 아니다.

많은 사람들이 MLE를 배우면서 가능도 함수에 대하여 착각한다. 가장 대표적인 것 중 하나가 바로 가능도 함수가 앞의 식에서 본 $x_i$에 대한 함수라고 생각하는 것이다. 분명히 해야할 것은 위의 식에서 우리가 변하게 만드는 값은 k$\theta$이지 $x_i$들이 아니다. \(x_i\)들은 필자가 여러분에게 준 데이터의 값들, 예를 들어 2.409058, 6.169255 같은 주어진 숫자들이며, 우리가 조절할 수 없는, 변하지 않는 값이라는 사실을 이해하도록 하자. 이러한 \(x_i\)값들이 주어졌다는 것을 나타내기 위하여 위의 식에서도 ;을 사용하고 있다.

다음의 주어진 sample에 대한 가능도 함수를 보면 필자가 이야기한 부분을 좀 더 이해할 수 있을 것이다.

loglik_gamma <- function(k, theta, given_sample = my_sample, log = TRUE){
  result <- dgamma(given_sample, shape = k, scale = theta, log = TRUE) %>% sum()
  if (log == FALSE) exp(result) else result
}
loglik_gamma(1, 3)
## [1] -91.36546
loglik_gamma(2, 3)
## [1] -75.88655

위의 두 값을 기준으로 판단해보았을 경우, 모수 (k = 1, theta = 3)의 조합보다는 (k = 2, theta = 3)의 조합이 데이터의 모분포로 훨씬 좋아보인다. 왜냐하면 두번째 모수 조합의 분포에서 my_sample이 나올 확률이 첫번째 모수 조합의 분포에서 나왔을 확률보다 높기 때문이다. 이런 방식을 사용하면 결국 우리는 가능한 모든 조합들 중에서 가장 그럴싸한(주어진 데이터와 가장 잘 어울리는) 모수의 조합을 선택하고 싶어진다. 결국 이 과정을 수행하기 위하여 가능도 함수의 maximum값을 찾아야만 한다. 그리고 그때의 입력값(모수의 조합)이 우리가 찾는 MLE(maximum likelihood estimator) 된다.

가능도 함수(likelihood function) 그리기

자, 이제 주어진 my_sample에 대한 감마 분포의 가능도 함수(likelihood function)를 그려보도록 하자. 또한 MLE를 구하기 앞서 이전 포스팅에서 구했던 MOME는 가능도 함수 어디쯤에 위치하고 있는지 그래프에 나타내어보도록 하자. 아래의 3D plot은 가능도 함수를 적절한 스케일링을 통하여 전체적인 모습을 나타낸 것이며, 빨간색의 구는 MOM 방법으로 구한 모수의 가능도 함수 값을 나타낸다.

library(rgl)
library(MASS)

# scaled likelihood function
lik_gamma <- function(k, theta, given_sample = my_sample){
  result <- dgamma(given_sample, shape = k, scale = theta, log = TRUE) %>% sum()
  10^32 * exp(result)
}
lik_f <- Vectorize(lik_gamma)

# draw likelihood function
persp3d(lik_f, 
        xlim = c(2, 5), ylim = c(1, 4), zlim = c(0, 1),
        n = 100)

# add MOME (http://www.issactoast.com/96)
mom_point <- mom_gamma(my_sample)

# 3d sphere add
rgl.spheres(x = mom_point$k_hat,
         y = mom_point$theta_hat,
         z = lik_gamma(mom_point$k_hat, mom_point$theta_hat),
         r = 0.05, color = 'red')
text3d(x = mom_point$k_hat,
       y = mom_point$theta_hat + 0.4,
       z = lik_gamma(mom_point$k_hat, mom_point$theta_hat),
       "MOME")
rglwidget()