#metroplis example: bivariate normal density with bivariate normal jumping kernel metropolis_function(theta=c(1,1),mu=c(0,0),Sigma=diag(1,2)){ thetastar_simulate.multnorm(theta,Sigma) Prec_solve(Sigma) ratio _ exp(-1/2*(t(thetastar-mu )%*%Prec%*%(thetastar-mu) - t(theta-mu)%*%Prec%*%(theta-mu))) u _ runif(1) if(u <= ratio) {theta_thetastar} test_ (u <= ratio) return(theta) } iterate_function(theta0=c(0,0),NN){ theta_matrix(NA,2,NN) test_NULL theta[,1]_theta0 for(n in 2:(NN-1)){ theta[,n]_metropolis(theta=theta[,n-1],mu=c(0,0),Sigma=diag(1,2))} return(theta) } plotmcmc_function(w=0,NN){ if (w == 1) postscript("/home/biostats/fdominic/course/mcmc.ps") par(oma=c(0,0,2,0)) I_iterate(theta0=c(0,0),NN) I1_iterate(theta0=c(3,3),NN) I2_iterate(theta0=c(3,-3),NN) I3_iterate(theta0=c(-3,3),NN) I4_iterate(theta0=c(-3,-3),NN) plot(I[1,],I[2,],type="l",xlim=c(-4,4), ylim=c(-4,4),xlab="theta1",ylab="theta2") points(0,0,pch=0,cex=4) lines(I1[1,],I1[2,],lty=2) points(3,3,pch=0,cex=4) lines(I2[1,],I2[2,],lty=2) points(3,-3,pch=0,cex=4) lines(I3[1,],I3[2,],lty=2) points(-3,3,pch=0,cex=4) lines(I4[1,],I4[2,],lty=2) points(-3,-3,pch=0,cex=4) par(oma=c(0,0,0,0)) if (w == 1) dev.off() }