What is a method to calculate precisely $P(Y \geq X, Y\leq Z)$, given three independent random variables $X, Y$, and $Z$

Question

For three independent normally distributed continuous random variables X, Y, and Z (each with its own mean and standard deviation), I need a way to calculate

$P(Y \geq X, Y \leq Z)$

I know that I can do this by the following:

$P(Y \geq X, Y \leq Z) = P(Y \geq X) \cdot P(Y \leq Z | Y \geq X)$

I am able to calculate $P(Y \geq X)$ using the following relation: $P(Y \geq X) = P(Y - X \geq 0)$

However, I'm having trouble calculating precisely $P(Y \leq Z | Y \geq X)$

Answer 1

One relatively easy approach is to consider $X$, $Y$, and $Z$ as having a joint multivariate normal distribution.

$\left[\begin{array}{c}X\\Y\\Z\end{array}\right]\sim\mathrm{MVN}\left(\left[\begin{array}{c}\mu_{X}\\\mu_{Y}\\\mu_{Z}\end{array}\right],\left[\begin{array}{ccc}\sigma_{X}^{2} & 0 & 0\\0 & \sigma_{Y}^{2} & 0\\0 & 0 & \sigma_{Z}^{2}\end{array}\right]\right)$

Let

$\left[\begin{array}{c} U\\ V\end{array}\right]=\left[\begin{array}{c} X-Y\\ Z-Y\end{array}\right]=\left[\begin{array}{ccc} 1 & -1 & 0\\ 0 & -1 & 1\end{array}\right]\left[\begin{array}{c} X\\ Y\\ Z\end{array}\right]$

Then by standard results on affine transformations of multivariate normal distributions,

$\left[\begin{array}{c} U\\ V\end{array}\right]\sim\mathrm{MVN}\left(\left[\begin{array}{c} \mu_{X}-\mu_{Y}\\ \mu_{Z}-\mu_{Y}\end{array}\right],\left[\begin{array}{cc} \sigma_{X}^{2}+\sigma_{Y}^{2} & \sigma_{Y}^{2}\\ \sigma_{Y}^{2} & \sigma_{Z}^{2}+\sigma_{Y}^{2}\end{array}\right]\right)$

And since $P(Y \geq X, Y \leq Z) = P(U \leq 0, V \geq 0)$, you want the probability mass of this bivariate distribution in the second quadrant. This is not analytically solvable in general, but is easy to compute. If $\mu_X = \mu_Y = \mu_Z$, then there is an analytical expression (from equation 73 here):

$P(U \leq 0, V \geq 0) = \frac{1}{2} \cos^{-1}\left(\frac{\sigma^2_{Y}}{\sqrt{(\sigma^2_{X} + \sigma^2_{Y}) (\sigma^2_{Z} + \sigma^2_{Y})}}\right)$.

Added: Here's R code to compute the probability.

install.packages("mvtnorm")
library(mvtnorm)
mu_x <- -1.4
mu_y <- 2
mu_z <- 1.7
mu_vec <- c(mu_x- mu_y, mu_z - mu_y) 
var_x <- 9
var_y <- 9
var_z <- 16
Sigma <- var_y + matrix(c(var_x, 0, 0 , var_z), nrow = 2)
pmvnorm(lower = c(-Inf, 0), upper = c(0, Inf), mean = mu_vec, sigma = Sigma)

Answer 2

I might just make many draws from the distribution and calculate the rate that the event you are interested in occurs. In R:

N=10^7
 x=rnorm(N,mu_x,sig_x)
 y=rnorm(N,mu_y,sig_y)
 z=rnorm(N,mu_z,sig_z)
 sum(x<y & y >z )/N

It is just an estimation so maybe do it a couple times. Quick and dirty