probability of bivariate max (2)

r_to_theta <- function (r) {
   identity_solver <- function(theta, r) {
       target = 2/pi * asin(r)
       theta_guess = 1 - ((2 * ((1 - theta)^2 * log(1 - theta) + theta))/(3 * theta^2))
       return(abs(target - theta_guess)) }

    theta_approx <- optim(0.5, identity_solver, r=r, method="Brent", lower=0, upper=1)$par
    theta_approx
    }

p_max_bivariate <- function(n,r) {
  ## for the special-case r=0-0.5:
  if (r > 0.5 || r < 0.0) { stop('Formula valid only for 0 <= r <= 0.5') }
    else {
        theta <- r_to_theta(r)
        t <- theta / (theta - 1)
        library(gsl) # for hyperg_2F1

        t * (hyperg_2F1(1, n+1, n+2, t) / (n+1)) -
         2 * n * t^2 * (hyperg_2F1(1, n+2, n+3, t) / ((n+1)*(n+2))) -
         (2*n*t)/((n+1)^2) +
         (1/n)
        }
}


library(MASS)
max1IsMax2 <- function(n,r) { mean(replicate(10000000, {
                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.20088695
round(digits=3, sapply(2:10, function(n) { max1IsMax2(n, r) }))
# [1] 0.564 0.408 0.325 0.274 0.238 0.212 0.191 0.175 0.162
round(digits=3, sapply(2:10, function(n) { p_max_bivariate(n, r) }))
# [1] 0.564 0.403 0.316 0.260 0.221 0.193 0.171 0.153 0.139