본문으로 바로가기



감마(Gamma) 분포 MLE 추정 with R - 2. 적률법(MOM)을 통한 모수 추정

category Statistics/Inference 2017.06.28 10:40
감마(Gamma) 분포 MLE 추정 with R - 2. 적률법(MOM)을 통한 모수 추정

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

library(tidyverse)
library(RColorBrewer)

감마 분포(Gamma distribution)의 MOME와 MLE

우리가 이전 포스팅에서 우리는 감마분포의 모수들이 어떠한 의미를 가지고 있는지에 대하여 알아보았다. 오늘은 데이터가 감마분포에서 나왔다는 전제하에(사실 이 부분은 중요한 부분이다. “주어진 데이터를 가지고 감마분포의 MLE로 모수를 추정했습니다.” 라고 말하기 전에 사실 우리가 가진 데이터가 감마분포를 따르고 있는지에 대하여부터 판단을 해야한다. 하지만 실제 수리 통계 시간에 배우는 순서는 좀 더 근본적인 개념인 추정량을 먼저 배운다.), 과연 이 데이터는 감마분포 어떤 모수에서 나온 것인지에 대한 퀴즈를 푸는 방법(MOME와 MLE)에 대하여 알아보자.

감마분포에서 나온 데이터

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

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

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

set.seed(1234)
my_sample <- rgamma(50, shape = k, scale = theta)
my_sample
##  [1]  1.9114268  6.0435346  6.2007164  2.9568366  6.7283355  3.3476757
##  [7]  8.5933373  3.2253450  3.6048406  0.9962285  5.2059141  8.4945011
## [13]  4.6573281  5.0580900  7.4788598  3.1632630  5.4330140  3.5687013
## [19]  5.0158459  4.8817107  1.4690227  5.0360996  5.0312285  7.3749086
## [25]  8.7857642  4.3043800  6.4729458  4.7796368  2.6609239  3.5614977
## [31]  2.4062988  2.0879139  3.4807080  7.2663264  5.5315687  2.7274018
## [37]  2.3054823 17.2093711  4.9759806  6.1425408  4.4050257 10.0883243
## [43]  6.1206274  5.0218208  3.6634073  6.0811815  5.5390992  5.6314484
## [49]  2.9253241  2.9736188

여러분의 주어진 임무(?)는 위의 my_sample이 어떤 감마분포에서 왔는지, ktheta값을 맞추는 것이다. 만약 필자가 여러명에게 같은 질문을 했을 경우, 그 중 별 관심없는 학생 한 명은 다음과 같이 찍을 수도 있겠다.

‘안봐도 알죠, k = 1, theta = 3 입니다.’

위와 같은 대답은 사실 합리적인(?) 사고를 가진 사람들에겐 한낮 웃기려고 한 말에 불과하겠지만, 사실 위의 문장은 엄연한 감마 분포 모수에 대한 추정량(estimator)이라고 말할 수 있다. 실제로 필자가 정말 k = 1, theta = 3인 감마분포에서 위의 데이터를 발생시켰다면, 제일 정확한 추정량이 될 것이다.

하지만, 통계학자들은 이 문제를 푸는 방법에 있어서 다음의 대표적인 두 가지를 제안했는데, 첫번째 방법은 MOME(Method of moments estimator, 한국어로는 적률법이라 한다.)를 사용하는 방법이고, 두번째 방법은 MLE(Maximum likelihood estimator, 최대 가능도 추정량)를 이용하여 푸는 방법이다.

오늘은 이 두 가지 방법 중 첫번째 방법인 Method of moments에 대하여 공부하도록 하자.

감마 분포의 MOME

영어로 Moments라 하지만 통계에서는 그냥 확률변수의 자승들의 기대값을 통칭하는 말이라고 생각하면 된다. (물리학에서 비롯된 개념이라 이름에 moments가 들어가지만 사실 개인적으로는 이름을 잘못 지었다고 생각한다). MOM의 핵심 아이디어는 우리가 가진 데이터가 분포를 잘 반영한다면, 데이터를 가지고 계산한 표본 평균과 표본 분산값들이 이론값과 일치할 것이다라는 믿음에 기반하고 있다.

위의 주어진 pdf를 가지고 감마분포를 따르는 확률변수의 평균을 구해보면 두 모수의 곱으로 나타난다. 즉, 통계 기호를 좀 써보면 확률변수 \(X\sim Gamma\left(k,\theta\right)\)의 기대값은 다음과 같다.

\[ \mathbb{E}\left[X\right]=k\theta \] 따라서 만약 정말 데이터가 감마분포에서 나왔다면, 표본평균 값이 k*theta의 값과 비슷해야 하지 않을까?

\[ \overline{X}\approx k\theta \] 따라서 MOME에서는 주어진 데이터의 표본 평균을 두 모수의 곱과 같다고 추측한다.

