probability of bivariate max (error?)

r_to_theta <- function (r) { 2/pi * asin(r) }
round(digits=3, sapply(seq(0,0.5, by=0.05), function(r) { r_to_theta(r) }))
# [1] 0.000 0.032 0.064 0.096 0.128 0.161 0.194 0.228 0.262 0.297 0.333

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)
        }
}
r=0.5; sapply(2:10, function(n) { p_max_bivariate(n, r) })
# [1] 0.540620144 0.378139568 0.292465038 0.238966583 0.202214466 0.175349622 0.154830414 0.138634576 0.125520185