rejection sampling

겨울학교에서 Gibbs Sampling 방법론에 대한 내용이 나와서 이게 바이오정보학과목에서 본거 같은데 확실한 의미를 몰라서 찾아보다가 여러 샘플링 기법에 대해서 살펴보고 있다.

대부분 복잡한 확률 모델들은 확률값을 추론하기가 힘들다. 이 때문에 여러 샘플링 방법을 사용해서 확률값을 근사하는게 되는데, 책을 보다가 책에 있는 코드를 기반 R로 시뮬레이션을 해봤다.

함수가 있는데 이런식으로 시뮬레이션 샘플링을 하는 이유는 많은 함수가 적분이 힘들기 때문이다. 적분만 가능하면 복잡한 모델이더라도 확률값을 알 수 있을텐데, 현실의 많은 경우 적분이 불가능하기 때문이다.

set.seed(1000)
# rejection sampling algorithm

qsample <- function() {
    runif(1, 0, 4)
}

p <- function(x) {
    0.3 * exp(-(x - 0.3)^2) + 0.7 * exp(-(x - 2)^2/0.3)
}


rejection <- function(nsamples) {
    M <- 0.72
    samples <- vector(nsamples, mode = "double")
    count <- 0
    for (i in seq(nsamples)) {
        accept <- FALSE
        while (accept != TRUE) {
            x <- qsample()
            u <- runif(1, 0, M)
            if (u < p(x)) {
                accept <- TRUE
                samples[i] <- x
            } else {
                count <- count + 1
            }
        }
    }
    print(count)
    return(samples)
}

시각화

결과물 시각화인데, 10000개의 샘플을 뽑기 위해 약 18000개의 샘플이 rejection되었다(물론 이 횟수는 매번 달라진다).

library(ggplot2)

# real data

real <- data.frame(real_y = p(seq(-1, 5, by = 1e-04)), real_x = seq(-1, 5, by = 1e-04))
rej_sampling <- data.frame(rejection = rejection(10000))

## [1] 17354

ggplot(rej_sampling, aes(rejection)) + geom_histogram(aes(y = ..density.., fill = ..density..)) + 
    geom_line(aes(x = real_x, y = real_y), data = real, colour = "red2") + scale_fill_gradient2(guide = guide_legend(reverse = TRUE))

후기

요즘 확률관련된 샘플링이라든지 시뮬레이션에 관심이 많이 가서 관련 자료들을 보고 있는데, 재미있다. 게다가.. 컴퓨팅 연산도 꽤 많은 작업이라서 더더욱 맘에 들고….

Rejection Sampling 시뮬레이션 by from __future__ import dream is licensed under a Creative Commons Attribution-NonCommercial 4.0 International License.