metropolis2=function( n=1000 , rho=0.0 , start1=0 , start2=0 , tau=0.25 , sigma2=0.1 ) { x1 <- c(1:n) x2 <- c(1:n) x1[1] <- start1 x2[1] <- start2 scale <- -0.5 / ( sigma2* ( 1 - rho^2 ) ) breaks <- c( -20:(1+20) )/4 par( mfrow=c( 2,2 ) ) for ( i in 2:n ) { try <- x1[i-1] + rnorm( 1 , 0 , tau ) temp <- exp( scale * ( try^2 - 2*rho*try*x2[i-1] +x2[i-1]^2) ) + exp( scale * ( (try-1)^2 - 2*rho*(try-1)*(x2[i-1]-1) + (x2[i-1]-1)^2 ) ) temp <- temp / ( exp( scale * ( x1[i-1]^2 - 2*rho*x1[i-1]*x2[i-1] + x2[i-1]^2 ) ) + exp( scale * ( (x1[i-1]-1)^2 - 2*rho*(x1[i-1]-1)*(x2[i-1]-1) + (x2[i-1]-1)^2 ) ) ) if ( runif(1) < temp ) { x1[i] <- try } else { x1[i] <- x1[i-1] } try <- x2[i-1] + rnorm(1,0,tau ) temp <- exp( scale * ( try^2 - 2*rho*try*x1[i] + x1[i]^2) ) + exp( scale * ( (try-1)^2 - 2*rho*(try-1)*(x1[i]-1) + (x1[i]-1)^2 ) ) temp <- temp / ( exp( scale * ( x2[i-1]^2 - 2*rho*x1[i]*x2[i-1] + x1[i]^2 ) ) + exp( scale * ( (x2[i-1]-1)^2 - 2*rho*(x1[i]-1)*(x2[i-1]-1) +(x1[i]-1)^2 ) ) ) if ( runif(1) < temp ) { x2[i] <- try } else { x2[i] <- x2[i-1] } hist( x1[1:i] , breaks , main="x1" , xlab="x1") hist( x2[1:i] , breaks , main="x2" , xlab="x2" ) plot( x1[1:i] , x2[1:i] , type="p" , pch="." ,xlab="x1",ylab="x2") plot( x1[1:i] , x2[1:i] , type="l" ,col="red",xlab="x1",ylab="x2") } }