mean(my_sample) #(set) == k * theta
## [1] 5.052508

하지만 표본 평균 값 5.0525076 하나만을 가지고 우리가 원하는 k값과 theta값을 구하는 것은 불가능하다. 왜냐하면 알고싶은 변수의 숫자(2개)보다 알고있는 정보(곱해서 5.0525076)가 적기 때문이다. 따라서 정보를 하나 더 얻어야하는데, 보통 평균 다음으로 만만한 분산을 사용한다. 감마 분포의 분산은 다음과 같다.

\[ Var\left(X\right)=\mathbb{E}\left[X^{2}\right]-\left(\mathbb{E}\left[X\right]\right)^{2}=k\theta^{2} \] 분산의 정의식을 사용하지 않고, 확률변수의 제곱의 기대값과 기대값을 사용한 이유는 앞에서 언급한 moments, 즉, 기대값들로 나타내야하기 때문이다. (사실 분산의 정의와 표본 분산을 사용해도 무방하다. 이 경우 중심적률법이라 따로 명명하는 것으로 알고 있다.) 이 식 역시 이론상 성립해야 한다면, 우리는 다음의 값이 비슷할 것이라 충분히 생각할 수 있다.

\[ \overline{X^{2}}-\left(\overline{X}\right)^{2}\approx k\theta^{2} \]

따라서 앞에서와 마찬가지로 MOME에서는 다음의 식을 같다고 생각한다.

mean(my_sample^2) - mean(my_sample)^2 #(set) == k * theta^2
## [1] 6.959629

자, 이제 모수 두 개에 대한 정보를 사용하여, 위에서 설정한 두 식을 모두 만족하는 ktheta를 기대값으로 나타낼 수 있는데, 그 중 \(\theta\)를 나타내어보면 다음과 같다.

\[ \theta=\frac{\theta^{2}k}{\theta k}=\frac{Var\left(X\right)}{\mathbb{E}\left[X\right]}=\frac{\mathbb{E}\left[X^{2}\right]-\left(\mathbb{E}\left[X\right]\right)^{2}}{\mathbb{E}\left[X\right]} \] MOME는 바로 위의 식의 기대값 부분을 데이터(샘플)에서 얻는 정보를 사용하여 근사시키는 것이라고 생각할 수 있다. 즉, 결국 MOM은 기대값을 매치시켜서 모수를 추측하는 방법이라 말할 수 있다. 위를 바탕으로 구한 감마분포의 MOME(MOM Estimator, MOM방법을 사용한 추정량)는 다음과 같다.

\[ \hat{\theta}_{MOM}=\frac{\overline{X^{2}}-\left(\overline{X}\right)^{2}}{\overline{X}} \\ \hat{k}_{MOM}=\frac{\left(\overline{X}\right)^{2}}{\overline{X^{2}}-\left(\overline{X}\right)^{2}} \] 위의 식을 바탕으로 감마분포의 MOME를 구하는 함수를 정의하면 다음과 같다.

mom_gamma <- function(sample_g){
    x_bar  <- mean(sample_g)
    x2_bar <- mean(sample_g^2)
    
    mom_theta <- (x2_bar - x_bar^2) / x_bar
    mom_k <- x_bar / mom_theta
    
    list(k_hat = mom_k, theta_hat = mom_theta)
}

자, 이제 mom_gamma()를 사용하여 우리가 알고싶어 하던 정답을 맞혀보자.

mom_gamma(my_sample)
## $k_hat
## [1] 3.667988
## 
## $theta_hat
## [1] 1.37746

MOM을 사용하여 모수를 추측하는 사람의 경우 k는 3.840, theta는 1.782라고 대답을 하는 것이다. 그렇다면 실제 답은 어떨까? 필자가 설정해 놓은 모수의 값은 다음과 같다.

k; theta
## [1] 3
## [1] 2

독자들 중에는 결과값을 보고 실망을 하는 사람도 있을 것이고, 긍정적인 마인드의 독자들은 적어도 ‘안봐도 알죠, k = 1, theta = 3 입니다.’라고 찍었던 학생의 답보다는 훨씬 정답에 가깝다고 생각할 수 도 있다. 만약 데이터의 갯수가 늘어나서 모집단에 대한 정보가 늘어날 경우에는 어떨까? 다음의 코드로 이번 포스팅을 마치며, 다음 시간에는 MOM보다 더 나은 방법이라 여겨지는 MLE를 통한 문제풀기에 도전하자.

set.seed(1234)
my_sample2 <- rgamma(1000, shape = k, scale = theta)
mome <- mom_gamma(my_sample2)
mome
## $k_hat
## [1] 3.104176
## 
## $theta_hat
## [1] 1.985311



SHARE TO



티스토리 툴바