####################################################
#    Scripts for Chapter 14 of
#    P. Brandimarte
#    Handbook in Monte Carlo Simulation
#    Applications in Financial Engineering,
#    Risk Management, and Economics
#    Wiley, 2014

#  Required packages
require(car)

#  Source the function file
source("Chapter14_Functs.R")

####################################################
#   Chunk 1: An illustration of the Metropolis--Hastings algorithm
####################################################
# define mixture of two normals
muA <- 1; sigmaA <- 3; muB <- 10; sigmaB <- 2; p <- 0.7
# define the unscaled target density
target <- function(theta) 20*(p*exp(-(theta-muA)^2/
                                      sigmaA^2/2)/sigmaA 
                              + (1-p)*exp(-(theta-muB)^2/sigmaB^2/2)/sigmaB)
# set candidate density
muQ <- 0; sigmaQ <- 3
# set starting value
theta0 <- (muA+muB)/2
burninSteps <- 1000; numSteps <- 50000
set.seed(55555)
# generate random variates, including BurnIn
Zb <- rnorm(burninSteps, muQ, sigmaQ)
Ub <- runif(burninSteps)
Z <- rnorm(numSteps, muQ, sigmaQ)
U <- runif(numSteps)
# first Burn In
oldTheta <- theta0
for (i in 1:burninSteps){
  tryTheta <- oldTheta + (muQ+Zb[i]*sigmaQ)
  alpha <- min(1,target(tryTheta)/target(oldTheta))
  if (Ub[i] <= alpha) # accept
    oldTheta <- tryTheta
}
# now collect info about target sample
targetSample <- numeric(numSteps)
for (i in 1:numSteps){
  tryTheta <- oldTheta + (muQ+Z[i]*sigmaQ)
  alpha <- min(1,target(tryTheta)/target(oldTheta))
  if (U[i] <= alpha)# accept
    oldTheta <- tryTheta
  targetSample[i] <- oldTheta
}
# plot actual density and histogram
theta<-seq(from=min(muA-3*sigmaA, muB-3*sigmaB),
           to=max(muA+3*sigmaA, muB+3*sigmaB),by=0.001)
trueF <- p*dnorm(theta, muA, sigmaA)+ 
  (1-p)*dnorm(theta,muB,sigmaB)
h<-hist(targetSample, plot = FALSE, nclass = 50)
ymax<-max(c(h$density,trueF))*1.05
hist(targetSample, prob = TRUE, col = "light blue",
     xlim = range(theta), ylim = c(0,ymax),
     main = "Sample from target density",
     nclass = 50, xlab = 'x', ylab = 'Density')
lines(theta, trueF)


####################################################
#   Chunk 2: An illustration of a Gibbs sampler
####################################################

set.seed(55555)
out <- gibbsNormal(10000,500,0.85)
scatterplotMatrix(out,diagonal="histogram",
                  smooth=FALSE,reg.line=FALSE,pch=20)
