Stack Exchange answer PR(X<Y|min(X,Y))

script.R

library(MASS)
library(mvtnorm)

# Set parameters for generating simulated data
set.seed(1)
mu = c(.1,1.5)
Sigma = matrix(c(1.5,.5,.5,.7),2,2)

# Set minumum
m = 1

# Set half of conditional set size
eps = .01

storage = matrix(, nrow = 0, ncol = 2)

while(dim(storage)[1] < 1000){
  
  X = mvrnorm(100000,mu,Sigma)
  # Calculate pairwise minimums
  mins = pmin(X[,1],X[,2])
  
  # Create subset for conditional statement
  dex = mins < m + eps & mins > m - eps # min(x,y) = m 
  
  storage = rbind(storage,X[dex,])
  
}


# Calculate conditional distribution parameters
muxy = mu[1]+Sigma[1,2]*(m-mu[2])/Sigma[2,2]
muyx = mu[2]+Sigma[1,2]*(m-mu[1])/Sigma[1,1]

s2xy = (1 - Sigma[1,2]^2/(Sigma[2,2]*Sigma[1,1]))*Sigma[1,1]
s2yx = (1 - Sigma[1,2]^2/(Sigma[1,1]*Sigma[2,2]))*Sigma[2,2]


# Check it
out1 = mean(storage[,1]<storage[,2]) 

out2 = dnorm((m-mu[1])/sqrt(Sigma[1,1]))/sqrt(Sigma[1,1]) * (1-pnorm((m-muyx)/sqrt(s2yx))) / 
      (dnorm((m-mu[1])/sqrt(Sigma[1,1]))/sqrt(Sigma[1,1]) * (1-pnorm((m-muyx)/sqrt(s2yx))) + 
       dnorm((m-mu[2])/sqrt(Sigma[2,2]))/sqrt(Sigma[2,2]) * (1-pnorm((m-muxy)/sqrt(s2xy))))

c(out1,out2)