pdf(file="chap4-MH.pdf", paper="a4") mhmc = function(x1, T, f, rq, dq){ X = array(NA, T) X[1]=x1 naccept=0 for(t in 1:(T-1)){ y = rq(X[t]) x = X[t] if( ( runif(1) <= f(y)*dq(y,x)/(f(x)*dq(x,y)) ) ){ X[t+1] = y naccept = naccept+1 } else { X[t+1] = x } rm(x,y) } return(list(X=X, accept.prop=naccept/T)) } ################ MH for density N(0,1) with gaussian proposals main.title="MH for density N(0,1) with gaussian proposals" T=1e4 x1=10 f = function(x){ return(exp(-x^2/2)) } # par(mfrow=c(4, 3), mar=c(2, 2, 1, 1), oma=c(0,0,2,0), cex=0.6) # for(sd in c(.1, .5, 2.5, 10)){ rq = function(x){ x + sd*rnorm(1) } dq = function(x,y){ dnorm(y, mean=x, sd=sd) } # ans=mhmc(x1, T, f, rq, dq) # plot(ans$X, type='l', ylim=c(-2,10), main=paste("sd=", sd, ", acceptance ratio=", ans$accept.prop, sep="")) abline(h=c(-2,2), lty=2) lines(cumsum(ans$X)/(1:(T)), type='l', ylim=c(-2, 2), col="blue") # z=seq(-3,3, by=.1) plot(z, dnorm(z), col='red', type='l', xlim=c(-4,4)) hist(ans$X, breaks='FD', prob=TRUE, add=TRUE) # acf(ans$X) #abline(h=mean(X), col='blue') #abline(h=0, col='red') } ################ MH for density N(0,1) with asymmetric proposals main.title="MH for density N(0,1) with asymmetric proposals" T=1e4 x1=3 f = function(x){ return(exp(-x^2/2)) } # par(oma=c(0,0,2,0), mfrow=c(4, 3), mar=c(2, 2, 1, 1), cex=0.6) # for(shift in c(.1, .5, 2.5, 10)){ rq = function(x){ x + rchisq(1, df=3)- shift } dq = function(x,y){ dchisq(y - x + shift, df=3) } #rq = function(x){ x + rexp(1, rate=1)- shift } #dq = function(x,y){ dexp(y - x + shift, rate=1) } # ans=mhmc(x1, T, f, rq, dq) # plot(ans$X, type='l', ylim=c(-2,10), main=paste("shift=", shift, ", acceptance ratio=", ans$accept.prop, sep="")) abline(h=c(-2,2), lty=2) lines(cumsum(ans$X)/(1:(T)), type='l', ylim=c(-2, 2), col="blue") # z=seq(-3,3, by=.1) plot(z, dnorm(z), col='red', type='l', xlim=c(-4,4)) hist(ans$X, breaks='FD', prob=TRUE, add=TRUE) # z=seq(-10,10, by=.01) plot(z, dq(0,z), type='l') #acf(ans$X) #abline(h=mean(X), col='blue') #abline(h=0, col='red') } ################ MH for density N(0,1) with uniform proposals main.title="MH for density N(0,1) with uniform proposals" T=1e4 x1=10 f = function(x){ return(exp(-x^2/2)) } # par(oma=c(0,0,2,0), mfrow=c(4, 3), mar=c(2, 2, 1, 1), cex=0.6) # for(sd in c(.1, 1, 2.5, 20)){ rq = function(x){ runif(1, min=x-sd, max=x+sd) } dq = function(x,y){ dunif(y, min=x-sd, max=x+sd) } # ans=mhmc(x1, T, f, rq, dq) # plot(ans$X, type='l', ylim=c(-2,10), main=paste("sd=", sd, ", acceptance ratio=", ans$accept.prop, sep="")) abline(h=c(-2,2), lty=2) lines(cumsum(ans$X)/(1:(T)), type='l', ylim=c(-2, 2), col="blue") # z=seq(-3,3, by=.1) plot(z, dnorm(z), col='red', type='l', xlim=c(-4,4)) #lines(density(ans$X), col='blue', type='l', xlim=c(-4,4)) hist(ans$X, breaks='FD', prob=TRUE, add=TRUE) # acf(ans$X) #abline(h=mean(X), col='blue') #abline(h=0, col='red') } log.mhmc = function(x1, T, log.f, rq, log.dq){ X = array(NA, T) X[1]=x1 naccept=0 for(t in 1:(T-1)){ y = rq(X[t]) x = X[t] if( log(runif(1)) <= log.f(y) - log.f(x) + log.dq(y,x) - log.dq(x,y)){ X[t+1] = y naccept = naccept+1 } else { X[t+1] = x } rm(x,y) } return(list(X=X, accept.prop=naccept/T)) } ################ MH for Cauchy density with Gaussian proposal main.title="MH for Cauchy density with Gaussian proposal" T=1e4 x1=10 log.f = function(x){ -log(1 + x^2) } # par(oma=c(0,0,2,0), mfrow=c(4, 3), mar=c(2, 2, 1, 1), cex=0.6) # for(sd in c(.1, 1, 2.5, 10)){ rq = function(x){ rnorm(1, mean=x, sd=sd) } log.dq = function(x,y){ - log(sd) - (y-x)^2/sd } # ans=log.mhmc(x1, T, log.f, rq, log.dq) # plot(ans$X, type='l', ylim=c(-30,30), main=paste("sd=", sd, ", acceptance ratio=", ans$accept.prop, sep="")) abline(h=c(-2,2), lty=2) lines(cumsum(ans$X)/(1:(T)), type='l', ylim=c(-2, 2), col="blue") # z=seq(-3,3, by=.1) plot(z, dcauchy(z), col='red', type='l', xlim=c(-4,4)) #lines(density(ans$X), col='blue', type='l', xlim=c(-4,4)) hist(ans$X, breaks='FD', prob=TRUE, add=TRUE) # acf(ans$X) #abline(h=mean(X), col='blue') #abline(h=0, col='red') } mhmc.indep = function(x1, T, f, rq, dq){ X = array(NA, T) X[1]=x1 naccept=0 for(t in 1:(T-1)){ y = rq() x = X[t] if( runif(1) <= f(y)*dq(x)/(f(x)*dq(y)) ){ X[t+1] = y naccept = naccept+1 } else { X[t+1] = x } rm(x,y) } return(list(X=X, accept.prop=naccept/T)) } log.mhmc.indep = function(x1, T, log.f, rq, log.dq){ X = array(NA, T) X[1]=x1 naccept=0 for(t in 1:(T-1)){ y = rq() x = X[t] logr = log.f(y) - log.f(x) + log.dq(x) - log.dq(y) if(logr >= 0 || runif(1) <= exp(logr) ){ X[t+1] = y naccept = naccept+1 } else { X[t+1] = x } rm(x,y) } return(list(X=X, accept.prop=naccept/T)) } ################ independent MH for Gaussian(0,1) density with Gaussian proposal main.title="independent MH for Gaussian(0,1) density with Gaussian proposal" T=1e4 log.f = function(x){ -x^2/2 } f = function(x){ dnorm(x) } # par(oma=c(0,0,2,0), mfrow=c(4, 2), mar=c(2, 2, 1, 1), cex=0.6, oma=c(0,0,2,0)) # for(sd in c( .1,.5, 1, 3)){ rq = function(){ rnorm(1, mean=0, sd=sd) } log.dq = function(y){ - log(sd) - (y)^2/(2*sd^2) } dq = function(y){ dnorm(y, sd=sd) } x1=1+rq() # ans=log.mhmc.indep(x1, T, log.f, rq, log.dq) #ans=mhmc.indep(x1, T, f, rq, dq) # plot(ans$X, type='l', ylim=c(-5,5), main=paste("sd=", sd, ", acceptance ratio=", ans$accept.prop, sep="")) abline(h=c(-2,2), lty=2) lines(cumsum(ans$X)/(1:(T)), type='l', ylim=c(-2, 2), col="blue") # z=seq(-3,3, by=.1) plot(z, f(z), col='red', type='l', xlim=c(-3,3)) #lines(density(ans$X), col='blue', type='l', xlim=c(-4,4)) hist(ans$X, breaks='FD', prob=TRUE, add=TRUE) # #acf(ans$X) #abline(h=mean(X), col='blue') #abline(h=0, col='red') } set.seed(1) ################ independent MH for Cauchy(0,1) density with Gaussian proposal main.title="independent MH, target distn = Cauchy(0,1), proposal N(0, sd^2)" T=1e4 log.f = function(x){ -log(1 + x^2) } f = function(x){ dt(x, df=1) } # par(mfrow=c(4, 3), mar=c(2, 2, 1, 1), cex=0.6, oma=c(0,0,2,0)) # for(sd in c( .1,.5, 2, 5)){ rq = function(){ rnorm(1, mean=0, sd=sd) } log.dq = function(y){ - log(sd) - (y)^2/(2*sd^2) } dq = function(y){ dnorm(y, sd=sd) } x1=1+rq() # ans=log.mhmc.indep(x1, T, log.f, rq, log.dq) #ans=mhmc.indep(x1, T, f, rq, dq) # plot(ans$X, type='l', main=paste("sd=", sd, ", acceptance ratio=", ans$accept.prop, sep="")) #, ylim=c(-5,5)) abline(h=c(-2,2), lty=2) lines(cumsum(ans$X)/(1:(T)), type='l', ylim=c(-2, 2), col="blue") # z=seq(-20,20, by=.1) plot(z, f(z), col='red', type='l', xlim=range(z)) abline(v=qnorm(T^(-1)/2)*c(-sd,sd), lty=2) #lines(density(ans$X), col='blue', type='l', xlim=c(-4,4)) hist(ans$X, breaks='FD', prob=TRUE, add=TRUE) # acf(ans$X) #abline(h=mean(X), col='blue') #abline(h=0, col='red') } title(main=main.title, outer=TRUE) ################ independent MH for Cauchy density with Cauchy proposal main.title="independent MH, target distn = Cauchy(0,1), proposal student t_df" T=1e4 x1=4 log.f = function(x){ -log(1 + x^2) } f = function(x){ dt(x,df=1) } # par(mfrow=c(4, 3), mar=c(2, 2, 1, 1), cex=0.6, oma=c(0,0,2,0)) # for(df in c( 3,2,.5,.2)){ rq = function(){ rt(1,df=df ) } dq = function(y){ dt(y, df=df) } log.dq = function(y){ log(dt(y, df=df)) } # ans=log.mhmc.indep(x1, T, log.f, rq, log.dq) #ans=mhmc.indep(x1, T, f, rq, dq) # plot(ans$X, type='l', ylim=c(-50,50), main=paste("df=", df, ", acceptance ratio=", ans$accept.prop, sep="")) abline(h=c(-2,2), lty=2) lines(cumsum(ans$X)/(1:(T)), type='l', ylim=c(-2, 2), col="blue") # z=seq(-10,10, by=.1) plot(z, dcauchy(z), col='red', type='l') #lines(density(ans$X), col='blue', type='l', xlim=c(-4,4)) hist(ans$X, breaks='FD', prob=TRUE, add=TRUE) # acf(ans$X) #abline(h=mean(X), col='blue') #abline(h=0, col='red') } title(main=main.title, outer=TRUE) rw.mhmc = function(x1, T, f, rq){ X = array(NA, T) X[1]=x1 naccept=0 for(t in 1:(T-1)){ y = rq(X[t]) x = X[t] if( ( runif(1) <= f(y)/(f(x)) ) ){ X[t+1] = y naccept = naccept+1 } else { X[t+1] = x } rm(x,y) } return(list(X=X, accept.prop=naccept/T)) } ################ random walk MH for complex density with Gaussian proposal set.seed(1) main.title="random walk MH for complex density with Gaussian proposal" T=1e4 x1=5 #f = function(x){ (exp(-x^2)*(( cos(2*x))^2 + (sin(pi*x))^2*(1+cos(5*x)^3) ) ) } f = function(x){ (exp(-x^2/10)*(sin(x/20)^2)*(1 + (sin(2*x))^3 + (sin(pi*x))^2*(1+cos(5*x)^3) ) ) } # par(mfrow=c(4, 3), mar=c(2, 2, 1, 1), cex=0.6, oma=c(0,0,2,0)) # for(sd in c(.4, 1, 2, 5)){ rq = function(x){ x + rnorm(1,mean=0, sd=sd)} # ans=rw.mhmc(x1, T, f, rq) # plot(ans$X, type='l', ylim=c(-7,7), main=paste("sd=", sd, ", acceptance ratio=", ans$accept.prop, sep="")) z=seq(-7,7, by=.0001) lines(f(z)*T/3, z, col=2) abline(h=c(-2,2), lty=2) lines(cumsum(ans$X)/(1:(T)), type='l', ylim=c(-2, 2), col="blue") # z=seq(-10,10, by=.0001) #plot(z, f(z)/3*5/7, col='red', type='l') plot(z, 10*f(z), col='red', type='l') #lines(density(ans$X), col='blue', type='l', xlim=c(-4,4)) hist(ans$X, breaks='FD', prob=TRUE, add=TRUE) # acf(ans$X) #abline(h=mean(X), col='blue') #abline(h=0, col='red') } title(main=main.title, outer=TRUE) ##### same example, but multiple run of the MC nruns=10 T=1e4 sd=.4 ### this is too low, but good for illustration purposes x1=4 rq = function(x){ x + rnorm(1,mean=0, sd=sd)} # X=array(NA, c(T, nruns)) for(i in 1:nruns){ ans=(rw.mhmc(rt(1, df=1), T, f, rq)) X[,i]=ans$X print(paste("acceptance rate ", ans$accept.prop)) } par(mfrow=c(1,1)) matplot(X[], type="l", lwd=.6, ylim=c(-10,10)) z=seq(-10,10, by=.0001) lines(f(z)*T*30, z, col=1) plot.ts(X, lwd=.5) ############################## MHMC for multidimensional target densities mhmc.rw.multidim = function(dim.pi, x1, T, f, rq){ X = array(NA, c(T,dim.pi)) X[1,]=x1 naccept=0 for(t in 1:(T-1)){ y = rq(X[t,]) x = X[t,] if( ( runif(1) <= f(y)/(f(x)) ) ){ X[t+1,] = y naccept = naccept+1 } else { X[t+1,] = x } rm(x,y) } return(list(X=X, accept.prop=naccept/T)) } ############################## Random Walk MH for bivariate mixture of Gaussians T=1e3 mu1=c(2,2) var1=.3 mu2=c(-2,-2) var2=.4 prob1=.2 library(mvtnorm) f = function(x){ return( prob1*dmvnorm(x, mean=mu1, sigma=diag(rep(var1, 2)) ) + (1-prob1)*dmvnorm(x, mean=mu2, sigma=diag(rep(var2,2))) ) } # par(mfrow=c(2,2), mar=c(2, 2, 1, 1), cex=0.6, oma=c(0,0,2,0)) for( sd in c(.1, .4, 1, 3 ) ){ rq=function(x){ x + rmvnorm(1, mean=c(0,0), sigma=sd^2*diag(2)) } dq=function(x,y){ dmvnorm(y, mean=x, sigma=sd^2*diag(2)) } # ans=mhmc.rw.multidim(2, mu1, T, f, rq) x=seq(-4,4, len=1e2) z = array(NA, c(length(x), length(x))) for(i in 1:length(x)){ for(j in 1:length(x)){ z[i,j] = log(f(c(x[i], x[j])), base=10) }} contour(x,x, z, method = "edge", nlevels=20) lines(ans$X, pch=20, type='o') title(main=paste("sd = ", sd , ", acceptance ratio=", ans$accept.prop, sep="")) } ## ############################## Random Walk MH for mixture of multivariate Gaussians (p=3) ## T=1e4 ## mu1=c(2,2,2) ## var1=.2 ## mu2=c(-2,-2, -2) ## var2=.2 ## mu3=c(4,4, 4) ## var3=.1 ## p1=.2 ## p2=.3 ## ## library(mvtnorm) ## f = function(x){ ## return( p1*dmvnorm(x, mean=mu1, sigma=diag(rep(var1, 3)) ) + p2*dmvnorm(x, mean=mu2, sigma=diag(rep(var2,3))) ## + (1-p1-p2)*dmvnorm(x, mean=mu3, sigma=diag(rep(var3,3))) ) ## } ## #par(mfrow=c(2,2), mar=c(2, 2, 1, 1), cex=0.6, oma=c(0,0,2,0)) ## sd=.4 ## rq=function(x){ ## rmvnorm(1, mean=x, sigma=sd^2*diag(3)) ### Gaussian proposal ## #x + rmvt(1, sigma=sd^2*diag(3), df=1) ### multivariate t proposal ## } ## # ## ans=mhmc.rw.multidim(3, mu3, T, f, rq) ## plot.ts(ans$X) ## ## acf(ans$X) ## ## library(lattice) ## pairs(ans$X) ## print(ans$a) ## ## par(mfrow=c(1,1)) ## library(scatterplot3d) ## scatterplot3d(ans$X) ############################## Random Walk MH for mixture of multivariate Gaussians (p=3) T=3e4 mu1=c(2,2,2) var1=.2 mu2=c(-2,-2, -2) var2=.2 mu3=c(4,4, 4) var3=.1 p1=.2 p2=.3 library(mvtnorm) f = function(x){ return( p1*dmvnorm(x, mean=mu1, sigma=diag(rep(var1, 3)) ) + p2*dmvnorm(x, mean=mu2, sigma=diag(rep(var2,3))) + (1-p1-p2)*dmvnorm(x, mean=mu3, sigma=diag(rep(var3,3))) ) } #par(mfrow=c(2,2), mar=c(2, 2, 1, 1), cex=0.6, oma=c(0,0,2,0)) sd=1 rq=function(x){ # rmvnorm(1, mean=x, sigma=sd^2*diag(3)) ### Gaussian proposal x + rmvt(1, sigma=sd^2*diag(3), df=1) ### multivariate t proposal } # ans=mhmc.rw.multidim(3, mu3, T, f, rq) plot.ts(ans$X) acf(ans$X) library(lattice) pairs(ans$X) print(ans$a) par(mfrow=c(1,1)) library(scatterplot3d) scatterplot3d(ans$X) dev.off()