Understanding MCMC - 1
MCMC의 기본 아이디어를 이해 해 보자.
난수 생성을 통해 원하는 정보를 알아내는 Monte Carlo method, 현재의 상태는 바로 전 상태에만 dependent 하다는 가정의 Markov Chain, 이 둘을 결합한 것이 MCMC 라고 볼 수 있다.
모델의 모수 추정 등의 상황에서 analytic method로 solution을 도출하기 어려울 경우 MCMC를 고려하게 된다. 특히 Bayesian method로 모델링을 할 경우 MCMC를 통해 모수의 candidates를 sampling하고 확률분포를 알아내는 것이 일반적이다.
본 포스팅에서는 2개의 설명변수와 intercept를 가지는 linear regression model 의 coefficient를 찾는 과정을 MCMC 컨셉으로 구현 해 보았다.
MCMC 알고리즘에는 Metropolis, Metropolis Hastings, Gibbs sampling, Hamiltonian MC 등이 있으며 본 포스팅에서는 Metropolis 알고리즘을 참고하여 실험 해 보도록 한다. 알고리즘을 엄밀하게 구현 한 것은 아니며 작동 과정을 나름의 방식으로 이해 해 보기 위해 약식으로 구현 해 본 것임을 먼저 밝힌다.
알고리즘을 구체적으로 구현해 보기 앞서서 MCMC의 기본 아이디어를 이용해서 문제의 답을 찾는 과정을 먼저 실험 해 보자.
아래와 같이 테스트를 위한 데이터를 생성하였다.
vars <- 2
obs <- 100
x <- replicate(vars, rnorm(obs))
X <- cbind(1, x)
coefs <- c(3,5,7)
mu <- X %*% coefs
sigma <- 5
y <- rnorm(obs, mu, sigma)이렇게 만들어진 데이터로 기본적인 linear model fitting을 해 보면 다음과 같다.
mod.lm <- lm(y~., data=data.frame(X[,-1]))
summary(mod.lm)##
## Call:
## lm(formula = y ~ ., data = data.frame(X[, -1]))
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.4075 -2.1986 0.5259 2.6423 9.5323
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.0374 0.4346 6.989 3.51e-10 ***
## X1 5.0136 0.3946 12.704 < 2e-16 ***
## X2 7.5366 0.4596 16.399 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.342 on 97 degrees of freedom
## Multiple R-squared: 0.8268, Adjusted R-squared: 0.8233
## F-statistic: 231.6 on 2 and 97 DF, p-value: < 2.2e-16이제 MCMC 컨셉을 활용해서 coefficient를 찾는 과정을 간단하게 구현 해 보자.
아래 작성한 코드는 대략 다음과 같이 동작한다.
임의의 coefficient를 생성하여 시작.
임의 샘플링한 coefficient로 모델의 log likelihood를 계산하고 이전 iteration의 log likelihood와 비교해서 개선이 있었으면 해당 sample 채택, 그렇지 않으면 버리고 이전 sample을 유지.
2에서 채택 혹은 유지 된 sample을 기준으로 proposal distribution을 활용하여 새로운 sample을 생성.
2,3 의 과정을 충분히 반복.
init_beta <- rep(0.5,vars+1) # 임의의 시작점 : 0.5
init_yhat <- X %*% init_beta
init_error <- init_yhat - y
ll.old <- sum(dnorm(init_error, mean=0, sd=5,log=T)) # 시작점에서의 log likelihood
iter <- 2000
betaset <- matrix(0, nrow=iter, ncol=vars+1)
betaset[1,] <- init_beta
acceptance <- rep(FALSE, iter)
sigma.hat <- diag(1,nrow=3) # proposal distribution 으로 mvnorm을 사용했으며 이를 위한 cov-matrix
for(i in c(2:iter)){
beta <- MASS::mvrnorm(1, mu=betaset[i-1,], Sigma = sigma.hat)
yhat <- X %*% beta
error <- yhat - y
ll.new <- sum(dnorm(error, mean=0, sd=5,log=T))
success <- round(ll.new,1) >= round(ll.old, 1)
acceptance[i] <- success
if (success==TRUE) {
betaset[i,] <- beta
ll.old <- ll.new
}
else {
betaset[i,] <- betaset[i-1,]
ll.old <- ll.old
}
if (sum(acceptance) > 5){
sigma.hat <- cov(betaset[acceptance,])
}
else {sigma.hat <- sigma.hat}
}구현 방식과 관련해서 생각해 볼 문제들이 있는데 아래 두 가지다.
이전 단계의 sample 로부터 현 단계의 sample의 improvement를 비교하는 과정에 log likelihood를 써서 코드를 작성 하였지만 필수 불가결한 것은 아니다. 더 간단히 SSE 등을 써도 목적 달성에는 문제가 없다.
proposal distribution으로 gaussian 분포를 사용했을 경우 covariance matrix가 성능에 중요한 역할을 한다. solution 탐색 범위를 지정하는 작용을 하기 때문이다. 범위가 너무 좁거나, 반대로 너무 넓게 지정될 경우 실제 solution에 수렴하지 못하거나 대단히 많은 iteration이 소요 될 수도 있을 것임을 예상할 수 있다. 위 코드에서는 임의의 covariance matrix를 시작에 사용하되, 적합한 sample이 어느정도 쌓이면 adaptive하게 갱신되도록 장치를 넣어 보았다.
이제 sampling 된 coefficient를 trace 시각화 해 보자.
par(mfrow=c(1,3))
plot(betaset[,1], type='l', xlab = 'iteration')
plot(betaset[,2], type='l', xlab = 'iteration')
plot(betaset[,3], type='l', xlab = 'iteration')
다행스럽게도 각 계수들이 정답을 찾아 수렴 하였다. (복잡도가 낮은 2변수 linear regression model로 수렴이 비교적 수월한 케이스로 볼 수 있다.)
1999번의 iteration 중에 실제로 이전 단계보다 개선을 이룬 iteration은 생각보다 많지 않았다. 빠르게 수렴했고 계수의 변동폭이 크지 않다는 얘기다. 이 역시 모델의 복잡도가 높지 않은데다가, 데이터에 포함된 노이즈가 적었기 때문으로 보인다.
sum(acceptance)## [1] 26이번엔 coefficient 2개씩 paring 해서 solution을 찾아온 경로를 살펴 보자.
par(mfrow=c(1,3))
plot(betaset[,c(1,2)], type="l")
plot(betaset[,c(2,3)], type="l")
plot(betaset[,c(1,3)], type="l")
데이터의 likelihood를 개선시키는 방향을 잡아 난수를 생성(sampling)하고 accept/reject을 반복하면서 점진적으로 model의 coefficient를 찾아가는 과정을 확인해 볼 수 있었다.