본문으로 바로가기



감마(Gamma) 분포 MLE 추정 with R - 4. MOME vs. MLE 의 비교 Bias와 MSE에 대하여

category Statistics/Inference 2017.07.10 11:09
감마(Gamma) 분포 MLE 추정 with R - 4. MOME vs. MLE 의 비교 Bias와 MSE에 대하여

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

library(magrittr)

이전 세 편의 연재를 통하여 우리는 기대값 매칭을 통하여 모수를 찾는 적률법(Method of moments)과 데이터의 가능도 함수를 최대화 시키는 모수를 찾는 최대 가능도법(Maximum likelihood )에 대하여 알아보았다. 오늘은 이 방법을 통하여 찾은 추정량이 얼마나 좋은지, 즉 두 방법의 성능을 비교해보는 포스팅이 되겠다.

추정량(estimator)은 확률변수

통계 공부를 하는 많은 사람들이 처음에 추정량(estimator)과 추정값(estimate)을 구분하는데 어려움을 겪는다. 대부분의 경우, 데이터가 주어진 상태에서 모델을 설정하고, 모수에 대한 추정량의 추정값을 구하는 과정을 추정량을 구한 것이라 착각한다. 좀 더 구체적인 예를 들어, 생각해보자.

이전 포스팅에서 우리는 감마분포의 모수에 대한 추정값을 MOME(Method of moments estimator), 즉, 적률추정량을 통하여 구하는 방법에 대하여 알아보았다.

# devtools::install_github("issactoast/r4issactoast")
library(r4issactoast)
set.seed(2017)
my_sample <- rgamma(30, shape = 3, scale = 2)
mom_gamma(my_sample)
##     k_hat theta_hat 
##  3.289676  1.869870

이 과정에서 우리는 mom_gamma 함수를 정의한 후 그 함수에 샘플을 입력값으로 넣어서 k는 3.29, $\theta$는 1.87과 같은 모수들에 대한 추정값을 얻었다. 위의 경우 통계학에서 말하는 추정량(estimator)은 우리가 정의한 mom_gamma() 함수를 의미하게 되고, my_sample 데이터를 기반으로 계산이 되어 나온 mom_gamma 함수의 결과값을 추정값(estimate)이라고 부른다. 여기서 확실히 짚고 넘어가야 할 사실은 주어진 데이터가 변하는 경우, 추정량(mom_gamma 함수)은 매번 다른 추정값을 만들어낸다는 사실을 이해한다면, 추정량(estimator)이 확률변수라는 사실 이해하기 쉬울 것이다.

MOME와 MLE의 성능비교

자 그렇다면 우리가 배웠던 모수의 추정량인 MOME와 MLE 중 어떤 추정량을 사용해야 할까? 어느 한쪽이 다른 한쪽보다 우수할까? 당연히 MLE가 우수한거 아닌가요? 추정량이 확률변수라는 점을 고려하면, 위의 대답은 사실 그리 간단하지 않다. 왜냐하면 샘플 A를 사용해서 MOME와 MLE의 추정값을 구할 경우와 샘플 B를 사용해서 MOME와 MLE의 추정값(이제부터는 MOMe(MOM estimate), MLe(ML estimate)이라 사용하겠다.)을 구하는 경우에 따라 좋은 쪽과 나쁜 쪽이 바뀔 수 있기 때문이다. 간단한 실험을 한번 해보도록 하자.

set.seed(2017) # for verification

# Draw 100 time samples from gamma dist. size of 50
my_samples <- replicate(1000, rgamma(1000, shape = 3, scale = 2))

# Calculate MOMe and MLe
result_mom <- apply(my_samples, 2, mom_gamma)
result_mle <- apply(my_samples, 2, mle_gamma)

위의 코드는 모수 k가 3, theta가 2인 감마분포에서 사이즈가 50인 샘플을 100번 뽑아서 각 샘플을 기준으로 MOM와 ML방법으로 모수의 값을 추정한 것이다. 이렇게 계산된 추정값들을 가지고 모수의 참값에서 얼마나 떨어져있는지를 유클리디안 거리를 계산하여 측정할 수 있다.

