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?
 A: 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.
