In a testing, ranking, or selection scenario, we have samples of size n where a measurement is correlated 0<r<1 with some second variable of interest; they are bivariate normally distributed. We'd like to select the maximum. If we pick the maximum on the first variable, we have a better chance of picking the maximum on the second variable as well, but of course, with less than 100% probability. What exactly is the probability of selecting the true maximum for a given n and r?

This can be easily simulated (R):

max1IsMax2 <- function(n,r) { mean(replicate(100000, {
                sample <- mvrnorm(n, mu=c(0,0), Sigma=matrix(c(1, r, r, 1), nrow=2))
                })) }

sapply(2:10, function(n) { max1IsMax2(n, r) })
# [1] 0.66925 0.53466 0.46197 0.41044 0.37627 0.34344 0.32760 0.30732 0.29484

But I have failed to figure out any way to calculate it exactly.

I spent a few hours trying to read through the order statistics literature, without little success (most of the discussions of bivariate or multivariate order statistics don't seem to be considering a scenario like this, or I completely failed to understand them) and also following one idea I had: given the expectation on the first variable (the expected maximum of n), the expected second variable is regressed to the mean by r, eg

exactMax <- function(n, mean=0, sd=1) { expect.max.ostat(n, para=vec2par(c(mean, sd), type="nor"), cdf=cdfnor, pdf=pdfnor) }
# [1] 1.53875273
exactMax(10) * r
# [1] 0.769376365

The winner on the test is +1.53SD, but regress to the mean down to +0.76SD on the true latent variable. For it to be the maximum on the second variable as well, that is another way of saying that all the other 9 datapoints in the sample must all be smaller than it. Each datapoint has a certain p of being greater than +0.76SD, so the probability of p failing 9 times is then simply $(1-p)^9$. So to estimate the repeat probability, you'd simply take the expected maximum of n, deflate by r, calculate the CDF for that, and raise it to n-1. p is easily estimated by simulation or exactly for the normal distribution:

1 - mean(rnorm(1000000) > (exactMax(10) * r))
# [1] 0.77932
pnorm(exactMax(10) * r)
# [1] 0.779165043

Except... $0.779^9=0.10$, which is way off from 0.29, and this is true of all the other values I've tried. I wondered if I was specifying the SDs wrong for the second variable, but looking at the simulation, it is indeed $N(0,1)$ as it should be, the maximum on variable 1 does indeed regress by r in variable 2, and I can't find any formula which gives max-probabilities consistent with the simulation estimates even fiddling around with the SDs.

So what am I doing wrong here and what is the right answer? And are there any generalizations discussing how much slippage there is for given n and r?

  • $\begingroup$ Is it correct if i rewrite your problem as follows? Given a sample of size $n$ from a $Normal(\mu,\Sigma)$ with known parameters. Let's call each sample $S_i$ and let's write $S_i = (X_i, Y_i)$. We want to find the probability that $P(X_i \ge X_j | Y_i \ge Y_j)$ for all $j$. $\endgroup$ Commented Sep 15, 2019 at 8:01
  • $\begingroup$ I think you could generalize it to ask about each order statistic, not just the max/min, yes. I am primarily interested in the maximum because most real-world problems take the form of either selecting a threshold (which I already have exact solutions for), or the single biggest (max); relatively few natural problems would be 'select the fifth', eg, and I'd rather not bite off more than I can chew. :) $\endgroup$
    – gwern
    Commented Sep 15, 2019 at 13:59
  • $\begingroup$ Some progress: apparently for n = 2, it's $\frac{1}{2} + \frac{1}{\pi}\arcsin(r)$. And Sibuya 1960 proves that asymptotically in n when r < 1, the two distributions are independent (ie the probability of the max in one being the max in the other just approaches 1 / n - so the correlation is basically a constant factor gain). And Dos Anjos et al 2005 might help but I can't understand it. $\endgroup$
    – gwern
    Commented Dec 24, 2019 at 15:45
  • 1
    $\begingroup$ Finally, The Asymptotic Theory of Extreme Order Statistics, Galambos 1987 might, from some references, cover exact formulas or approximations good for small n, at least, but I can't find any digital copies or books for sale to scan, so I've set up an alert on Abebooks in case any copies surface. $\endgroup$
    – gwern
    Commented Dec 24, 2019 at 15:49

1 Answer 1


There is an analytic answer for a slightly different distribution, the Ali-Mikhail-Haq copula.

If $0\le r \le \frac12$, we can choose this copula to have

  • the same standard normal distribution of $X$'s as the bivariate normal
  • the same standard normal distribution of $Y$'s as the bivariate normal
  • the same Kendall's tau measure of correlation between the two variables.

First we choose the parameter $\theta$ to get Kendall's tau to agree for the two distributions: $$1 - \frac{2((1-\theta)^2\log(1-\theta)+\theta)}{3\theta^2}=\frac{2}{\pi}\arcsin(r)$$

Then we can calculate the probability of $X$ and $Y$ being maximized at the same observation:

  • In the limiting case where $r=\frac{1}{2}$, we have $\theta=1$, and the probability is $\frac{2}{1+n}$

  • In the limiting case of high $n$, the probability tends to $\frac{1+\theta}{n}$ and more precisely $\lim_{n\rightarrow \infty}np(n,\theta)=1+\theta$.

  • In general, the probability is $$p(n,\theta)=t\frac{\, _2F_1\left(1,n+1;n+2;t\right)}{n+1} - 2nt^2\frac{\, _2F_1\left(1,n+2;n+3;t\right)}{(n+1)(n+2)} -\frac{2nt}{(n+1)^2}+\frac{1}{n}$$ where $t=\theta/(\theta-1)$ and $\,_2F_1$ is the hypergeometric function.

The following plot of pdfs shows the closeness of this approximation: the orange is the AMH copula with $\theta=1$, and the blue is the standard bivariate normal with $r=\frac{1}{2}$.

enter image description here

The copula is defined by the formula $$P_{copula}[X\le u,\ Y\le v]=\frac{\Phi(u)\Phi(v)}{1-\theta(1-\Phi(u))(1-\Phi(v))}$$

The advantage of such a copula for order statistics is that we can also do the analysis on the simpler distribution where $X$ and $Y$ are uniform variables between 0 and 1. Then $$P[X\le x,\ Y\le y]=F[x,y]=\frac{xy}{1-\theta(1-x)(1-y)}$$ The pdf for this distribution is: $$f(x,y)=\frac{1-\theta(2-x-y-xy)+\theta^2(1-x-y+xy)}{(1+\theta(1-x)(1-y))^3}$$

If there are $n$ samples from the distribution, the pdf for the maximal $X$ is just $$nx^{n-1}$$ If the maximum of the $X$'s is $x$, then the maximal $Y$ occurs at the same observation with probability $$q(x)=\int_0^1 \frac{F[x,y]^{n-1}}{F[x,1]^{n-1}}f(x,y) dy = \frac{1+n-\theta+n\theta(2x-1)}{n(1+n)(1+\theta x-\theta)}$$

So the overall probability that the maximal $X$ and maximal $Y$ occur at the same observation is $\int_0^1 nx^{n-1}q(x)dx$, which gives the expression for $p(n,\theta)$ at the beginning.

  • $\begingroup$ I tried implementing your formula to compare with the simulation, and I think there's something I'm misunderstanding because the results disagree with the simulation results in my question: pastebin.com/8QKH6Aqm In particular, I don't understand your r to theta transformation: you say the limiting case of r=0.5 should give theta=1, but (2/pi)*asin(0.5) doesn't equal 1, it equals 0.333? $\endgroup$
    – gwern
    Commented Mar 11, 2020 at 16:47
  • $\begingroup$ When $r=1/2$, the right hand side of that equation is indeed $1/3$, and when you solve for the left hand side of the equation equal to $1/3$ you get $\theta=1$. $\endgroup$
    – Matt F.
    Commented Mar 11, 2020 at 17:10
  • $\begingroup$ Oh, you were simply giving an identity, you weren't implicitly giving a definition of theta in terms of r... Rewriting it is a bit tricky for me, but I can use one of R's built-in optimizers to find the right theta. Redoing it that way, and increasing the Monte Carlo iterations to ensure accuracy, I find that the hypergeometric function seems to explode with NaNs for r > 0.20 (the library's fault?) and that while it closely matches the simulation initially, your formula underestimates increasingly more as n increases: pastebin.com/0fx4WTTR Is that level of error expected for this? $\endgroup$
    – gwern
    Commented Mar 11, 2020 at 20:35
  • 1
    $\begingroup$ The formulas are a good match for the empirical results from an AMH copula, even though that is not quite the same as a binormal distribution. You can see this with library(copula) and : p_max_amh_montecarlo <- function(n,r) { theta <- r_to_theta(r) mean(replicate(10000, { sample <- rCopula(n, amhCopula(theta)) which.max(sample[,1])==which.max(sample[,2]) }))} $\endgroup$
    – Matt F.
    Commented Mar 12, 2020 at 2:03
  • 1
    $\begingroup$ OK. I was worried I might still have an implementation error. But if that approximation error is expected, then that's fine. It doesn't look like anyone else will provide a better formula anytime soon, and this formula does provide insight of the sort I was hoping for, in explicitly showing even in finite-samples that the correlation is equivalent to a constant factor boost and fades out with n, so I will accept this answer. Thanks. $\endgroup$
    – gwern
    Commented Mar 12, 2020 at 2:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.