# Calculate Euclidian distance from true parameter c(3, 2)
mom_dist <- sweep(result_mom, 1, c(3, 2), "-")^2 %>%
              colSums() %>%
              sqrt()
mle_dist <- sweep(result_mle, 1, c(3, 2), "-")^2 %>%
              colSums() %>%
              sqrt()

# percentage of MOM beating MLE
mean(mom_dist < mle_dist)
## [1] 0.389

위의 결과가 의미하는 것은 40% 정도는 MOMe가 MLe 보다 실제 모수값에 더 가깝게 추정을 했다는 의미이다. 실제 샘플값을 몇개 체크해보자.

# Take a look at some cases:
# MOME beats MLE w.r.t distance from true parameter
momwin_index <- which(c(mom_dist < mle_dist) == TRUE)
result <- rbind(result_mom, result_mle)
row.names(result) <- c("mom_k", "mom_theta", "mle_k", "mle_theta")
head(t(result[, momwin_index]))
##         mom_k mom_theta    mle_k mle_theta
## [1,] 2.981216  2.035823 2.933311  2.069071
## [2,] 2.931271  2.123398 2.864626  2.172799
## [3,] 3.186029  1.879386 3.311943  1.807935
## [4,] 3.009847  1.994634 2.908903  2.063851
## [5,] 3.200642  1.919525 3.205993  1.916321
## [6,] 2.990417  1.999269 2.891512  2.067654

6번째 줄에 위치한 샘플의 경우에는 MOME의 추정값이 MLE의 추정값보다 훨씬 정확한 경우도 있다는 것을 알 수 있다. 그렇다면 도대체 어떠한 방식을 선택해야하는 것일까?

좋은 추정량의 판단 기준

