Using R: Calculating Probability for a Bivariate Normal Random Variable

This a useful topic to know. Ultimately we are going to use this for calculating power and sample size for the multivariate sign test. Help this channel to remain great! Donating to Patreon or Paypal can do this!

[statisticsmatt]

###
### Calculating probability for a bivariate normal random
### variable.
###


library("MASS")
rm(list=ls())

### Bivariate Normal Density
h <- function(x, y, mx, my, sx, sy, p) { (1/(2* pi * - 2*p*((x-mx)/sx)*((y-my)/sy) + (y-my)^2/sy^2})}

### Function for calculating probabilities
bvn.prob<-function(mx, my, sx, sy, p, x.llim, x.ulim, y.llim, y.ulim) {
integrate(function(y){ sapply(y, function(y){
integrate(function(x) h(x, y, mx, my, sx, sy, p), x.llim,x.ulim)$value})},y.llim,y.ulim)$value }

### Picking values for a simulation
0.2
### Calculating the Variance/Covariance Matrix
(Sigma <- matrix(c(sdx^2, sdx*sdy*rho, sdx*sdy*rho, sdy^2), 2, 2))
### Mean vector
Mu<- c(mux, muy)

### Sample Size of One million to estimate probabilities for
### each quadrant.

y<-mvrnorm(n, Mu, Sigma)
tmp1<-sum((y[, 1]>0)*(y[, 2]>0))/n
tmp2<-sum((y[, 1]<0)*(y[, 2]>0))/n
tmp3<-sum((y[, 1]<0)*(y[, 2]<0))/n
tmp4<-sum((y[, 1]>0)*(y[, 2]<0))/n

cbind(
Simulation=c(tmp1, tmp2, tmp3, tmp4),
Calculation=c(
bvn.prob(mux, muy, sdx, sdy, rho, 0, Inf, 0, Inf),
bvn.prob(mux, muy, sdx, sdy, rho, -Inf, 0, 0, Inf),
bvn.prob(mux, muy, sdx, sdy, rho, -Inf, 0, -Inf, 0),
bvn.prob(mux, muy, sdx, sdy, rho, 0, Inf, -Inf, 0)))

### Integrate out one variable and then compare with
### univarite built in function
### (1)
pnorm(1.2)
bvn.prob(2, 0, 5, 1, .4, -Inf, Inf, -Inf, 1.2)
bvn.prob(0, 0, 1, 1, 0, -Inf, 1.2, -Inf, Inf)
### (2)
pnorm(1.5)-pnorm(-.5)
bvn.prob(2, 0, 5, 1, .4, -Inf, Inf, -.5, 1.5)

### Plot the density.
x<-seq(-3, 3, .1)
y<-seq(-3, 3, .1)
f<-function(x, y) h(x, y, 0, 0, 1, 1, -.7)
res<-persp(x, y, outer(x, y, f), theta = 45, phi = 20)
lines (trans3d(x, y = -1, z = f(x, -1), pmat = res), col = 2)