cauchy.one.par.MCMC.R

albert — Oct 7, 2013, 10:38 AM

# define log posterior of Cauchy(mu, 1) sampling, uniform prior

post.cauchy <- function(theta, y){
  sum(dcauchy(y, location=theta, scale=1, log=TRUE))
}

# some data

y1 <- c(2.5, 4, 3, 4, 6, 4, 7.6) + 1
y2 <- c(12, 13, 12.5, 14, 13, 11)
y <- c(y1, y2)

# plot posterior

mu <- seq(0, 15, length=200)
plot(mu, exp(sapply(mu, post.cauchy, y)), type="l")
rug(y)

# laplace fitting -- we find one of the modes

library(LearnBayes)
laplace(post.cauchy, 4, y)

$mode
[1] 5.2

$var
       [,1]
[1,] 0.2311

$int
[1] -45.04

$converge
[1] TRUE

# get that posterior is approx N(5.2, 0.23)

# illustrate two random walk chains
# one using (var=0.22, scale=1) and one using (var=0.22, scale=10)

rw.fit1 <- rwmetrop(post.cauchy, 
                   proposal=list(var=0.22, scale=1),
                   start=4, 
                   m=1000, 
                   y)

rw.fit2 <- rwmetrop(post.cauchy, 
                    proposal=list(var=0.22, scale=10),
                    start=4, 
                    m=1000, 
                    y)

# use the coda package to compare the streams of simulated draws


library(coda)

Loading required package: lattice

plot(mcmc(rw.fit1$par))

plot(mcmc(rw.fit2$par))

# rerun the random walk using scale=10, collecting 100,000 iterates

rw.fit <- rwmetrop(post.cauchy, 
                   proposal=list(var=0.22, scale=10),
                   start=4, 
                   m=100000, 
                   y)

str(rw.fit)  # acceptance rate is 37%

List of 2
 $ par   : num [1:100000, 1] 4 4.74 6.79 6.79 6.65 ...
 $ accept: num 0.33

# compare density plot of simulated draws with exact posterior

mu <- seq(0, 15, by=0.1)
f.mu <- sapply(mu, post.cauchy, y)
f.mu <- exp(f.mu) / sum(exp(f.mu)) / 0.1

plot(mu, f.mu, type="l", lwd=3)

lines(density(rw.fit$par), lwd=2, col="red")