Advertisement
Guest User

probability of bivariate max (2)

a guest
Mar 11th, 2020
101
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
Advertisement
Advertisement
Advertisement
RAW Paste Data Copied
Advertisement