Metropolis-Hasting sampling

References

Load packages

library(magrittr)
library(ggplot2)

Examine the target distribution

$f(x) \propto \exp \left(- \frac{(x - \mu)^2}{2\sigma^2}\right)x^2 I(x \ge 0)$

For simplicity, $\mu = 0, \sigma^2 = 1$ are used. The normalizing constant is the inverse of $.

Reference: http://math.stackexchange.com/questions/66084/integral-int-infty-infty-x2-e-x2-mathrm-dx

## Target density constructor (normalized)
construct_target_density <- function(mu, sigma) {
    function(x) {
        ifelse(x >= 0,
               exp(-(x - mu)^2 / (2 * sigma^2)) * x^2,
               0) / (sqrt(2*pi) / 2)
    }
}

## Simple case with mu = 0, sigma = 1
target_density0 <- construct_target_density(mu = 0, sigma = 1)

## Examine visually to examine peak height
g1 <- ggplot(data = data.frame(x = 0), mapping = aes(x = x)) +
    geom_path(stat = "function", fun = target_density0,
              n = 10001, mapping = aes(color = "1")) +
    scale_x_continuous(limits = c(-10,10)) +
    scale_color_manual(name = "Functions",
                       values = c("blue", "red", "green"), # Color specification
                       labels = c("Target", "N(0,1^2)", "10*N(0,3^2)")) +
    theme_bw() + theme(legend.key = element_blank())
print(g1)

Implement Metropolis-Hasting sampler

The proposal distribution proposing a move to $x^{(j+1)}$ from the current value $x^{(j)}$ should be a distribution centered at $x^{(j+1)}$ and should respect the nonnegative support. Thus, we will use $Expo\left(- \frac{1}{x^{(j)}}\right)$ as the proposal distribution. Importantly, this proposal distribution is not symmetric around the current value, therefore, $q(x^{(j+1 )} | x^{(j)}) = q(x^{(j)} | x^{(j+1)})$ does not hold in general. Thus, the acceptance probability is $\alpha(x^{(j+1)} | x^{(j)}) = \min \left(1, \frac{\pi(x^{(j+1)}) q(x^{(j)} | x^{(j+1)})}{ \pi(x^{(j)}) q(x^{(j+1)} | x^{(j)})}\right)$.

MHSampling <- function(target_density, proposal_density, proposal_generator, init, R) {

    ## Initialize a vector
    x <- numeric(length = R)
    ## Acceptance indicator
    a <- logical(length = R)
    ## Initial value
    x[1] <- init

    for (r in seq_len(R)[-1]) {
        ## Obtain a proposed value conditioning on the previous value
        proposal <- proposal_generator(x[r-1])
        ## Evaluate both densities at the proposed value
        target_at_proposal <- target_density(proposal)
        target_at_previous <- target_density(x[r-1])

        ## This is assured to be [0,1]
        acceptance_ratio <- min((target_at_proposal * proposal_density(x[r-1], proposal)) /
                                (target_at_previous * proposal_density(proposal, x[r-1])),
                                1)

        if (runif(n = 1) <= acceptance_ratio) {
            ## If U(0,1) is less than or equal to acceptance ratio
            ## Make a move
            x[r] <- proposal
            a[r] <- TRUE
        } else {
            ## Otherwise, don't move
            x[r] <- x[r-1]
        }
    }

    data.frame(x = x,
               a = a)
}

out1 <- MHSampling(target_density = target_density0,
                   proposal_density = function(x_next, x) {dexp(x = x_next, rate = 1 / x)},
                   proposal_generator = function(x) {rexp(n = 1, rate = 1 / x)},
                   init = 3,
                   R = 10^5)
out1$i <- seq_along(out1$x)
## Acceptance ratio
mean(out1$a)

## [1] 0.40905

## Plot of chain
ggplot(data = out1, mapping = aes(x = i, y = x)) +
    geom_point(alpha = 1 / 5) +
    geom_line(alpha = 1 / 5) +
    theme_bw() + theme(legend.key = element_blank())

## Check auto-correlation (slicing by 10 looks ok)
acf(out1$x)

## Plot of chain (every k-th point)
ggplot(data = subset(out1, i %% 10 == 0), mapping = aes(x = i, y = x)) +
    geom_point(alpha = 1 / 5) +
    geom_line(alpha = 1 / 5) +
    theme_bw() + theme(legend.key = element_blank())

## Density plot after slicing
g1 + geom_density(data = subset(out1, i %% 10 == 0),
                  mapping = aes(x = x))

Metropolis sampler with asymmetric proposal (error)

When the proposal is not symmetric around the current value, the Metropolis algorithm cannot be used. It gives an incorrect sample.

MHSampling <- function(target_density, proposal_density, proposal_generator, init, R) {

    ## Initialize a vector
    x <- numeric(length = R)
    ## Acceptance indicator
    a <- logical(length = R)
    ## Initial value
    x[1] <- init

    for (r in seq_len(R)[-1]) {
        ## Obtain a proposed value conditioning on the previous value
        proposal <- proposal_generator(x[r-1])
        ## Evaluate both densities at the proposed value
        target_at_proposal <- target_density(proposal)
        target_at_previous <- target_density(x[r-1])

        ## This is assured to be [0,1]
        acceptance_ratio <- min((target_at_proposal) / (target_at_previous),
                                1)

        if (runif(n = 1) <= acceptance_ratio) {
            ## If U(0,1) is less than or equal to acceptance ratio
            ## Make a move
            x[r] <- proposal
            a[r] <- TRUE
        } else {
            ## Otherwise, don't move
            x[r] <- x[r-1]
        }
    }

    data.frame(x = x,
               a = a)
}

out1 <- MHSampling(target_density = target_density0,
                   proposal_density = function(x_next, x) {dexp(x = x_next, rate = 1 / x)},
                   proposal_generator = function(x) {rexp(n = 1, rate = 1 / x)},
                   init = 3,
                   R = 10^5)
out1$i <- seq_along(out1$x)
## Acceptance ratio
mean(out1$a)

## [1] 0.49602

## Plot of chain
ggplot(data = out1, mapping = aes(x = i, y = x)) +
    geom_point(alpha = 1 / 5) +
    geom_line(alpha = 1 / 5) +
    theme_bw() + theme(legend.key = element_blank())

## Check auto-correlation (slicing by 10 looks ok)
acf(out1$x)

## Plot of chain (every k-th point)
ggplot(data = subset(out1, i %% 10 == 0), mapping = aes(x = i, y = x)) +
    geom_point(alpha = 1 / 5) +
    geom_line(alpha = 1 / 5) +
    theme_bw() + theme(legend.key = element_blank())

## Density plot after slicing
g1 + geom_density(data = subset(out1, i %% 10 == 0),
                  mapping = aes(x = x))

• Top Page: http://rpubs.com/kaz_yos/
• github: https://github.com/kaz-yos