Understanding MCMC

Metropolis 알고리즘을 구현해 보자

앞서,

임의 sample 생성 -> 모델의 적합도 개선여부 체크하여 sample 채택 -> 채택된 sample을 기준으로 proposal distribution을 통해 새로운 sample 생성

하는 과정을 반복하여 model의 coefficient를 찾아가는 과정을 실험해 보았다.

이번엔 Metropolis 알고리즘을 구체적으로 살펴본다.

앞 포스팅에서 작성한 코드의 몇 부분을 수정해서 Metropolis 알고리즘을 구현해 보자. 이해를 위한, 약식의 구현임을 감안하자.

먼저 이전 포스팅에서와 동일하게 실험을 위한 데이터를 생성하였다.

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)

Metropolis 알고리즘을 구현하는 코드는 다음과 같다.

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))

# prior 반영
lpriors<-function(beta){
  lp1<-dnorm(beta[1], mean=3, sd=1, log=T)
  lp2<-dnorm(beta[2], mean=5, sd=1, log=T)
  lp3<-dnorm(beta[3], mean=7, sd=1, log=T)
  return(lp1+lp2+lp3)
}

postp.old <- ll.old + lpriors(init_beta)   # 초기 likelihood와 prior를 결합하여 초기 posterior 계산

iter <- 1000   #  짧게 1000 iteration 으로 trial       
betaset <- matrix(0, nrow=iter, ncol=vars+1) 
betaset[1,] <- init_beta
acceptance <- rep(FALSE, iter)
sigma.hat <- vcov(lm(y~., data=data.frame(X[,-1]))) 

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))
  postp.new <- ll.new + lpriors(beta)
  success <- (runif(1) < min(1,exp(postp.new-postp.old)))
  
  acceptance[i] <- success
  
  if (success==TRUE) {
    betaset[i,] <- beta
    postp.old <- postp.new
  }
  else {
    betaset[i,] <- betaset[i-1,]
    postp.old <- postp.old
  }
  
}

이전 포스트에서의 코드와 달라진 점은 두 가지다.

coefficient 별 prior를 설정하고 likelihood와 결합하여 posterior를 구성.
이전 sample의 posterior와 현 단계 sample의 posterior를 비교하여 acceptance를 결정하되, 두 posterior probability의 ratio를 확률로 하여 현 단계 sample 채택.

posterior를 개선 시키지 못한 sample 일지라도 낮은 확률로 채택한다는 의미다. 이렇게 sampling 한 결과는 이어지는 결과에서 확인 해 보자.

먼저 trace plot.

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')

burnin period를 거쳐 본래 값으로 수렴하였음을 알 수 있는데, 본래 값 주변의 값들도 드문 빈도로 함께 채택 되었음을 알 수 있다.

다음은 sampling 된 coefficient의 histogram 이다. burnin period를 초기 200개로 보고, 이를 제외하고 그려 보자.

par(mfrow=c(1,3))
hist(betaset[-c(1:200),1], 50)
hist(betaset[-c(1:200),2], 50)
hist(betaset[-c(1:200),3], 50)

coefficient의 확률분포 형상을 알 수 있다. prior를 normal 분포를 사용하여 posterior 역시 normal distribution의 모습을 띈다. (conjugate)

이와 같이 MCMC를 활용하여 모수의 확률 분포를 확인 할 수 있다. 복잡한 모델의 모수를 알고자 할 때, Bayesian method 로 모델을 추정하고자 할 때 빛을 발하게 된다.

대략의 과정과 원리를 알았다. 이번엔 sampling 갯수를 늘려보자(1000 -> 5000).

sampling 갯수를 늘렸더니 확률 분포의 형상이 한 층 뚜렷해 졌다.

갑툭튀 Prior?

갑자기 등장한 Prior가 궁금해진다. Prior는 Bayes theorem에 따라 posterior를 구성함에 있어서 추정하고자 하는 파라메터에 관한 사전적 정보를 반영하게 하는 역할을 한다. Prior에 대한 정보를 모른다면, 혹은 배제하고자 한다면 uniform 분포와 같이 noninformative prior를 사용하던가, simulation 과정에서 prior 부분을 삭제한다. Prior를 배제한다면 모델은 주어진 데이터, likelihood 정보만으로 추정 되는 셈이다.

아래 그림은 이전 실험에서 prior setting을 배제하고 sampling 한 결과다.

현재의 실험 데이터와 모델로는 큰 차이가 나지 않지만, 모델의 복잡도가 높거나 데이터가 충분치 않은 경우 prior의 설정과 활용은 매우 중요한 역할을 하게 된다.

Understanding MCMC - 2

Metropolis 알고리즘을 구현해 보자

갑툭튀 Prior?