Sampling-Importance-Resampling algorithm

Sampling-Importance-Resampling algorithm

어제에 이어 Sampling-Importance-Resampling algorithm 코드를 올려본다. 물론 책의 코드를 참고했는데, R코드를 구현하다 보니 거의 pseudo코드에 근접하는 것을 알 수 있었다. 이는 시뮬레이션과 데이터에 최적화된 R의 특징 때문이라 생각된다.

특히나 가중치 값에 따라 resampling하는 부분이 함수 하나로 치환되었고 컴퓨팅 속도도 무척 빨라졌다.

weight값에 따라 재 샘플링하는 부분이 어떻게 동작하는지 comment된 부분을 참고하는 것도 좋을거 같고, 이것 이외에 다른 가중샘플링 기법을 생각해 보는 것도 좋을듯 하다.

set.seed(1000)
# The Sampling-Importance-Resampling algorithm

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

q <- function(x) {
    4
}

sir <- function(n) {
    w <- vector(n, mode = "numeric")
    # sample2 <- vector(n, mode='numeric')

    # sample from q
    sample1 <- runif(n, 0, 4)

    # 특정 점에서 나올 가중치를계산
    w <- p(sample1)/q(sample1)
    # 가중치 정규화 값
    w <- w/sum(w)

    # 가중치에 따른 재 샘플링 cumw <<- cumsum(w)

    # 정규화 되었으니 0,4까지 구할 필요가 없다.  u <- runif(n) for(i in
    # seq(n)){ indices <- which(u < cumw[i]) sample2[indices] <- sample1[i]
    # u[indices] <- 100 }
    sample2 <- sample(sample1, n, prob = w, replace = T)
    return(sample2)
}

시각화

rejection sampling과 같은 시각화 결과를 보여주는 것을 볼 수 있다.

library(ggplot2)

# real data

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

ggplot(sir_sampling, aes(sir)) + 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))

## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust
## this.