What is a method to calculate precisely $P(Y \geq X, Y\leq Z)$, given three independent random variables $X, Y$, and $Z$ 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)$
 A: 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)

A: 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
