```

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())