MCMC-Zipf 분포 시뮬레이션

Summary

본 내용은 Introduction to Probability-2nd, Blitzstein의 538-539 페이지 내용을 정리했습니다.

\(M \ge 2\) 이며 정수라고 하자. 확률 변수 X 는 파라미터 \(a \gt 0\) 를 갖는 Zipf 분포라고 하면, PMF는 아래와 같다.

\(P(X = k) = {1/k^a \over \sum_{j=1}^{M}(1/j^a)}\) , for k=1,2,…,M

이 분포는 단어 빈도를 연구하는 언어학에서 폭넓게 사용된다. 여기서는 MCMC(Markov Chain Monte Carlo) 를 이용한 Zipf 분포를 시뮬레이션 한다.

제안분포는 {1, 2, …, M} 상에서의 random walk 를 따르게 한다. 상태 i 가 \(i \ne 1\), \(i \ne M\) 일 때 확률 1/2 로 상태 i - 1 또는 i + 1로 이동한다. 상태 1에서는 머물러 있거나 상태 2로 이동한다. 상태 M 에서는 머물러 있거나 M - 1 상태로 이동한다.

set.seed(2021)

#install.packages("VGAM")
library(VGAM)

M <- 10 # 상태 갯수, ex) M: 1~10
a <- 2 # a > 0, a 가 1보다 작으면 완만한 그래프 형성됨.

# 타겟분포 설정
f_x <- function (x) {
  j <- 1:M
  1/(x^a) / sum(1 / j^a)
} 

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

# 반복 횟수 설정
N <- 10000

# 제안 분포 설정
g_x <- function (x) {
  # random walk, 확률 1/2을 가지고 실행한다.
  move <- sample(c(-1, 1), 1)

  x_next <- 0  
  # 양 끝의 상태인 1, M 에서는 예외처리 한다.
  if (x == 1 && move == -1) {
    x_next <- x
  } else if (x == M && move == 1) {
    x_next <- x
  } else {
    x_next <- x + move
  }
  
  return (x_next)
}

# 테스트
print(g_x(20))

## [1] 19

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

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

for (i in 1:N) {
  # 0과 1사이의 균등분포 샘플링을 한번 수행해서, 제안된 값을 받아들일지 거부할지 결정한다.
  u <- runif(1) 
  
  # 제안분포에서 샘플링을 한번 수행한다. 제안분포: random walk
  x_next <- g_x(x_current)
  
  # 현재 상태와 다음 상태에서의 함수값을 계산하고 그 비율을 계산한다.
  # 그 중 가장 작은 값을 선택한다.
  a_ij <- min(1, (x_current)^a / (x_next)^a)
  
  # 제안된 값 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)
lines(x, f_x(x), col='red')

# rbeta 를 이용해서 결과를 비교한다.
hist(rzipf(N, M, a), freq = FALSE, breaks = 100)

MCMC-Zipf 분포 시뮬레이션

Seokkyu Kong

2021-10-22

Summary

References