#!/usr/bin/env Rscript

# Probability mass function

p <- function(x) {
  if ( x == 0 ) return(1/6)
  if ( x == 1 ) return(2/6)
  if ( x == 2 ) return(3/6)
  return(0)
}

# Implementation of the Metropolis algorithm

number.of.accepts <- 0
state <- 0
draws <- rep(NA,10000)
for ( i in 1:length(draws) ) {
  cat("Iteration:                ",i,"\n")
  cat("Current state:            ",state,"\n")
  proposal <- state + ( 2*rbinom(1,n=1,p=0.99) - 1 )
  cat("Proposed state:           ",proposal,"\n")
  metropolis.ratio <- p(proposal)/p(state)
  cat("Metropolis ratio:         ",metropolis.ratio,"\n")
  u <- runif(1)
  cat("Random uniform:           ",u,"\n")
  if ( u < metropolis.ratio ) {
    cat("Accepted\n")
    number.of.accepts <- number.of.accepts + 1
    state <- proposal
  } else {
    cat("Rejected\n")
  }
  cat("Acceptance rate:          ",number.of.accepts/i,"\n")
  draws[i] <- state
  print(table(draws)/i)
  #system("sleep 1")
  cat("\n\n")
}


# Theoretical exercise

# Update the probabilities
run <- function(A,x,nreps) {
  for ( i in 1:nreps ) {
    cat(i,": ",paste(round(x,5),collapse=" "),"\n")
    x <- x %*% A
  }
}

# Define the transition kernel...
A <- matrix(NA,nrow=3,ncol=3)
A[1,] <- c(1/2,  1/2,  0)
A[2,] <- c(1/4,  1/4,  1/2)
A[3,] <- c(0,    1/3,  2/3)
A

# Notice that it doesn't matter where you start from!
run(A,c(1,0,0),200)
run(A,c(0,1,0),200)
run(A,c(0,0,1),200)

run(A,c(1/3,1/3,1/3),20)