simple_MCMC.R

library(plyr)

num_data <- 100
true_lambda <- 1

sample_from_distribution <- function(candidates, distribution, threshold) {
  while (TRUE) {
    x <- sample(candidates, size = 1)
    dice <- runif(n = 1, min = 0, max = 1)
    if (dice < distribution(x) / threshold) {
      return(x)
    }
  }
}

likelihood <- function(y, lambda) {
  dpois(y, lambda = lambda, log = TRUE)
}

y <- rpois(n = num_data, lambda = true_lambda)
lambda_candidates <- seq(from = 0.0, to = 3.0, by = 0.01)
dist <- ldply(lambda_candidates, function(lambda) {
  cbind(
    data.frame(lambda = lambda, likelihood = likelihood(y, lambda)),
    data.frame(y = y)
  )
})
lambda_to_likelihood <- aggregate(dist$likelihood, by=list(lambda=dist$lambda), FUN = function(x) exp(sum(x)))
lambda_density <- function(lambda) approxfun(lambda_to_likelihood$lambda, lambda_to_likelihood$x)(lambda)
threshold <- max(lambda_to_likelihood$x) * 1.1
num_samples <- 1024
hist(sapply(seq(num_samples), function(x) {
  sample_from_distribution(lambda_candidates, lambda_density, threshold)
}), probability = TRUE)

poisson_dist <- function(x) {
  dpois(x, 1)
}

# validate sample_from_distribution
numerical_distribution <- sapply(seq(1024), function(x) {
  sample_from_distribution(0:4, poisson_dist)
})
hist(numerical_distribution, probability = TRUE)
hist(rpois(1024, lambda = 1), probability = TRUE)