bayesian_MCMC.R

#!/usr/bin/Rscript
rm(list=ls())
library(plyr)


# Observed data. Taken from a normal distribution with
# mean 0, and s.d. 2.5
observations = read.csv('distribution.txt',row.names=1)


# Number of MCMC generations to explore the
# posterior distribution
generations = 1:1000

# Calculates the sum of log likelihoods
# for the probability that each observation
# is taken from the normal distribution
estimate_likelihood <- function(a,b){
  sum(apply(observations,1,function(x){
    x = pnorm(x,mean=a,sd=b)
    if(is.nan(x)){ x = 10^-100}
    if(x <= 0)   { x = 10^-100}
    log(x)
  }))
}

# Calculate the probability that the observed
# value is taken from the uniform distribution
# in the given limits
estimate_probability <- function(a,b){
  # Prior probabilities. Both uniform but with different limits
  (punif(a,min=-5,max=5) + punif(b,min=0,max=5)) / 2
}

# Randomly mutate a passed variable half the time
mutate <- function(x){
  if(runif(1) < 0.5) {
    x * runif(1,min=0.5,max=1.5)
  } else {
    x
  }
}

# Starting values
alpha     = runif(1)
beta      = runif(1)

result = adply(generations,1,function(x){

  # likelihood for initial values
  likelihood      = estimate_likelihood(alpha,beta)
  probability     = estimate_probability(alpha,beta)

  # Create two new variables to test and calculate the likelihood
  # probability for each
  new_alpha       = mutate(alpha)
  new_beta        = mutate(beta)
  new_likelihood  = estimate_likelihood(new_alpha,new_beta)
  new_probability = estimate_probability(new_alpha,new_beta)

  # Determine the proposal ratio between the original and
  # test values then determine if to accept the values
  proposal.ratio  = (likelihood * probability)/
    (new_likelihood * new_probability)
  accept          = proposal.ratio >= 0.975

  # Collect the results
  out = data.frame(
    alpha           = alpha,
    beta            = beta,
    new_alpha       = new_alpha,
    new_beta        = new_beta,
    new.likelihood  = new_likelihood,
    new.probability = new_probability,
    likelihood      = likelihood,
    probability     = probability,
    proposal.ratio  = proposal.ratio,
    accept          = accept
  )

  # Update the values if the proposal ratio is accepted
  if(accept == TRUE){
    alpha       <<- new_alpha
    beta        <<- new_beta
  }
  out
})

create_distribution.R

write.csv(file='distribution.txt',data.frame(observation=rnorm(200,sd=2.5)))