rbvnorm<-function(x,y,m1=0,m2=0,s1=1,s2=1,rho=.5){
exp(-((x-m1)^2/s1^2-2*rho*((x-m1)/s1)*((y-m2)/s2)+(y-m2)^2/s2^2)/(2*(1-rho^2)))
}