본문으로 바로가기



감마(Gamma) 분포 최대우도추정 with R - 1. 감마 분포의 pdf에 관하여

category Statistics/Inference 2017.06.27 01:07
감마(Gamma) 분포 최대우도추정 with R - 1. 감마 분포의 pdf에 관하여

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

library(tidyverse)
library(RColorBrewer)

감마 분포(Gamma distribution)에 대하여

감마분포는 통계시간에 한번쯤 들어는 봤을 법한 분포이면서, 약간은 만만치 않은 분포일 것이다(하지만 가장 만만한 지수 분포에서 비롯되었다.). 오늘은 이 감마 분포의 모수를 최대우도 추정법을 사용해서 구해보는 것의 첫 걸음으로 pdf에 대하여 알아보도록 하겠다.

먼저 감마분포의 pdf는 두가지 버젼이 있는데, 이 글에서 사용할 확률밀도함수(pdf)의 형태는 다음과 같다.

\[ f(x;k,\theta) = {\frac {1}{\Gamma (k)\theta ^{k}}}x^{k\,-\,1}e^{-{\frac {x}{\theta }}} \] 위의 식에서 \(\Gamma (k)\)는 감마 함수(gamma function, gamma dist.와는 다르다.)를 의미하고, \(k\)\(\theta\)는 감마 분포의 형태를 결정하는 shape 모수와 범위를 결정하는 scale모수를 의미한다.

위의 식을 R에서 정의하고 Gamma분포의 pdf로 사용해도 되지만, 통계 계산에 최적화된 프로그래밍 언어임을 자처하는 R에는 이미 정의된 dgamma()라는 함수가 있으므로 이것을 이용하도록 하자.

감마분포의 shape 모수

먼저 shape 모수의 의미를 알아보기 위하여, scale 모수를 고정시켜놓고, shape 모수의 값을 변화시켜가면서 pdf의 움직임을 관찰해보자.

my_color <- brewer.pal(5, "Accent")
shape_pr <- seq(2,10, by = 2)
p <- ggplot(data = data.frame(x = 0), mapping = aes(x = x)) +
         geom_hline(yintercept=0) + xlim(0, 20) + ylim(0, 0.5)
for (i in shape_pr){
  p <- p + stat_function(fun = dgamma,
                        args = list(shape = i, scale = 1.0),
                        color= my_color[i/2],
                        size = 1.0)
}
p  + ggtitle("pdf of Gamma dist. w/ various shapes when scale = 1.0") +
     annotate("text",
              x = (shape_pr - 1) + 0.5,
              y = map2_dbl((shape_pr-1), shape_pr, dgamma, scale = 1.0) + 0.02,
              label = paste("shape =", shape_pr))


위의 결과를 보면 shape 모수가 커질 수록 pdf의 형태가 종모양에 가까워지는 것을 알 수 있다. 만약 수리통계를 들었던 적이 있다면, 감마분포를 따르는 확률변수는 여러 지수분포를 따르는 확률변수의 합으로 볼 수 있다는 것을 떠올려보자. 감마분포의 shape 모수가 나타내는 것은 이러한 지수확률 변수가 몇개가 더해졌냐를 나타내는 모수이기도 한데, 그렇다면 shape모수가 크다는 의미는 여러개의 지수 확률변수가 더해졌다는 의미이므로, 중심극한 정리에 의하여 pdf는 종모양에 가까워져야 하는 것이다.

위의 코드에서 annotate() 함수는 그래프 위에 자유롭게 글씨를 첨가할 수 있는 함수인데, 감마 분포의 mode(최빈값)가 \((k-1)\theta\)라는 것을 이용하였다.

감마분포의 scale 모수

이번에는 반대로 shape 모수를 고정시켜놓고, scale 모수의 값을 변화시켜가면서 pdf의 움직임을 관찰해보자.

scale_pr <- seq(1,9, by = 2)
p <- ggplot(data = data.frame(x = 0), mapping = aes(x = x)) +
         geom_hline(yintercept=0) + xlim(0, 20) + ylim(0, 0.5)
for (i in scale_pr){
  p <- p + stat_function(fun = dgamma,
                        args = list(shape = 2.0, scale = i),
                        color= my_color[(i+1)/2],
                        size = 1.0)
}
p  + ggtitle("pdf of Gamma dist. w/ various scale when shape = 2.0") +
     annotate("text",
              x = scale_pr + 0.5,
              y = map2_dbl(scale_pr, 1/scale_pr, dgamma, shape = 2) + 0.02,
              label = paste("scale =", scale_pr))

위의 결과에서 우리는 shape모수가 고정이 되어있을 경우 scale 모수가 커질 수록 pdf의 도메인, 즉 함수값이 0보다 큰 지역의 크기가 넓어지는 것을 확인 할 수 있다. 이것은 앞에서 언급한 shape 모수가 더해진 지수확률 변수의 갯수를 나타내는 모수라면, scale 모수는 각 더해지는 지수분포의 scale이 얼마나 큰가를 나타내는 모수이기 때문이다. (scale 모수 \(\theta\)인 지수 확률변수의 평균은 \(\theta\)가 된다.)

따라서 우리가 그린 예제의 경우, 두 개의 지수분포를 더한다는 가정하에, 각 지수분포의 퍼짐 정도(scale)가 커짐에 따라서 두 확률변수의 합으로 나타내어지는 감마분포의 pdf 도메인 영역이 넓어진다는 것을 확인 할 수 있다.

다음 시간에는 본격적으로 감마 분포의 모수에 대한 MOME(적률 추정량)과 MLE(최대우도 추정량)를 구하는 방법에 대하여 알아보도록 하자.


SHARE TO



티스토리 툴바