Libraries to load

1 Compute chance of a coin flip

1.1 Simulate some data

set.seed(5003)
theta = 0.7  # P(Head)
N = 200
# simulate some coin flips
coin.flips <- rbinom(n = N, size = 1, prob = theta)

1.2 Determine how to compute the Likelihood of the parameter

likelihood <- function(theta, x) prod(dbinom(x, size = 1, prob = theta, log = TRUE))

1.3 Naive Metropolis Hastings

# In this example, we set the prior of P(Head) to be a uniform distribution
niters <- 1000
curr.theta <- 0.5  # This is the parameter we want to find; initialise it 

theta.chain <- numeric(niters)  # Initialise the vector to be created
for (i in 0:(niters - 1)) {
    proposal.theta <- curr.theta + rnorm(1, mean = 0, sd = 0.05)
    
    if (proposal.theta < 0 || proposal.theta > 1) {
        theta.chain[i + 1] <- curr.theta
        next
    }
    
    curr.likeli <- likelihood(curr.theta, coin.flips)
    proposal.likeli <- likelihood(proposal.theta, coin.flips)
    
    alpha <- min(proposal.likeli/curr.likeli, 1)
    
    # Randomly decide whether to accept our new proposal
    if (runif(1) < alpha) 
        curr.theta <- proposal.theta
    theta.chain[i + 1] <- curr.theta
}
plot(theta.chain)

1.4 Log Likelihood computation

log.likelihood <- function(theta, x) sum(dbinom(x, size = 1, prob = theta, log = TRUE))

1.5 Log likelihood computed approach to the Metropolis-Hasting algorithm

# In this example, we set the prior of P(Head) to be a uniform distribution
niters <- 1000
curr.theta <- 0.5  # This is the parameter we want to find; initialise it 

theta.chain <- numeric(niters)
for (i in 1:niters) {
    proposal.theta <- curr.theta + rnorm(1, mean = 0, sd = 0.05)
    
    # If proposed theta leaves the feasible range of values
    if (proposal.theta > 1 || proposal.theta < 0) {
        theta.chain[i + 1] <- curr.theta
        next
    }
    
    # It is easier to work with log of probabilities
    curr.log.likeli <- log.likelihood(curr.theta, coin.flips)
    proposal.log.likeli <- log.likelihood(proposal.theta, coin.flips)
    
    # Because we have taken log(prob), now we need a slightly different formula to
    # get the acceptance probability. Follows from exp(log(a/b)) = exp(log(a) -
    # log(b))
    alpha <- min(exp(proposal.log.likeli - curr.log.likeli), 1)
    
    # Randomly decide whether to accept our new proposal
    if (runif(1) < alpha) 
        curr.theta <- proposal.theta
    theta.chain[i + 1] <- curr.theta
}

plot(theta.chain, type = "l")

hist(theta.chain[100:1000], n = 30)

print(mean(theta.chain[100:1000]))
## [1] 0.6522441
print(sd(theta.chain[100:1000]))
## [1] 0.03574393

2 Generate Samples from a non-trivial distribution

p <- 0.4
mn <- c(-1, 2)
std <- c(0.2, 1.5)
f <- function(x, mu = mn, sd = std) p * dnorm(x, mu[1], sd[1]) + (1 - p) * dnorm(x, 
    mu[2], sd[2])
curve(f(x), col = "blue", from = -4, to = 8, n = 300, las = 1)

2.1 Use Metropolis-Hastings algorithm to draw samples from the normal mixture above

q <- function(x, sigma) rnorm(1, mean = x, sd = sigma)

2.2 Function to do a random walk Metropolis-Hasting algorithm

mc.algorithm <- function(x_0, n, sigm) {
    x <- numeric(n)
    x[1] <- x_0
    for (i in 1:(n - 1)) {
        prop.x <- q(x[i], sigma = sigm)
        alpha <- min(f(prop.x)/f(x[i]), 1)
        if (runif(1) < alpha) 
            x[i + 1] <- prop.x else x[i + 1] <- x[i]
    }
    x
}

2.3 Explore the impact of choices of x_0 and sigma

set.seed(5003)
x_0.vals <- c(-1, 1, 5)
sigma.vals <- c(0.01, 0.5, 20)
n <- 1000
choices <- expand.grid(x_0.vals, n, sigma.vals)
choices_name <- paste0("x_0=", choices[, 1], " sigma=", choices[, 3])
result <- do.call(mapply, c(mc.algorithm, unname(choices)))
result <- data.frame(result)
colnames(result) <- choices_name

library(reshape)
hist_result <- melt(result)
colnames(hist_result) <- c("Parameter_Choices", "Value")


hist_result$line_x <- rep(seq(-2, 7.99, 0.01), times = 9)
hist_result$line_y <- rep(f(seq(-2, 7.99, 0.01), mn, std), times = 9)

Sample the target distribution using the random walk Metropolis-Hastings sampling with different choices of x_0 and sigma(the red line is the true target density):

ggplot(data = hist_result) + geom_histogram(aes(x = Value, color = Parameter_Choices, 
    y = ..density..), fill = "white") + geom_line(aes(x = line_x, y = line_y, color = "red")) + 
    facet_wrap(~Parameter_Choices, scales = "free") + theme(legend.position = "none") + 
    labs(title = "Explore The Impact Of Choices Of `x_0` And `sigma`", x = "x", y = "")

According to the result, we will see that x_0=-1, sigma=0.5 can best fit the true target density.