Fig1. Bias와 Variation (photo from https://home.ubalt.edu/ntsbarsh)

Fig1. Bias와 Variation (photo from ‘https://home.ubalt.edu/ntsbarsh’)

치우침(Bias)

통계학에서는 여러 추정량들을 비교하기 위하여 여러가지 기준을 선정해 놓았는데, 그 중 아마도 독자들에게 가장 유명한 기준은 바로 Bias, ’치우침’이다. 어떠한 추정량이 치우침이 없다는 말은 추정량의 기대값이 모수와 일치한다는 말이다. 위의 그림에서 왼쪽 위의 그림과 오른쪽 아래 그림을 보면, 같은 variantion을 가지고 있을 경우, 그 분포의 중심이 실제값과 일치하는 경우가 훨씬 좋은 추정량이 된다.

즉, 어떤 추정량이 좋은 추정량이 되려면 적어도 그 중심이 실제값과 일치하는 것이 좋을 것이다. 이러한 정도를 측정하는 기준이 바로 치우침(Bias)인데, 이것을 수학적으로 나타내면 다음과 같다.

\[ \operatorname {Bias} _{\theta }[\,{\hat {\theta }}\,]=\operatorname {E} _{\theta }[\,{\hat {\theta }}\,]-\theta \]

즉, 추정량의 기대값에서 실제 모수값을 뺀 값이다. 이런 측면에서 보았을 때, MOM 방법은 아주 좋은 추정량이다. 왜냐하면 MOM 자체가 unbiased가 되도록 추정량을 계산하기 때문에(MOME 계산 과정은 이전 포스팅을 참고하자), 이론상의 치우침값은 0이기 때문이다. 이러한 치우침이 없는 추정량을 불편(unbiased) 추정량(estimator)라고 한다. 우리의 예제에서, MOM의 치우침은 다음의 코드로 확인해볼 수 있다.

result_mom %>% apply(1, mean) - c(3, 2)
##        k_hat    theta_hat 
##  0.013535005 -0.002239513

예상과 같이 실제 모수의 값과 거의 일치하는 것을 알 수 있다. 그렇다면 MLE은 어떨까? 일반적으로 MLE은 불편성을 만족하지 않는 추정량으로 알려져있다. 하지만 위의 계산을 mle_gamma를 사용하여 똑같이 적용시켜보면 다음과 같다.

result_mle %>% apply(1, mean) - c(3, 2)
##         k_hat     theta_hat 
##  0.0078523602 -0.0002244214

MLE도 실제 모수에 가까우며, 오히려 MOME보다 더 좋은 성능을 보여준다. 그렇다면 MLE는 불편추정량도 아닌데 어째서 더 좋은 성능을 내는 것일까?

추정량을 판단하는 기준 - MSE(mean squared error)

사실 많은 사람들이 알고있는 소문난 잔치집 격인 bias라는 기준은 좋은 추정량을 판단하는데에 좋은 기준이 아니다(먹을게 없다). 이것의 근본적인 이유는 사실 추정량이 확률변수라는 사실에서 부터 출발하는데, 불편성을 가지고 추정량을 판단하는 방식은 어떤 확률변수를 성질을 판단할 때 평균 하나만을 가지고 이야기하는 것과 마찬가지이기 때문이다.

앞에서 본 Fig1 의 오른쪽에 있는 두 경우를 생각해보자. 추정량 A는 치우침이 없는 반면 분산이 큰 경우와, 추정량 B는 치우침이 살짝 있지만, 분산이 작은 경우에는 추정량 B가 훨씬 안정적으로 추정을 할 것이다.

따라서 추정량의 좋고 나쁨을 판단하는 기준은 추정값의 평균과 분산을 동시에 고려할 수 있는 측정도구를 선택하는게 좋은데, 이에 대표적인 도구가 MSE(mean squared error)이다.

\[ \begin{align} \operatorname{MSE} ({\hat {\theta }}) &= \mathbb {E} \left[({\hat {\theta }}-\theta )^{2}\right] \\ &= \operatorname {Var} ({\hat {\theta }})+\operatorname {Bias} ({\hat {\theta }},\theta )^{2}. \end{align} \]

위의 식을 바탕으로 위에 주어진 데이터를 사용하여 감마함수 모수에 대한 MOME와 MLE의 MSE를 계산해보면 다음과 같다.

sweep(result_mom, 1, c(3, 2), "-")^2 %>% rowMeans()
##      k_hat  theta_hat 
## 0.02362443 0.01238715
sweep(result_mle, 1, c(3, 2), "-")^2 %>% rowMeans()
##       k_hat   theta_hat 
## 0.015940933 0.008660383

MSE의 측면에서 MLE는 MOME보다 더 우위에 있는 추정량으로 보인다. 이것을 시각적으로 나타내어 좀 더 잘 이해해보도록 하자.

언제나 데이터의 시각화는 생각하지 못했던 것들을 보여준다. 사실 필자도 그래프를 그려보기 전에는 k와 theta의 추정값이 서로 상관성이 있다것을 인지하지 못하고 있었지만, 그래프를 그려본 후, 구하는 과정을 생각해보고 결과값이 당연하다고 생각하였다. 이유는 추정과정에서 모수 theta에 대한 추정값을 구한 후, 그것을 바탕으로 k에 대한 추정값을 구하도록 되어있으므로, 당연히 둘 사이에는 상관성이 존재해야한다. 좀 더 통계적으로 말하면, k의 추정량은 theta의 추정량의 반비례 함수 꼴로 이루어져 있으므로, 두 추정량은 음의 상관관계를 갖게 되는 것이다.

위의 그래프에서 우리가 주목해야할 것은 각 추정값들의 퍼짐 정도이다. 두 추정량 모두 실제 패러미터를 중심으로 퍼져있지만 MLE의 퍼짐 정도가 MOME보다 작은 것을 알 수 있다. 따라서 모수를 추정하는 것에 있어서 MLE를 사용하는게 MOME를 사용하는 것보다 좀 더 안정적인 추정을 할 수 있는 것이다.


SHARE TO



티스토리 툴바