CS计算机代考程序代写 chain # Metropolis-Hastings toy problem

# 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)