Compute probability of union set of four continuous random variables in R

Question

I have four independent variables X1, X2, X3, X4 with a standard normal distribution. I want to calculate the probability of P(X1 > A | X2 > A | X3 > A | X4 > A). I wrote a function in R that correctly calculates the probability for any value of A that's greater 0, but for any value smaller 0 the results are to small:

prob_union_greater <- function(x){ 
   (1 - pnorm(x))*4 - 6*((1 - pnorm(x))^2) + 3*((1 - pnorm(x))^3) - ((1 - pnorm(x))^4)  
} 

I tried to write a similar function for the case of P(X1 < A | X2 < A | X3 < A | X4 < A) and here I'm faced with the opposite problem: for negative values of A it works, for positive values it doesn't.

prob_union_smaller <- function(x){
   pnorm(x)*4 - 6*(pnorm(x)^2) + 3*(pnorm(x)^3) - (pnorm(x)^4)  
}

What am I missing here?

Answer 1

You could try the code below

prob_union_greater <- function(x) {
  p <- pnorm(x, lower.tail = FALSE)
  4 * p - 6 * p^2 + 4 * p^3 - p^4
}

prob_union_smaller <- function(x) {
  p <- pnorm(x)
  4 * p - 6 * p^2 + 4 * p^3 - p^4
}

and you will get

> prob_union_greater(1)
[1] 0.4989328

> prob_union_smaller(1)
[1] 0.9993664