#pdf for bivariate normal distribution N(0,S) 
# where S = (  1  rho )
#           ( rho  1  )
# used in plotting contours in gibbs.update
f=function(x,y,rho) {
exp(-(x^2 - 2*rho * x*y +y^2)/(2*(1-rho^2)))/(2*pi*sqrt(1-rho^2))
}