Simulate a bivariate normal

Cut and past the following code to define a function called rbivariate that generates a random sample from a bivariate normal pdf.


rbivariate <- function(mean.x = 70, sd.x=3, mean.y=162, sd.y=14, r=.50, iter=100) {
   z1 <- rnorm(iter)
   z2 <- rnorm(iter)
   x <- sqrt(1-r^2)*sd.x*z1 + r*sd.x*z2 + mean.x
   y <- sd.y*z2 + mean.y
   return(list(x,y))
}

Now run the following code to see that it works.


data <- rbivariate(iter=1000)
mean(data[[1]])
sd(data[[1]])
mean(data[[2]])
sd(data[[2]])
plot(data[[1]],data[[2]])