# Random Walk Metropolis algorithm for MacKay's (2003, Fig. 29.12) 'toy' problem
# MacKay wants "to illustrate how slowly a random walk explores a state space"(p.369)
# R-Code by PCM 2016/02/10
library(lattice)
## Warning: package 'lattice' was built under R version 3.2.3
linear_density_graph <- function(r=-0.5, s=20.5) {
# uniform density, MacKay, p.369
curve((x-x+1)*(1/21), from=r, to=s, type="l", xlab="x", ylab="f(x)", lwd=3,
main="Uniform Proposal Density p(x) = 1/21", col="steelblue")
# vertical lines from (xi1, yi1=0) to (xi2, yi2=1/21)
for(i in r:s) {lines(c(i,i),c(0,1/21),type="l",col="steelblue",lwd=1)}
legend(x=2,y=0.15, legend="area under density f(x) is 1 ")
}
linear_density_graph()
P <- function(x) # uniform target pmf P(x) (29.33)
{if (x %in% c(0:20)) 1.0/21.0 else 0}
flip <- function(p=0.5) rbinom(1,1,p) # one Bernoulli trial (inspired by WebChurch)
rQ <- function(x) # proposal function
{x + flip()*2-1} # rQ ~ Q(xp | x) (29.34)
pQ <- function(xp,x) # proposal density
{if(abs(xp-x)==1) 0.5 else 0} # pQ(xp,x) = Q(rQ | x) (29.34)
# lower bound on number of iterations of a Metropolis method (29.32)
T <- function(len, eps) (len/eps)^2 # MacKay's "Rule-of-Thumb", 2003, p.367
cat("\nRandom walk Metropolis-Hastings (MH) algorithm, MacKay, 2003, Fig.29.12\n")
##
## Random walk Metropolis-Hastings (MH) algorithm, MacKay, 2003, Fig.29.12
cat(" R-Code by PCM 2016/02/10\n")
## R-Code by PCM 2016/02/10
nt <- 480000 # total nr of steps
nthin <- 300 # thin-factor
maxsteps <- 100000 # max MH-steps for trace-plot
cat("total number of MH-steps = ", nt)
## total number of MH-steps = 480000
T_est <- T(20,1)
print(sprintf("MacKay's rule of thumb %f iterations for one independent sample", T_est))
## [1] "MacKay's rule of thumb 400.000000 iterations for one independent sample"
print(sprintf("MacKay's rule of thumb %f independent samples within total of %i samples", nt/T_est, nt))
## [1] "MacKay's rule of thumb 1200.000000 independent samples within total of 480000 samples"
steps <- c(1:nt) # abcissa in plot
flips <- matrix(0,nt)
for(i in 1:nt) {flips[i] <- flip()} # Bernouilli trials
hist(flips,breaks=seq(-0.5,1.5,1),main="Distribution of Steps (left = -1 and right = +1)",col="steelblue")
x <- matrix(10,nt) # states x[t]; start state x[1] = 10
xpt <- matrix(0,nt) # proposals xp[t] = x'[t]
importance_ratio <- matrix(0,nt) # initials
proposal_ratio <- matrix(0,nt) # initials
alpha <- matrix(0,nt) # initial alpha
accept <- matrix(0,nt) # initially no acceptance
p_acc <- matrix(0,nt) # initial P(acceptance)
#
# MH-loop
#
for (t in 1:nt) {
xt <- x[t]
xp <- rQ(xt) # xp = x' = proposal
xpt[t] <- xp
importance_ratio[t] <- P(xp)/P(xt)
proposal_ratio[t] <- pQ(xp,xt)/pQ(xt,xp)
alpha[t] <- importance_ratio[t]/proposal_ratio[t]
p_acc[t] <- min(1, alpha[t]) # P(Acceptance)
x[t+1] <-
if(flip(p_acc[t]) == TRUE)
{accept[t] <- 1; xp; } # accept ...
else x[t] # ... or reject proposal xp = x'
}
if(nt <= maxsteps) plot(steps,xpt,type="l",col="steelblue",main="Random walk MH: trace of proposals xp (MacKay, p.368)")
hist(xpt,main=sprintf("Distribution of %i proposed states xp",nt),breaks=c(-2.5:22.5),col="steelblue")
if(nt <= maxsteps) plot(steps,x[-(nt+1)],type="l",main="Random walk MH: trace of states x (MacKay, p.368)",col="steelblue")
hist(x,main=sprintf("Distribution of States (%i samples)",nt),breaks=c(-0.5:20.5),col="steelblue")
acf(x,lag.max=min(nt*.10,1000),main="ACF(x)", col="steelblue")
if(nt <= maxsteps) plot(steps,importance_ratio,type="l",main="Trace: importance_ratio[t] <- P(xp) / P(xt)",col="steelblue")
if(nt <= maxsteps) densityplot(importance_ratio,main="Density plot: importance_ratio[t] <- P(xp) / P(xt)",col="steelblue",lwd=2)
if(nt <= maxsteps) plot(steps,proposal_ratio,type="l",main="Trace: proposal_ratio[t] <- pQ(xp,xt) / pQ(xt,xp)",col="steelblue")
if(nt <= maxsteps) densityplot(proposal_ratio,main="Density plot: proposal_ratio[t] <- pQ(xp,xt) / pQ(xt,xp)",col="steelblue",lwd=2)
if(nt <= maxsteps) plot(steps,alpha,type="l",main="Trace: alpha <- importance_ratio / proposal_ratio",col="steelblue")
if(nt <= maxsteps) densityplot(alpha,main="Density plot: alpha <- importance_ratio / proposal_ratio",col="steelblue",lwd=2)
if(nt <= maxsteps) plot(steps,p_acc,type="l",main="Trace: P(Acceptance) <- min(1,alpha)",col="steelblue")
percentage_of_acceptance <- mean(accept)*100.0
acceptance_rate <- mean(p_acc[])
if(nt <= maxsteps) legend(x=0,y=0.4,legend=sprintf("PercentageOfTotalAcceptance = %f",percentage_of_acceptance))
if(nt <= maxsteps) legend(x=0,y=0.2,legend=sprintf("Acceptance Rate = %f",acceptance_rate))
if(nt <= maxsteps) densityplot(p_acc,main="Density plot: P(Acceptance) <- min(1,alpha)",col="steelblue",lwd=2)
print(sprintf("Percentage of acceptance = %f", percentage_of_acceptance))
## [1] "Percentage of acceptance = 95.310000"
print(sprintf("acceptance-rate = %f", acceptance_rate))
## [1] "acceptance-rate = 0.953100"
# excluding samples by thinning
nthinned <- round(nt/nthin)
print(sprintf("Number of samples (after thinning) = %i", nthinned))
## [1] "Number of samples (after thinning) = 1600"
x <- x[seq(from=1,to=nt,by=nthin)]
accept <- accept[seq(from=1,to=nt,by=nthin)]
p_acc <- p_acc[seq(from=1,to=nt,by=nthin)]
print(sprintf("Percentage(AcceptanceDecisions) = %f (Thinned)", mean(accept)*100.0))
## [1] "Percentage(AcceptanceDecisions) = 95.562500 (Thinned)"
print(sprintf("AcceptanceRate = %f (Thinned)", mean(p_acc)))
## [1] "AcceptanceRate = 0.955625 (Thinned)"
if(nt <= maxsteps) plot(x, type="l", main="MH-Trace for x (Thinned)", col="steelblue")
hist(x,main=sprintf("Distribution of (thinned) States (%i samples)",nthinned),breaks=c(-0.5:20.5),col="steelblue")
acf(x, main="ACF(m) (Thinned)")
x_est <- mean(x)
sdx_est <- sd(x)
print(sprintf("mean(x) = %f (Thinned)", x_est))
## [1] "mean(x) = 9.790625 (Thinned)"
print(sprintf("sd(x) = %f (Thinned)", sdx_est))
## [1] "sd(x) = 6.032685 (Thinned)"
#