Correlated bivariate normal distribution: finding percentage of of data which is 2 standard deviations above the mean?

Question

Background

Background: in an essay of mine, I point out that if a selection process (like higher education) requires successful applicants to be above the mean on 2 different variables, there will necessarily be fewer successful applicants than when successful applicants just had to be above the mean on just 1 variable. (The 2 different variables in this case are the Big Five personality factor of Conscientiousness, and IQ.) This reduction is most dramatic when the 2 variables have 0 correlation, but even a large correlation will still result in many applicants being filtered out. How many are filtered out, exactly?

Simple questions

Well, if the filter is for 2 standard deviations above the mean and the variables are correlated with 1, then 2.3% of the population will pass; if the variables are uncorrelated with 0, then 2.3% * 2.3% or 5.29e-4% of the population will pass.

Correlated

But what about intermediate values? For example, the psychology literature has reported a correlation of -0.21 between Conscientiousness & IQ.

I consulted Wikipedia on bivariate normal distributions, but I didn't understand much of it. The closest I found was sum of correlated normal random variables, but in this case what I want is closer to a min function.

Simulation

I wa able to work up a R simulation to see how that worked, and it seemed in line with my intuitions:

install.packages("fMultivar")
library ("fMultivar")

x <- rnorm2d(10000000, rho=0.5)

xgreater <- length(subset(x, x[,1] > mean(x[,1])+2*sd(x[,1])))
xandygreater <- length(subset(x, x[,1] > mean(x[,1])+2*sd(x[,1]) & x[,2] > mean(x[,2])+2*sd(x[,2])))

c(xgreater, xandygreater); c(xgreater / length(x), xandygreater / length(x), xgreater / xandygreater) * 100

# example results for different values of 'rho='
0.1
[1] 454,664  17,570
[1] 2.273e+00 8.785e-02 2.588e+03

0.2
[1] 458,284  82,552
[1]   2.2914   0.4128 555.1458

0.5
[1] 454,484  80,872
[1]   2.2724   0.4044 561.9794

0.9
[1] 455,242 267,912
[1]   2.276   1.340 169.922

0.95
[1] 455,162 321,024
[1]   2.276   1.605 141.784

0.99
[1] 455,260 394,448
[1]   2.276   1.972 115.417

Exact pdf calculation?

I really was hoping for more of a precise analytic solution, so some more searching eventually turned up a paper, "Exact Distribution of the Max/Min of Two Gaussian Random Variables", which gives a definition for the min of 2 correlated normal variables. This seems to be what I want; top of pg1, second column:

...where $\phi(.)$ and $\Phi(.)$ are, respectively, the pdf and the cumulative distribution function (cdf) of the standard normal distribution. It is known that the pdf of $Y = \min(X_1, X_2)$ is $f(y) = f_1(y) + f_2(y)$, where

(3) $f_1(y) = \frac{1}{\sigma_1} \phi (\frac{y-\mu_1}{\sigma_1}) \times \Phi (\frac{p(y - \mu_1)}{\sigma_1 \sqrt{1 - p^2}} - \frac{y - \mu_2}{\sigma_2 \sqrt{1 - p^2}})$ (4) $f_2(y) = \frac{1}{\sigma_2} \phi (\frac{y-\mu_2}{\sigma_2}) \times \Phi (\frac{p(y - \mu_2)}{\sigma_2 \sqrt{1 - p^2}} - \frac{y - \mu_1}{\sigma_1 \sqrt{1 - p^2}})$

They give an R implementation on pg6 (first column); it seems to have a pnorm typo, but I fixed that. Once it was working, I tried generating a slightly (0.1) correlated bivariate distribution, which look OK:

fmin<-function (y,mu1,mu2,sigma1,sigma2,rho)
     {t1<-dnorm(y,mean=mu1,sd=sigma1)
     tt<-rho*(y-mu1)/(sigma1*sqrt(1-rho*rho))
     tt<-tt-(y-mu2)/(sigma2*sqrt(1-rho*rho))
     t1<-t1*pnorm(tt)
     t2<-dnorm(y,mean=mu2,sd=sigma2)
     tt<-rho*(y-mu2)/(sigma2*sqrt(1-rho*rho))
     tt<-tt-(y-mu1)/(sigma1*sqrt(1-rho*rho))
     t2<-t2*pnorm(tt)
     return(t1+t2)}
