# Markov Chain Monte Carlo (MCMC) Examples -- Ref: Ravenzwaaif et al 2016 (MCMC-intro.pdf)
# STAT 753
# March 3, 2020


#########################################
# Metropolis-Hastings algorithm for MCMC
#########################################

# Draw samples from a target distribution via MH algorithm. 
# Define target and proposal distributions.


####################################################################
### Example 1 uses target distribution = Normal(mean=100,sd=15) and
### proposal distribution = Normal(mean=0,sd=5).

n = 500;              # Define number of samples
samples = numeric(n)  
samples[1] = 110      # Initial guess (reasonable)
#samples[1] = 250      # Initial guess (starting in tail of distribution)
#samples[1] = 600      # Initial guess (starting with a value far from true distribution)

for (i in 2:n){
  proposal = samples[i-1] + rnorm(1,0,5)  # Proposal value
  if ((dnorm(proposal,100,15) / dnorm(samples[i-1],100,15)) > runif(1))
    samples[i] = proposal                 # Accept proposal value
  else (samples[i] = samples[i-1])        # Reject proposal value
}

# Plot sampled values
plot(1:n, samples, type="b", pch=19, col="blue", xlab="Iteration", ylab="Sampled Values",
     main="Metropolis-Hastings MCMC: Normal(100,15) Target Distribution")

# Plot density (histogram) of sampled values
hist(samples, freq=FALSE, col="wheat", xlab="Sampled Values", main="Histogram of Sampled Values")

# Calculate sample mean
sample_mean = mean(samples); sample_mean



####################################################################
### Example 2 uses target distribution = Normal(mean=0,sd=1) and
### proposal distribution = Uniform(-alpha,alpha).

MH_norm <- function(n, alpha) 
{
  samples <- vector("numeric", n)
  x <- 0
  #x <- 10  # What happens if your initial value is far from the true distribution?
  samples[1] <- x
  for (i in 2:n) {
    increment <- runif(1, -alpha, alpha)
    proposal <- x + increment
    aprob <- min(1, dnorm(proposal,0,1)/dnorm(x,0,1))
    u <- runif(1)
    if (u < aprob){
      x <- proposal
    } 
    samples[i] <- x
  }
  samples
}


# Run the MH function above
run1 <- MH_norm(1000,1)

# Plot samples as a time series (partition the plot window into 2 panels)
par(mfrow=c(2,1))
plot(ts(run1), xlab="Iteration", ylab="Sampled Values", col="darkgreen",
     main="Metropolis-Hastings MCMC: Normal Target Distribution")

# Plot a histogram of sampled values
hist(run1, breaks=30, freq=FALSE, col="wheat", main="Histogram of Sampled Values", xlab="Sampled Values")
par(mfrow=c(1,1))



#######################################################################################
### Exercise 1 
### For Example 2 above, see how different choices of alpha affect the "mixing" of the chain.
### Try alpha = 0.1, 200, anything in between.

# Run the MH function above: alpha=0.1
#######################################


# Run the MH function above: alpha=200
#######################################



#######################################################################################
### Exercise 2 
### Modify the R function (Example 2) so that it records and then prints the overall acceptance rate of the chain, 
### as well as returning the vector of simulated values. "Tune" the MH sampler by finding a value of alpha which 
### results in an acceptance probability of about 0.3.


#######################################################################################
### Exercise 3 
### Modify Example 2 to use target distribution = Gamma(alpha,beta) and proposal distribution = Normal with the
### same mean and SD as the desired gamma distribution, i.e. Normal(a/b, a/(b*b)).


# Metropolis-Hastings sampler for a gamma target distribution based on normal candidates with 
# the same mean and variance as the gamma distribution

MH_gamm <- function(n, a, b) 
{
  # fill this in using Example 2 as a template!
}


# This shows a well mixing chain and a reasonable gamma-like distribution of sampled values
run <- MH_gamm() # fill in arguments

# Plot
par(mfrow=c(2,1))
plot(ts(run), xlab="Iteration", ylab="Sampled Values", col="tomato", 
     main="Metropolis-Hastings MCMC: Gamma Target Distribution")
hist(run, breaks=30, freq=FALSE, col="wheat", main="Histogram of Sampled Values", xlab="Sampled Values")
par(mfrow=c(1,1))


sample_mean = mean(run); sample_mean