# 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?
• I take it "correlated with 1" and "uncorrelated with 0" mean r = 1.0 and r = 0.0, respectively. – rolando2 Jan 17 '17 at 0:35

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.

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.

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.

• I'm not sure that's helpful because the PDF of the min of 2 variables isn't the same thing as the PDF of 2 variables. – gwern Jan 18 '17 at 1:52