# 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
#######################################

# This yields a chain which is far too "cold", meaning that most candidates are accepted, but they represent 
# conservative moves within the sample space. 
run2 <- MH_norm(1000,0.1)

# Plot sampled values as a time series and histogram (partition the plot window into 2 panels)
par(mfrow=c(2,1))
plot(ts(run2), xlab="Iteration", ylab="Sampled Values", col="darkgreen",
     main="Metropolis-Hastings MCMC: Normal Target Distribution")
hist(run2, breaks=30, freq=FALSE, col="wheat", main="Histogram of Sampled Values", xlab="Sampled Values")
par(mfrow=c(1,1))


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

# This yields a chain which is far too "hot". Very few candidate points are accepted, but when they are 
# they often represent quite considerable moves within the sample space. 
run3 <- MH_norm(5000,200)

# Plot sampled values as a time series and histogram (partition the plot window into 2 panels)
par(mfrow=c(2,1))
plot(ts(run3), xlab="Iteration", ylab="Sampled Values", col="darkgreen",
     main="Metropolis-Hastings MCMC: Normal Target Distribution")
hist(run3, breaks=30, freq=FALSE, col="wheat", main="Histogram of Sampled Values", xlab="Sampled Values")
par(mfrow=c(1,1))


#######################################################################################
### 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(0,1)

# Metropolis-Hastings sampler for a gamma target distribution based on normal proposal distribution

MH_gamm <- function(n, a, b) 
{
  mu <- a/b
  sigma <- sqrt(a/(b*b))
  samples <- vector("numeric", n)
  x <- a/b
  samples[1] <- x  # initial value is the mean of the gamma distribution
  for (i in 2:n) {
    proposal <- x + rnorm(1,0,1)  
    aprob <- min(1, (dgamma(proposal,a,b)/dgamma(x,a,b)))
    u <- runif(1)
    if (u < aprob){
      x <- proposal
    } 
    samples[i] <- x
  }
  samples
}


# This shows a well mixing chain and a reasonable gamma-like distribution of sampled values
run4 <- MH_gamm(1000,2.3,2.7)

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


sample_mean = mean(run4); sample_mean
2.3/2.7



#######################################################################################
### Exercise 4 - Independence Chain
### 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 INDEPENDENCE sampler for a gamma target distribution based on normal candidates with 
# the same mean and variance as the gamma distribution. Independence chain version means that the proposed 
# transition is formed independently of the previous position of the chain, so Q_ij = f(x) for some density
# function f - in this case, Normal(a/b, a/(b*b)).

MH_gamm_indep <- function(n, a, b) 
{
  mu <- a/b
  sigma <- sqrt(a/(b*b))
  samples <- vector("numeric", n)
  x <- a/b
  samples[1] <- x  # initial value is the mean of the gamma distribution
  for (i in 2:n) {
    proposal <- rnorm(1, mu, sigma)
    aprob <- min(1, (dgamma(proposal,a,b)/dgamma(x,a,b))/(dnorm(proposal,mu,sigma)/dnorm(x,mu,sigma)))
    u <- runif(1)
    if (u < aprob){
      x <- proposal
    } 
    samples[i] <- x
  }
  samples
}


# This shows a well mixing chain and a reasonable gamma-like distribution of sampled values
run5 <- MH_gamm_indep(1000,2.3,2.7)

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


sample_mean = mean(run5); sample_mean
2.3/2.7