MCMC - Beta 분포
Kong, Seokkyu
2021-10-16
Summary
본 내용은 Introduction to Probability-2nd, Blitzstein의 539~540 페이지 내용을 정리했습니다.
\(W \sim Beta(a, b)\) 를 따르는 랜덤 샘플링을 MCMC(Markov Chain Monte Carlo)를 이용해서 구현해보자. R 에서는 rbeta 를 사용하면 쉽게 구할 수 있는데, 여기서 목표는 이 함수를 흉내 내는 것이다.
먼저, Beta 분포에 대해서 균등분포의 보편성(universality of the Uniform)을 직접 적용하는 것은 어렵다. 왜냐하면 Beta CDF의 역함수를 구하기 어렵기 때문이다. 그래서 만약 \(X \sim Gamma(a, 1)\) 이고 \(Y \sim Gamma(b, 1)\) 이 서로 독립이라면, \(\frac{X}{X+Y} \sim Beta(a, b)\) 를 따르게 된다. 그래서 Gamma 분포를 시뮬레이션 할 수 있다면, Beta 분포를 시뮬레이션 할 수 있다.
\(X \sim Gamma(a, 1)\) 을 시뮬레이션 하기 위해서는 \(X_1 + X_2 + ... + X_a\) 를 사용할 수 있는데 여기서 \(X_j \sim i.i.d. Expo(1)\) 이다. 마찬가지로 \(Y_j \sim i.i.d. Expo(1)\) 이다. Expo(1) CDF 의 역함수를 사용하면, \(-log(1-U) \sim Expo(1)\) 이 된다.
MCMC를 이용한 Beta분포 랜덤 샘플링
\(W \sim Beta(a, b), a>0, b>0\) 일때, Beta 분포의 PDF 는 \(f(x) = c * x^{a-1} * (1-x)^{b-1}, 0<x<1\) 이고 여기서 c 는 정규화 상수이다. [2]
결과로 생기는 Metropolis-Hastings chain을 Independence sampler 라고 부른다.
set.seed(2021)
# 타겟분포 설정
a <- 5
b <- 2
f_x <- function (x, a, b) {
# 정규화 상수 c 는 어차피 약분되기 때문에 사용하지 않아도 되지만,
# 여기서는 그래프를 그리기 위해서 임의로 설정했다.
10 * x^(a-1) * (1-x)^(b-1)
}
# 타겟 분포를 그려본다. 그래프 이하의 면적은 당연히 1 이 아니기 때문에 PMF 가 아니다.
x <- seq(0, 1, by=0.01)
plot(x, f_x(x, a, b), type='l', col='red')
# 반복 횟수 설정
N <- 10000
# x 의 현재 상태를 초기화한다. 랜덤하게 생성해도 좋고, 고정으로 해도 상관없다.
x_current <- runif(1)
result <- c() # 샘플링된 결과 수집
for (i in 1:N) {
# 제안분포를 통해서 0과 1사이의 값을 1개 제안을 받는다.
u <- runif(1)
# 현재 상태와 다음 상태에서의 함수값을 계산하고 그 비율을 계산한다.
# 그 중 가장 작은 값을 선택한다.
a_ij <- min(1, f_x(u, a, b) / f_x(x_current, a, b))
# 제안된 값 u 를 accept 할지 reject 할지 결정한다.
# 위에서 a_ij 는 0보다 크고 1보다 같거나 작은 값을 가진다.
# 만약 1 이라면, 무조건 제안된 값을 선택하고
# 1 이하라면, 선택할지 말지를 확률적으로 결정한다.
if (runif(1) < a_ij) {
# 만일 a_ij 가 크다면, 다음 상태값이 현재 상태로 변경된다.
x_current <- u
} else {
# 만일 다음 상태에서의 함수값이 현재 상태보다도 작고, 균등분포로 구한 랜덤값보다도 작으면
# 현재 상태를 유지한다.
}
result <- c(result, x_current)
}
hist(result, freq = FALSE, breaks = 100)
lines(x, f_x(x, a, b), col='red')
# rbeta 를 이용해서 결과를 비교한다.
hist(rbeta(N, a, b), freq = FALSE, breaks = 100)
References
[1] Joseph K. Blitzstein and Jessica Hwang, 2019, Introduction to Probability 2nd, CRC Press, 539-540