R function to plot contours of a bivariate normal distribution

plot_bivar_norm.R

plot_bivar_norm <- function(mx, my, sx, sy, r){
  
  require(ggplot2)
  
  logbinorm <- function (x, y, par) {
    m <- par$m
    v <- par$v
    zx <- (x - m[1]) / sqrt(v[1, 1])
    zy <- (y - m[2]) / sqrt(v[2, 2])
    r <- v[1, 2] / sqrt(v[1, 1] * v[2, 2])
    return(-0.5 / (1 - r ^ 2) * (zx ^ 2 - 
              2 * r * zx * zy + zy ^ 2))
  }

  set_up <- function(mx, my, sx, sy, r){
    v <- matrix(c(sx ^ 2, r * sx * sy, r * sx * sy, sy ^ 2),
            2, 2)
    d <- expand.grid(X = seq(mx - 4 * sx, mx + 4 * sx, 
                         length.out=50),
                 Y = seq(my - 4 * sy, my + 4 * sy, 
                         length.out=50))
    d$Z <- logbinorm(d$X, d$Y, list(m = c(mx, my), v = v))
   d
  }
  
  d1 <- set_up(mx, my, sx, sy, r)
  
  ggplot(d1, aes(x=X, y=Y, z=Z)) +
    stat_contour(breaks = c(-6.9, -4.6, -2.3)) + 
    coord_fixed() 
}