fmin(c(1:200),100,100,15,15,0.1)
  [1] 1.849e-11 2.864e-11 4.418e-11 6.784e-11 1.037e-10 1.578e-10 2.392e-10 3.608e-10 5.418e-10
 [10] 8.101e-10 1.206e-09 1.787e-09 2.636e-09 3.872e-09 5.663e-09 8.243e-09 1.195e-08 1.724e-08
 [19] 2.476e-08 3.542e-08 5.043e-08 7.148e-08 1.009e-07 1.417e-07 1.982e-07 2.760e-07 3.827e-07
 [28] 5.282e-07 7.257e-07 9.928e-07 1.352e-06 1.833e-06 2.475e-06 3.326e-06 4.449e-06 5.926e-06
 [37] 7.859e-06 1.037e-05 1.364e-05 1.784e-05 2.324e-05 3.014e-05 3.891e-05 5.002e-05 6.401e-05
 [46] 8.154e-05 1.034e-04 1.306e-04 1.641e-04 2.054e-04 2.558e-04 3.173e-04 3.917e-04 4.814e-04
 [55] 5.889e-04 7.173e-04 8.696e-04 1.049e-03 1.261e-03 1.507e-03 1.794e-03 2.125e-03 2.506e-03
 [64] 2.941e-03 3.435e-03 3.993e-03 4.620e-03 5.318e-03 6.093e-03 6.945e-03 7.878e-03 8.890e-03
 [73] 9.982e-03 1.115e-02 1.239e-02 1.370e-02 1.506e-02 1.647e-02 1.791e-02 1.938e-02 2.084e-02
 [82] 2.230e-02 2.371e-02 2.508e-02 2.636e-02 2.755e-02 2.863e-02 2.956e-02 3.034e-02 3.095e-02
 [91] 3.138e-02 3.162e-02 3.165e-02 3.149e-02 3.112e-02 3.056e-02 2.981e-02 2.889e-02 2.781e-02
[100] 2.660e-02 2.526e-02 2.383e-02 2.233e-02 2.077e-02 1.920e-02 1.762e-02 1.605e-02 1.452e-02
[109] 1.305e-02 1.164e-02 1.031e-02 9.063e-03 7.912e-03 6.857e-03 5.899e-03 5.039e-03 4.272e-03
[118] 3.595e-03 3.004e-03 2.491e-03 2.050e-03 1.675e-03 1.358e-03 1.093e-03 8.732e-04 6.923e-04
[127] 5.447e-04 4.254e-04 3.297e-04 2.535e-04 1.935e-04 1.466e-04 1.102e-04 8.220e-05 6.085e-05
[136] 4.470e-05 3.258e-05 2.357e-05 1.692e-05 1.205e-05 8.517e-06 5.973e-06 4.157e-06 2.870e-06
[145] 1.966e-06 1.337e-06 9.018e-07 6.036e-07 4.008e-07 2.641e-07 1.727e-07 1.120e-07 7.210e-08
[154] 4.604e-08 2.917e-08 1.834e-08 1.144e-08 7.078e-09 4.346e-09 2.647e-09 1.600e-09 9.594e-10
[163] 5.708e-10 3.369e-10 1.973e-10 1.146e-10 6.607e-11 3.778e-11 2.143e-11 1.207e-11 6.737e-12
[172] 3.733e-12 2.052e-12 1.119e-12 6.052e-13 3.248e-13 1.730e-13 9.137e-14 4.788e-14 2.489e-14
[181] 1.284e-14 6.571e-15 3.336e-15 1.680e-15 8.393e-16 4.160e-16 2.046e-16 9.979e-17 4.830e-17
[190] 2.319e-17 1.104e-17 5.218e-18 2.446e-18 1.137e-18 5.248e-19 2.402e-19 1.090e-19 4.911e-20
[199] 2.194e-20 9.727e-21

Now, I understand the PDF to be "a function that describes the relative likelihood for this random variable to take on a given value. The probability for the random variable to fall within a particular region is given by the integral of this variable’s density over the region". So I suppose I should sum up every point in the pdf >130 (since 130 is 2 standard deviations up, by construction when I specified SD=15) and that's my probability that a random deviate will be min(130,130). What's the total probability someone will be over 130 on both variables? I think that would be:

