내 맘대로 만드는 MCMC 간단 예제

Summary

MCMC(Markov chain Monte Carlo) 는 샘플링 방법이다. 여기서는 계산가능한 함수 f(x) 를 타겟분포로 주어졌을 때 이 함수로부터 샘플링하는 간단한 예제를 작성해본다. 자세한 내용은 아래 참고 사이트를 참조한다.

참고 사이트:

# 타겟분포 설정
f_x <- function (x) {
  0.3 * exp(-0.2*x^2) + 0.7*exp(-0.2*(x-10)^2) # 타겟 분포
} 

# 타겟 분포를 그려본다. 그래프 이하의 면적은 당연히 1 이 아니기 때문에 PMF 가 아니다.
x <- seq(-10, 20, by=0.01)
plot(x, f_x(x), type='l', col='red')

# 제안분포 설정
g_x <- function (x) {
  # 제안분포로 정규분포를 사용한다. 표준편차는 임의로 설정한다.
  rnorm(1, mean = x, sd = 2)
}

# 반복 횟수 설정
N <- 10000

# x 의 현재 상태를 초기화한다. 랜덤하게 생성해도 좋고, 고정으로 해도 상관없다.
x_current <- 0

result <- c() # 샘플링된 결과 수집

for (i in 1:N) {
  # 0과 1사이의 균등분포 샘플링을 한번 수행해서, 제안된 값을 받아들일지 거부할지 결정한다.
  u <- runif(1) 
  
  # 제안분포에서 샘플링을 한번 수행한다. 제안분포: 정규분포
  x_next <- g_x(x_current)
  
  # 현재 상태와 다음 상태에서의 함수값을 계산하고 그 비율을 계산한다.
  # 그 중 가장 작은 값을 선택한다.
  a_ij <- min(1, f_x(x_next) / f_x(x_current))
  # 제안분포가 대칭이기 때문에 아래에서 dnorm 에 해당하는 부분은 약분되고 결국 위의 식과 동일해진다.
  #a_ij <- min(1, (f_x(x_next) * dnorm(x_current, x_next, sd=2)) / (f_x(x_current) * dnorm(x_next, x_current, sd=2)) )
  
  
  # 제안된 값 x_next 를 accept 할지 reject 할지 결정한다.
  # 위에서 a_ij 는 0보다 크고 1보다 같거나 작은 값을 가진다. 
  # 만약 1 이라면, 무조건 제안된 값을 선택하고
  # 1 이하라면, 선택할지 말지를 확률적으로 결정한다.
  if (u < a_ij) {
    # 만일 a_ij 가 크다면, 다음 상태값이 현재 상태로 변경된다.
    x_current <- x_next
    
  } else {
    # 만일 다음 상태에서의 함수값이 현재 상태보다도 작고, 균등분포로 구한 랜덤값보다도 작으면 
    # 현재 상태를 유지한다.
    
  }
  
  result <- c(result, x_current)
}
hist(result, freq = FALSE, breaks = 100, ylim = c(0, 1))
lines(x, f_x(x), col='red')

# 함수값 기대치 계산
mean(f_x(x))

## [1] 0.1320669

# 위에서 얻은 샘플링 데이터를 이용해서, 테스트로 복원 샘플링을 해본다.
hist(sample(result, 1000, replace = T), breaks = 100)

내 맘대로 만드는 MCMC 간단 예제

Seokkyu Kong

2021-10-12

Summary