normal.mh=function(mu,sigma, sigmaq=c(sqrt(.6),sqrt(.4)), istart=100, istop = 2000) { 
# I wasn't sure I wanted as little flexibility in sigmaq as I show here.  I considered using a default matrix and 
# computing the cholesky decomposition.  I also thought about generate a bivariate normal using mvrnom in 
# MASS.  Both choices seem to defeat the purpose of a simple choice for q, so I stuck with the choice you 
# see here.  There's nothing magic about using .6 and .4 for variances; it's just the same values I had used in
# class notes.  
isig=solve(sigma)
iestop=istop+istart
x <- matrix(0,2*iestop,nrow=2)
x[,1]=mu
for(i in 2:iestop) { 
y=c(rnorm(1,x[1,(i-1)],sigmaq[1]),rnorm(1,x[2,(i-1)],sigmaq[2]))
sigy=t(y-mu)%*%isig%*%(y-mu)
sigx=t(x[,(i-1)]-mu)%*%isig%*%(x[,(i-1)]-mu)
alph=min(exp(-.5*(sigy-sigx)),1)
if(runif(1)<alph) x[,i]=y else x[,i]=x[,(i-1)]
} 
out=x[,(istart+1):iestop]
rownames(x)=c("X","Y")
colnames(out)=1:istop
return(out) 
} 

# This is the target distribution we discussed in class
mu=c(1,2)
s2=matrix(c(1,.9,.9,1),ncol=2)
x=normal.mh(mu,s2)
# Plot with contour overlay
plot(x[1,],x[2,],xlab="X",ylab="Y")
cxseq=seq(-2,4,by=.1)
cyseq=seq(-1,5,by=.1)
xmat=cxseq%o%rep(1,length(cxseq))
ymat=rep(1,length(cyseq))%o%cyseq
# I couldn't figure out a good way to stack individual elements of xmat and ymat in vectors so that
# could compute the density over a grid using vector/matrix notation.  So I used the element-wise
# approach below instead.
is2=solve(s2)
sd1=sqrt(s2[1,1])
sd2=sqrt(s2[2,2])
rho=s2[1,2]/(sd1*sd2)
xmat=(xmat-mu[1])/sd1
ymat=(ymat-mu[2])/sd2
zmat=exp(-.5*(xmat^2+ymat^2-2*rho*xmat*ymat ))/(2*pi*sd1*sd1*sqrt((1-rho)^2))
# add=T overlays the contours on the current plot
contour(cxseq,cyseq,zmat,add=T)
# Let's look at the individual histograms
win.graph()
par(mfrow=c(1,2))
hist(x[1,],freq=F,xlab="",main=paste("Histogram of ","X"))
lines(cxseq,dnorm(cxseq,mu[1],sd1))
hist(x[2,],freq=F,xlab="",main=paste("Histogram of", "Y"))
lines(cyseq,dnorm(cyseq,mu[2],sd2))