# Metropolis-Hastings toy problem
# Target density is Exponential(1)
rm(list=ls()); set.seed(20210310);
S <- 1000
MHIS.update <- function(theta0, lambda)
{
theta.star <- rexp(1, lambda)
r <- exp( (lambda - 1) * (theta.star - theta0) )
if(runif(1) < r){ theta0 <- theta.star }
theta0
}
MRW.update <- function(theta0, sigma)
{
theta.star <- rnorm(1, mean=theta0, sd=sigma)
if(theta.star > 0)
{
if(runif(1) < exp( -(theta.star - theta0)))
{
theta0 <- theta.star
}}
theta0
}
# Now run a few!
op <- par(mfrow=c(2,3))
theta.chain <- rep(NA, S)
theta0 <- 1
for(s in 1:S)
{
theta0 <- MHIS.update(theta0, lambda=1/2)
theta.chain[s] <- theta0
}
plot.ts(theta.chain, ylim=c(0,9), main="MHIS with lambda=.5")
mean(abs(theta.chain[-1] - theta.chain[-S]) > 0)
chain1 <- theta.chain theta.chain <- rep(NA, S) theta0 <- 1 for(s in 1:S) { theta0 <- MHIS.update(theta0, lambda=1) theta.chain[s] <- theta0 } plot.ts(theta.chain, ylim=c(0,9), main="MHIS with lambda=1") mean(abs(theta.chain[-1] - theta.chain[-S]) > 0)
chain2 <- theta.chain theta.chain <- rep(NA, S) theta0 <- 1 for(s in 1:S) { theta0 <- MHIS.update(theta0, lambda=2) theta.chain[s] <- theta0 } plot.ts(theta.chain, ylim=c(0,9), main="MHIS with lambda=2") mean(abs(theta.chain[-1] - theta.chain[-S]) > 0)
chain3 <- theta.chain theta.chain <- rep(NA, S) theta0 <- 1 for(s in 1:S) { theta0 <- MRW.update(theta0, sigma=.2) theta.chain[s] <- theta0 } plot.ts(theta.chain, ylim=c(0,9), main="MRW with sigma=.2") mean(abs(theta.chain[-1] - theta.chain[-S]) > 0)
chain4 <- theta.chain theta.chain <- rep(NA, S) theta0 <- 1 for(s in 1:S) { theta0 <- MRW.update(theta0, sigma=1) theta.chain[s] <- theta0 } plot.ts(theta.chain, ylim=c(0,9), main="MRW with sigma=1") mean(abs(theta.chain[-1] - theta.chain[-S]) > 0)
chain5 <- theta.chain theta.chain <- rep(NA, S) theta0 <- 1 for(s in 1:S) { theta0 <- MRW.update(theta0, sigma=5) theta.chain[s] <- theta0 } plot.ts(theta.chain, ylim=c(0,9), main="MRW with sigma=5") mean(abs(theta.chain[-1] - theta.chain[-S]) > 0)
chain6 <- theta.chain rm(theta.chain) par(op)