미분 없이 컴퓨터로 최대가능도 추정하기

지난 포스팅에서 최대가능도법의 개념에 대해 간략히 설명했습니다. 반복하자면, 최대가능도법은 자료가 있을 때 그것을 가장 잘 설명하는 모수치의 값을 찾는 것을 목표로 하고, 구체적으로 가능도함수를 최대화하는 모수치의 값을 구하는 방식으로 이루어진다고 했습니다. 보통 미분을 통해 이 최대화 문제를 해결한다고 설명하는 경우가 많은데, 이 “미분”이라는 말이 주는 압박감 때문에 더이상의 이해를 포기하신 분들도 더러 있으리라 생각합니다. 이 글에서는 그래서 미분 없이 최대가능도법으로 모수치를 추정하는 방법을 간단한 예시를 통해 보여드리려 합니다. 여기서 컴퓨터의 도움을 좀 받아야 하는데, 실제로도 최대가능도법은 컴퓨터의 도움을 받는 경우가 많습니다.

대부분의 독자들에게 가장 친숙한 예시는 역시 정규분포일 것입니다. 정규분포에서 n=100의 자료를 추출하여 모평균을 추정하는 문제를 생각해 봅시다. 여기서 모분산은 1임을 안다고 가정하고, X_i ~ N(5, 1^2) 라고 해 봅시다.는 (5는 실제 모평균의 값이긴 하지만 분석가는 모른다고 칩시다.) R에서는 다음과 같이 할 수 있습니다:

data <- rnorm(100, 5, 1) rnorm은 정규분포에서 랜덤 샘플링을 하는 함수, 100은 샘플 크기, 5는 모평균, 1은 모표준편차입니다. 이제 모분산은 알고 모평균을 모르는 상태에서 최대가능도 추정을 한다고 생각해 봅시다. 그러면 가능도함수는 모평균의 함수로 나타낼 수 있습니다. 원래대로라면 가능도함수는 각각의 자료의 확률밀도함수의 곱으로 표현할 수 있지만, 최대가능도법을 컴퓨터로 구현할 때는 보통 이것을 직접 다루지 않고 로그를 취해 "로그-가능도" log-likelihood로 만듭니다. 그 이유는 확률밀도함수의 곱은 매우 빠르게 아주 작은 값이 되어 버려서 어느 정도 데이터가 커지면 컴퓨터가 허용하는 소숫점 연산 정밀도로는 다룰 수 없을 정도가 되어 버리며, 따라서 연산의 정밀도가 떨어지기 때문입니다. 그리고 어떤 숫자들의 곱에 로그를 취하면 합으로 바꿀 수 있는데, 보통 곱보다는 합을 다루는 게 쉽고 더 정밀하죠. 그래서 보통은 원 가능도함수 대신 로그-가능도함수를 최대화합니다. (로그함수는 증가함수이기 때문에, 원 함수를 최대화하는 것은 그 함수의 로그 값을 최대화하는 것과 같습니다.) 일단 데이터와 모평균을 받아서 가능도함수를 계산해주는 함수를 하나 짜 봅시다. 다음과 같이 하면 됩니다: log_likelihood <- function(mu, data) -sum(log(dnorm(data, mu, 1))) 모분산(표준편차)은 안다고 가정하기 때문에 1로 고정시켰습니다. 그런데 한 가지 눈여겨볼 것은 로그가능도함수 대신 그것에 마이너스 부호를 붙인 것을 함수값으로 지정했다는 것인데, 이는 앞으로 사용할 최적화 명령어는 함수를 최대화 대신 최소화하는 기능만 있기 때문입니다. 하지만 함수의 음의 값을 최소화하는 것은 원래의 값을 최대화하는 것과 같기 때문에, 이렇게 하면 원 로그-가능도함수는 최대화됩니다. 이제 로그-가능도함수를 최대화해 봅시다. R은 미분 없이도 목적함수를 (수치적으로) 최소화하는 명령어들을 제공하는데, 가장 간단한 것은 optim() 함수입니다. 다음 명령어를 실행하면 최대화 결과를 fit 이라는 변수에 저장해 줍니다: fit <- optim(0, log_likelihood, data=data) 0은 알고리즘이 시작하는 값인데 지금 당장은 크게 중요하지 않습니다. (대체로 어떤 값을 넣든 같은 결과를 산출해 줍니다.) 우리는 데이터가 주어진 상태에서 가능도함수를 최대화하는 mu의 값을 찾고 싶은 것이므로, 함수의 입력값에서 data의 값은 data=data로 고정시켜 줍니다. 실행하면 뭐라뭐라 경고 메세지가 나오는데 지금은 별로 안 중요하므로 무시합시다. optim이 찾은, 로그-가능도함수를 최대화하는 mu의 값은 fit의 par라는 변수에 저장되어 있습니다. 저는 다음과 같은 결과를 얻었습니다 (할 때마다 조금씩 달라질 수 있습니다): fit$par [1] 5.1125 최대가능도법에 의해 찾은 mu의 값은 5.1125입니다. 이 값을 최대가능도 추정치 maximum likelihood estimator 라 부릅니다. 그런데 정규분포의 경우, 평균에 대한 최대가능도 추정치의 이론적 값은 표본평균과 일치한다고 우리는 이미 알고 있습니다. 자료의 표본평균을 계산해 보니 다음과 같았습니다: mean(data) [1] 5.112447 수치적 방법을 썼기 때문에 최대가능도 추정치와 약간의 차이는 나지만, 표본평균과 대체로 일치하는 것을 확인할 수 있습니다. 이렇게 가능도함수의 값만 계산할 수 있으면 복잡한 미분 없이도 R 등 통계 소프트웨어가 제공하는 최적화 함수를 이용하면 최대가능도 추정이 가능합니다. 그리고 저도 이 방법을 매우 자주 사용합니다. 해본 김에 이번에는 모표준편차까지도 추정해 봅시다. 로그가능도 함수를 다시 짜야겠습니다: log_likelihood <- function(par, data) -sum(log(dnorm(data, par[1], exp(par[2])))) 함수의 표준편차 자리에 지수함수 - exp() - 를 쓰는 이유는 최적화 과정에서 모수치의 값을 이래저래 바꾸게 되는데, 그 과정에서 분산이 음의 값이 돼버리면 함수값 계산이 안 돼서 (정규분포에서 표준편차는 음의 값이 될 수 없죠) 루틴이 멈추기 때문입니다. (최대가능도법의 좋은 점은 변환을 해서 추정한 값도 역변환을 하면 원래 값이 제대로 추정된다는 것인데, 이를 invariance라 부릅니다.) 이제 이 문제는 모평균 뿐 아니라 모분산도 추정해야 하기 때문에 2차원 최적화 문제입니다. 가능도함수를 다시 최대화해 봅시다: fit <- optim(c(0,0), log_likelihood, data=data) fit$par[1] [1] 5.1125566 exp(fit$par[2]) [1] 1.184361 이 값들을 표본 추정치와 비교해 봅시다: mean(data) [1] 5.112447 sd(data) [1] 1.19038 보니까 자료의 평균은 최대가능도 추정치와 대체로 맞는데, 표본표준편차와 그 최대가능도 추정치는 약간 안 맞습니다. 그런데 기억을 되살려보면, 표본분산은 sum of squares를 (n-1)로 나눈 것이고 분산의 최대가능도 추정치는 n으로 나눈 것입니다. 애초에 분모가 다르니까 값이 차이가 날 수밖에 없습니다. 다시 제대로 계산해보면 sqrt(var(data)*99/100) [1] 1.184413 이 값은 최대가능도 추정치인 1.184361과 거의 비슷합니다. 이렇게 컴퓨터를 잘 활용하면 최대가능도법은 알고보면 그렇게까지 어렵지는 않고, 그냥 가능도함수가 주어졌을 때 그 로그 값을 최대화시켜 주는 모수치 값을 어떤 방법으로든 찾는 과정이라고 생각하시면 편합니다. 그리고 그 과정은 - 이 글에서 묘사한 것과 같이 - 대개 컴퓨터가 알아서 합니다. 때에 따라 기술적 이슈가 생길 때도 있긴 하지만요. 이 문제에 대해서는 나중에 다룰 기회가 있으면 다시 설명하도록 하겠습니다.

About The Author

Related posts

Leave a Reply

Your email address will not be published. Required fields are marked *