sum(fmin(c(1:200),100,100,15,15,0.1)[130:200])
[1] 0.001004

If I increase the r to 0.9, the result is 0.01455 which is satisfyingly larger.

A sanity check - as the correlation goes to 1.0, there should be no decrease. So we do the same question for a single normal distribution defined the same way:

sum(dnorm(c(1:200), 100, 15)[130:200])
[1] 0.02459

# the function blows NaN chunks on 1.0, so we'll try a lot of 9s:
sum(fmin(c(1:200),100,100,15,15,0.9999999999)[130:200])
[1] 0.02459

I guess that works too.

Problems

So my questions are:

1. Is my R simulation right?
2. Is my version and use of fmin right?
3. Is there some more direct, possibly even pen-and-paper, avenue of calculating the answer to my original question?

Answer 1

If I have understood your question correctly, you want the CDF of the bivariate normal distribution. That is, for the standardized case:

$$ \Phi(\mathrm{\pmb{b}},\rho) = \frac{1}{2\pi\sqrt{1-\rho^{2}}}\int_{-\infty}^{b_{1}}{\int_{-\infty}^{b_{2}}}\exp\left[{-(x^{2}-2\rho xy+y^{2}})/(2(1-\rho^{2})\right]\mathrm{d}y\mathrm{d}x $$

This has no closed form solution and must be integrated numerically. With modern software, this is quite trivial.

Here is an example in R with perfectly correlated normal distributions (i.e. $\rho = 1$) with $\pmb{\mu}=(100, 40)^\intercal$ and covariance matrix $\Sigma = \begin{bmatrix} 225 & 75 \\ 75 & 25 \end{bmatrix}$. We calculate the probabiltiy of both variables being above $2$ SD:

library(mvtnorm)
library(MASS)

corr.mat <- matrix(c(1, 1, 1, 1), 2, 2, byrow = TRUE) # correlations
sd.mat <- matrix(c(15, 0, 0, 5), 2, 2, byrow = TRUE) # standard deviations

cov.mat <- sd.mat %*% corr.mat %*% sd.mat # covariance matrix

mu <- c(100, 40) # means

pmvnorm(lower = mu + 2*diag(sd.mat), upper = Inf, mean = mu, sigma = cov.mat)

[1] 0.02275013
attr(,"error")
[1] 2e-16
attr(,"msg")
[1] "Normal Completion"

As you rightly said: When they are perfectly correlated, the probability is about 2.3%.

What about a correlation of -0.21?

corr.mat <- matrix(c(1, -0.21, -0.21, 1), 2, 2, byrow = TRUE) # correlations
sd.mat <- matrix(c(15, 0, 0, 5), 2, 2, byrow = TRUE) # standard deviations

cov.mat <- sd.mat %*% corr.mat %*% sd.mat # covariance matrix

mu <- c(100, 40) # means

pmvnorm(lower = mu + 2*diag(sd.mat), upper = Inf, mean = mu, sigma = cov.mat)

[1] 0.0001228342
attr(,"error")
[1] 1e-15
attr(,"msg")
[1] "Normal Completion"

The probability is much lower, namely 0.0001228342.

We can verify our calculations by simulation. For the example above:

set.seed(21)
dat <- rmvnorm(1e7, mean = c(100, 40), sigma = cov.mat)

sum(dat[, 1] > (100+2*15) & dat[, 2] > (40+2*5))/dim(dat)[1]

[1] 0.0001261

This is very close to the result from numerical integration.

These calculations can easily be extended to multivariate normal distributions.

Answer 2

Good question. Your intuitions and approaches are right. There is, however, no simple analytic formula for the joint probability you're after. One could write a short program in R to compute the probability but it would necessarily use numerical integration.

Answer 3

Here is a useful explanation for #3: Multivariate Normal distribution

See 'bivariate normal distribution' in that section you can see the pdf for a bivariate normal distribution with the correlation coefficient. You can integrate the pdf between -infinity to $\mu$ + 2*$\sigma$. Or perhaps this is enough info to find an equation for the CDF, which can be evaluated at $\mu$ + 2*$\sigma$ to give you what you need.