#############################################################################################################
#
# Exercise 6.9 (page 236)  - Univariate version of Example 6.2 (pages 199 and 202-3) 
#
# Metropolis paths when the target is a mixture of two univariate normal densities
#
##############################################################################################################
#
# Author : Hedibert Freitas Lopes
#          Graduate School of Business
#          University of Chicago
#          5807 South Woodlawn Avenue
#          Chicago, Illinois, 60637
#          Email : hlopes@ChicagoGSB.edu
#
###############################################################################################################
# Univariate Mixture of two Normal densities
# ------------------------------------------
p    = function(x){w*dnorm(x,mu[1],sig[1])+(1-w)*dnorm(x,mu[2],sig[2])}
mu   = c(0,3.5)
sig  = c(1,0.5)
w    = 0.9
x2   = seq(-5,7,length=1000)
f    = w*dnorm(x2,mu[1],sig[1])+(1-w)*dnorm(x2,mu[2],sig[2])


# Random Walk Metropolis
# ----------------------
set.seed(124763)
inis = c(-5,7)
sds  = c(0.1,0.5)
M    = 10000
xx   = NULL
for (sd in sds){
  for (ini in inis){
    x  = ini
    xs = x
    for (iter in 1:M){
      x1 = rnorm(1,x,sd)
      accept = min(1,p(x1)/p(x))
      if (runif(1)<=accept){
        x = x1
      }
      xs = c(xs,x)
    }
    xx = cbind(xx,xs)
  }
}
# Figure 6.5
# ----------
x1 = seq(-5.3,7.2,length=30)
par(mfrow=c(4,2))
i = 0
for (sd in sds){
  for (ini in inis){
    i = i + 1
    plot(xx[,i],xlab="iterations",ylab="",main="",ylim=range(x1),axes=F,type="l")
    axis(1);axis(2)
    title(paste("SD=",sd,"  -  initial value=",xx[1,i],sep=""))
    hist(xx[,i],xlab="",ylab="",prob=T,breaks=x1,main="",ylim=c(0,0.4))
    lines(x2,f,col=2)
    title(paste("Acceptance rate=",round(100*length(unique(xx[,i]))/M),"%",sep=""))
  }
}


# Independent Metropolis-Hastings
# -------------------------------
set.seed(124763)
m    = 0
inis = c(-5,7)
sds  = c(1.0,3.0)
M    = 10000
xx   = NULL
for (sd in sds){
  for (ini in inis){
    x  = ini
    xs = x
    for (iter in 1:M){
      x1 = rnorm(1,m,sd)
      accept = min(1,p(x1)/p(x)*dnorm(x,m,sd)/dnorm(x1,m,sd))
      if (runif(1)<=accept){
        x = x1
      }
      xs = c(xs,x)
    }
    xx = cbind(xx,xs)
  }
}


# Figure 6.6
# ----------
x1 = seq(-5.3,7.2,length=30)
par(mfrow=c(4,2))
i = 0
for (sd in sds){
  for (ini in inis){
    i = i + 1
    plot(xx[,i],xlab="iterations",ylab="",main="",ylim=range(x1),axes=F,type="l")
    axis(1);axis(2)
    title(paste("SD=",sd,"  -  initial value=",xx[1,i],sep=""))
    hist(xx[,i],xlab="",ylab="",prob=T,breaks=x1,main="",ylim=c(0,0.4))
    lines(x2,f,col=2)
    title(paste("Acceptance rate=",round(100*length(unique(xx[,i]))/M),"%",sep=""))
  }
}