Load packages

library(tidyverse)
library(rockchalk)
library(MASS)

Defining the data

myR <- lazyCor(X = 0.7, d = 2)
mySD <- c(1, 1)
myCov <- lazyCov(Rho = myR, Sd = mySD)

myR
#>      [,1] [,2]
#> [1,]  1.0  0.7
#> [2,]  0.7  1.0
mySD
#> [1] 1 1
myCov
#>      [,1] [,2]
#> [1,]  1.0  0.7
#> [2,]  0.7  1.0

Drawing from the multivariate normal

Let’s draw 1000 cases. Met \(\mu\) be zero.

twogauss <- mvrnorm(n = 1000, 
                    mu = c(0,0), 
                    Sigma = myCov) %>% as_tibble()

Estimate densities

mydens <- kde2d(twogauss$V1, twogauss$V2)

Plot it

plotly::plot_ly(z = mydens$z, 
                x=mydens$x, 
                y=mydens$y, 
                type = "surface")

Bonus: Uncorrelated bivariate Gaussian

This is simpler.

First define a grid for the value range of \(x\) and \(y\).

x <- y <- seq(-6, 6, .1)

Define the density function (here Gaussian)

Say, \(\sigma=2\).

s <- 2
f <- function(x,y){exp(-(y^2+x^2)/(2*s^2))}

Compute density for each x-y combination. That is, apply function f on each x-y combination. outer will do that for us.

z <- outer(x,y,FUN = f)

Plot it:

plotly::plot_ly(z = z, x=x, y=y, type = "surface")

Simpler (and maybe more appropriate) alternative

Plotting in 3D may appear cool, but a simple heat map may be more appropriate.

Check this out:

First, let’s come up with some noisy data:

myR <- lazyCor(X = 0.3, d = 2)
mySD <- c(1, 1)
myCov <- lazyCov(Rho = myR, Sd = mySD)

twogauss <- mvrnorm(n = 1000, 
                    mu = c(0,0), 
                    Sigma = myCov) %>% as_tibble()

head(twogauss)
#> # A tibble: 6 x 2
#>        V1      V2
#>     <dbl>   <dbl>
#> 1  0.0412  0.975 
#> 2 -0.694  -0.735 
#> 3  0.761   0.0920
#> 4 -0.786   0.0233
#> 5  0.708   0.258 
#> 6 -0.923  -1.09

myplot1 <- ggplot(twogauss, aes(x = V1, y = V2))

myplot1 + geom_bin2d() + scale_fill_viridis_c()