# Probability of selecting maximum in bivariate correlated order statistics?

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):

library(MASS)
max1IsMax2 <- function(n,r) { mean(replicate(100000, {
sample <- mvrnorm(n, mu=c(0,0), Sigma=matrix(c(1, r, r, 1), nrow=2))
which.max(sample[,1])==which.max(sample[,2])
})) }

r=0.5
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

library(lmomco)
exactMax <- function(n, mean=0, sd=1) { expect.max.ostat(n, para=vec2par(c(mean, sd), type="nor"), cdf=cdfnor, pdf=pdfnor) }
exactMax(10)
# [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?

• 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$. Sep 15, 2019 at 8:01
• 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. :) Sep 15, 2019 at 13:59
• 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. Dec 24, 2019 at 15:45
• 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. Dec 24, 2019 at 15:49

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}$$.

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.

• 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? Mar 11, 2020 at 16:47
• 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$. Mar 11, 2020 at 17:10
• 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? Mar 11, 2020 at 20:35
• 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]) }))} Mar 12, 2020 at 2:03
• 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. Mar 12, 2020 at 2:09