Maximum Likelihood and Profile Likelihood for the SEIR model

통계학은 많은 부분 확률모형의 모수를 추정하는 (inferential statistics) 과정이고 모수 추정방법으로 가장 많이 사용되는 방법이 maximum likelihood (ML)이다. 이번 포스트는 2014년 출간된 Cole et al.의 Maximum Likelihood, Profile Likelihood, and Penalized Likelihood: A Primer을 차용하여 maximum likelihood (ML) 와 profile likelihood에 대하여 기술하여 보고자 한다.

최대 가능도 (ML)

잠복기를 추정하기 위해 증상 발현일과 기존 감염자와의 접촉일을 묻는 설문조사를 했다고 하자. \(n\) 명을 인터뷰하고 \(n\) 개의 관찰값 \(y_1, y_2, ..., y_n\) 을 얻었다고 하자. 잠복기의 분포에 대한 확률모형 \(f (y|\boldsymbol{\theta})\) 은 주어진 모수\(\boldsymbol{\theta}=(\theta_1, \theta_2, ..., \theta_j)\) 하에서 \(Y=y\)가 될 확률을 나타낸다.

최대 가능도 방법은 미지의 모수 하에서 관찰값의 확률의 나타낸다. 확률 모형 \(f(y|\boldsymbol{\theta})\)이 주어진 모수 \(\boldsymbol{\theta}\) 하에서 \(Y\)의 확률을 나타내는 반면 최대 가능도 방법은 \(Y\)를 관찰값에 고정한 채 \(\boldsymbol{\theta}\)의 함수로 표현하게 된다. 따라서 확률모형과는 다르게 다음과 같은 식 \(\mathcal{L}(\boldsymbol{\theta};y_i)\) 을 사용한다. 즉, 우리가 관심있어 하는 것은 확률모형 \(f (y|\boldsymbol{\theta})\)이 \(Y\) 가 아니고 \(\boldsymbol{\theta}\) 에 따라 어떻게 변하는가 하는 것이다. \(\mathcal{L}(\boldsymbol{\theta};y_i)\) 를 \(i\) 번째 관찰값이 가능도에 영향을 미치는 정도라 하고 관찰값이 상호독립적이라고 가정하면 관찰값 전체의 가능도는 아래와 같이 표현할 수 있다.

\[\mathcal{L}(\boldsymbol{\theta};\boldsymbol{y}) = \prod_{i=1}^{n} \mathcal{L}(\boldsymbol{\theta};y_i) = \prod_{i=1}^{n} f(y_i;\boldsymbol{\theta})\]

위 식에서 \(\boldsymbol{y}=(y_1, y_2, ..., y_n)\)을 나타낸다.

\(\mathcal{L}(\boldsymbol{\theta};\boldsymbol{y})\) 는 \(\boldsymbol{\theta}\)에 대한 확률을 알 수는 없기 때문에 확률모형이 아닌 가능도 (혹은 우도) 함수라고 한다. ML은 가능도 함수를 최대로 만들어 주는 \(\boldsymbol{\theta}\)로 모수에 대한 추정치를 정의한다.

\[\hat{\theta} = \textrm{argmax}_{\theta}\{{\mathrm{log} \mathcal{L}(\theta)}\}\]

위 식에서 \(\mathrm{log}\)를 사용한 이유는 가능도 값이 매우 작은 수가되는 경우가 많고 따라서 컴퓨터를 이용한 계산상의 안정성을 위해서 (i.e., arithmetic underflow 가 일어나지 않게 하기 위해) 실제로는 \(\mathrm{log} \mathcal{L}(\theta)\)를 사용하기 때문이다. 추가적으로 많은 최적화 알고리듬의 경우 최소화가 기본값으로 설정되어 있어 최대 가능도법을 구현할 때는 \(-\mathrm{log} \mathcal{L}(\theta)\)를 사용하는 경우가 많다.

최대 가능도법을 이용하여 푸아송 분포의 모수를 추정하는 과정을 살펴보자.

푸아송 분포 모수 추정

위에서 언급했던 잠복기의 예를 살펴보자. 잠복기는 Weibull, Gamma, 혹은 Lognormal 등 두 개의 모수를 가지는 확률모형이 많이 사용되는 데 아래 예에서는 계산상의 편의를 위해서 하나의 모수를 가지는 푸아송 분포를 사용하였다.

set.seed(1220)
n <- 50 # number of observations
lamb <- 23 # true parameter value
y <- rpois(n, lambda=lamb) # observations
nll_theta <- function(theta){
  - sum(dpois(y, lambda=theta, log=T)) # negative log likelihood
}
res = optimize(f=nll_theta, interval=c(0,1e6))
res$minimum #\hat{\theta} compare w/ lamb

[1] 24.2

exp(- res$minimum) # likelihood

[1] 3.090828e-11

다음 번 포스팅에는 ML로 추정된 모수의 신뢰구간을 구하는 방법을 샆펴보자.