# define the log-likelihood function 
log.lik <- function(theta, y){

  if((theta<1)&(theta>0)) {
    loglik <- y[1]*log(2+theta)+(y[2]+y[3])*log(1-theta)+y[4]*log(theta)
  }
  else {
    # if theta is not in the defined range return NA
    loglik <- NA
  }
  return(loglik)
}

mcmc_indep <- function(M, y, F){

  # store samples here
  xsamples <- rep(NA, M)

  # Idea: Normal Independence proposal with mean equal to the posterior
  # mode and standard deviation equal to the standard error or to 
  # a multiple of the standard error.
  mymean <- optimize(log.lik, y=y, lower = 0, upper = 1, maximum = TRUE)$maximum

  # negative curvature of the log-posterior at the mode
  a <- -1*(-y[1]/(2+mymean)^2 - (y[2]+y[3])/(1-mymean)^2 - y[4]/mymean^2)
  mystd <- sqrt(1/a)
  
  ##################################################################
  ## Alternative using optim instead of optimize. Optim directly returns the Hessian.
#     eps <- 1E-9
#     mycontrol <- list(fnscale=-1, maxit=100)
#   
#     ml <- optim(0.5, log.lik, y=data, control=mycontrol, hessian=T, 
#                 method="L-BFGS-B", lower=0+eps, upper=1-eps)
#     mymean <- ml$par
#     mystd <- sqrt(-1/ml$hessian)
  ##################################################################

  # factor blows up standarddeviation
  factor=F
  yes <- 0
  no <- 0

  # use as initial starting value mymean
  xsamples[1] <- mymean

  # MH-Iteration
  for(k in 2:M){
    # value of the past iteration
    old <- xsamples[k-1]

    # propose new value
    proposal <- rnorm(1, mean=mymean, sd=mystd*factor)

    # compute acceptance ratio
    posterior.ratio <-  exp(log.lik(proposal, data)-log.lik(old, data))
    if(is.na(posterior.ratio)){
      # happens when the proposal is not between 0 and 1 => acceptance probability will be 0
      posterior.ratio <- 0
    }
    proposal.ratio <- exp(dnorm(old, mymean, mystd*factor, log=T)-
      dnorm(proposal, mymean, mystd*factor, log=T))

    # get the acceptance probability
    alpha <- posterior.ratio*proposal.ratio

    # accept-reject step
    if(runif(1) <= alpha){	
      # accept the proposed value
      xsamples[k] <- proposal
      # increase counter of accepted values
      yes <- yes + 1
    }
    else{
      # stay with the old value
      xsamples[k] <- old
      no <- no + 1
    }
  }

  # acceptance rate
  cat("The acceptance rate is: ", round(yes/(yes+no)*100,2), "%\n", sep="")
  return(xsamples)
}


mcmc_rw <- function(M, y){

  # store samples here
  xsamples <- rep(NA, M)

  # Idea: Use a random walk proposal: U(theta^(k) - d, theta^(k) + d)
  d <- sqrt(12)/2*0.1

  yes <- 0
  no <- 0

  # specify a starting value
  xsamples[1] <- 0.5

  # MH-Iteration
  for(k in 2:M){
    # value of the past iteration
    old <- xsamples[k-1]

    # propose new value
    proposal <- runif(1, old-d, old+d)

    # compute acceptance ratio
    posterior.ratio <-  exp(log.lik(proposal, data)-log.lik(old, data))
    if(is.na(posterior.ratio)){
      # happens when the proposal is not between 0 and 1 => acceptance probability will be 0
      posterior.ratio <- 0
    }
    # the proposal ratio is equal to 1 as we have a symmetric proposal distribution
    proposal.ratio <- 1

    # get the acceptance probability
    alpha <- posterior.ratio*proposal.ratio

    # accept-reject step
    if(runif(1) <= alpha){	
      # accept the proposed value
      xsamples[k] <- proposal
      # increase counter of accepted values
      yes <- yes + 1
    }
    else{
      # stay with the old value
      xsamples[k] <- old
      no <- no + 1
    }
  }

  # acceptance rate
  cat("The acceptance rate is: ", round(yes/(yes+no)*100,2), "%\n", sep="")
  return(xsamples)
}
## data
data <- c(125, 18, 20, 34)
# number of iterations
M <- 10000
xsamplesF1 <- mcmc_indep(M=M, F=1, y=data)
xsamplesF10 <- mcmc_indep(M=M, F=10, y=data)
xsamplesRW <- mcmc_rw(M=M, y=data)

## some plots
# Independence proposal F=1
par(mfrow=c(3,3))
# traceplot
plot(xsamplesF1, type="l")
# autocorrelation plot
acf(xsamplesF1)
# histogram
hist(xsamplesF1, nclass=100, xlim=c(0.4,0.8), prob=T)

# Independence proposal F=10
plot(xsamplesF10, type="l")
acf(xsamplesF10)
hist(xsamplesF10, nclass=100, xlim=c(0.4,0.8), prob=T)

# RW proposal
plot(xsamplesRW, type="l")
acf(xsamplesRW)
hist(xsamplesRW, nclass=100, xlim=c(0.4,0.8), prob=T)