mix.exponentials.R

Jim — Oct 2, 2013, 8:00 PM

# Mixture of exponentials problem
# Chapter 5, exercises 6

log.exponential.mix <- function(theta, y){
  lambda.A <- exp(theta[1])
  lambda.B <- exp(theta[2])
  sum(log(.8 *dexp(y, 1/lambda.A) + 
            (1 - .8) * dexp(y, 1/lambda.B)))
}

y <- c(9.3,4.9,3.5,26.0,0.6,1.0,3.5,26.9,
    2.6,20.4,1.0,10.0,1.7,11.3,7.7,14.1,
    24.8,3.8,8.4,1.1,24.5,90.7,16.4,30.7,
    8.5,5.9,14.7,0.5,99.5,35.2)

# construct contour plot of posterior

library(LearnBayes)
mycontour(log.exponential.mix,c(1,4,-2, 8),y)

# since the posterior is bimodal, we'll get different estimates
# of the posterior mode depending on the starting value

mode.1 <- laplace(log.exponential.mix, c(3, 0), y)$mode
points(mode.1[1], mode.1[2], pch="x")

mode.2 <- laplace(log.exponential.mix, c(2, 4), y)$mode
points(mode.2[1], mode.2[2], pch="x")

# illustrate random walk Metropolis

sim1 <-  rwmetrop(log.exponential.mix,
         proposal=list(var=diag(c(4, 10)), scale=0.5),
         start=c(2, 4),
         m=10000, y)
str(sim1)

List of 2
 $ par   : num [1:10000, 1:2] 2 2 2 2 2 ...
 $ accept: num 0.183

library(coda)

Loading required package: lattice

plot(mcmc(sim1$par))

mycontour(log.exponential.mix,c(1,4,-2, 8),y)
points(sim1$par)

sim2 <-  rwmetrop(log.exponential.mix,
                  proposal=list(var=diag(c(4, 10)), scale=0.1),
                  start=c(2, 4),
                  m=10000, y)

plot(mcmc(sim2$par))