# probability of bivariate max (2)

a guest
Mar 11th, 2020
126
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. r_to_theta <- function (r) {
2.    identity_solver <- function(theta, r) {
3.        target = 2/pi * asin(r)
4.        theta_guess = 1 - ((2 * ((1 - theta)^2 * log(1 - theta) + theta))/(3 * theta^2))
5.        return(abs(target - theta_guess)) }
6.
7.     theta_approx <- optim(0.5, identity_solver, r=r, method="Brent", lower=0, upper=1)\$par
8.     theta_approx
9.     }
10.
11. p_max_bivariate <- function(n,r) {
12.   ## for the special-case r=0-0.5:
13.   if (r > 0.5 || r < 0.0) { stop('Formula valid only for 0 <= r <= 0.5') }
14.     else {
15.         theta <- r_to_theta(r)
16.         t <- theta / (theta - 1)
17.         library(gsl) # for hyperg_2F1
18.
19.         t * (hyperg_2F1(1, n+1, n+2, t) / (n+1)) -
20.          2 * n * t^2 * (hyperg_2F1(1, n+2, n+3, t) / ((n+1)*(n+2))) -
21.          (2*n*t)/((n+1)^2) +
22.          (1/n)
23.         }
24. }
25.
26.
27. library(MASS)
28. max1IsMax2 <- function(n,r) { mean(replicate(10000000, {
29.                 sample <- mvrnorm(n, mu=c(0,0), Sigma=matrix(c(1, r, r, 1), nrow=2))
30.                 which.max(sample[,1])==which.max(sample[,2])
31.                 })) }
32.
33. r=0.20088695
34. round(digits=3, sapply(2:10, function(n) { max1IsMax2(n, r) }))
35. # [1] 0.564 0.408 0.325 0.274 0.238 0.212 0.191 0.175 0.162
36. round(digits=3, sapply(2:10, function(n) { p_max_bivariate(n, r) }))
37. # [1] 0.564 0.403 0.316 0.260 0.221 0.193 0.171 0.153 0.139