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)$

  • 3
    $\begingroup$ Hint: Conditioned on $Y = y$, the event $\{Y \geq X, Y \leq Z\}$ has conditional probability $$P\{Y \geq X, Y \leq Z \mid Y = y\} = P\{X \leq y, Z \geq y \mid Y = y\} = \Phi((y-\mu_X)/\sigma_x)(1 - \Phi((y-\mu_z)/\sigma_Z)).$$ Multiply by the density of $Y$ and integrate. $\endgroup$ May 24, 2012 at 2:39
  • $\begingroup$ Is this a homework problem? Darn it. I have an answer I'm proud of, but I'm not sure if I ought to post it. $\endgroup$
    – Cyan
    May 24, 2012 at 3:33
  • $\begingroup$ Just to check - are you calculating $P(x \leq y \leq z)$ or do you already know that $x \leq z$ and are calculating the probability that $y$ is inbetween? $\endgroup$
    – jbowman
    May 24, 2012 at 3:56
  • $\begingroup$ Cyan, no it's not homework, it's just a problem that was driving me mad. $\endgroup$
    – Abey
    May 24, 2012 at 4:04
  • $\begingroup$ Hi jbowman. x is not necessarily <= z, so it would be P(x <= y <= z). Thanks. $\endgroup$
    – Abey
    May 24, 2012 at 4:05

2 Answers 2


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)$


$\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.

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)
  • $\begingroup$ Thank you Cyan. Although I'm not experienced on the level you post about, I will study upon it. Much appreciated. :) $\endgroup$
    – Abey
    May 24, 2012 at 4:20
  • $\begingroup$ What pieces are you missing? Maybe I can point you to some resources. $\endgroup$
    – Cyan
    May 24, 2012 at 4:26
  • $\begingroup$ Honestly, I am not very familiar with the standardized ways of expressing these concepts. Thus, the meaning of those first three parts is unclear to me. I understand the last line; however, the means of the three distributions are not necessarily the same. $\endgroup$
    – Abey
    May 24, 2012 at 4:34
  • $\begingroup$ The bivariate normal distribution might be a good starting point -- the linked page has nice animations. Then you'll need to know about vectors (especially in Cartesian space) and matrices. (cont'd) $\endgroup$
    – Cyan
    May 24, 2012 at 5:11
  • $\begingroup$ After that you'll be ready to tackle random vectors and the multivariate normal distribution. $\endgroup$
    – Cyan
    May 24, 2012 at 5:11

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

 sum(x<y & y >z )/N

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

  • 1
    $\begingroup$ The final line of code should be "sum(x<y & y<z )/N". $\endgroup$
    – Cyan
    May 31, 2012 at 13:04
  • $\begingroup$ This is a great way to check an answer. But did you notice the word "precisely" in the question? $\endgroup$
    – whuber
    May 31, 2012 at 14:13
  • $\begingroup$ Thanks Seth, and ditto on good idea to check the answer. $\endgroup$
    – Abey
    May 31, 2012 at 15:11
  • $\begingroup$ I noticed the word 'precisely', I just ignored it. Just kidding. Both cyan's and this method can be computed with any level of precision. $\endgroup$
    – Seth
    May 31, 2012 at 15:15
  • $\begingroup$ Up to a point Seth: the precision is inherently limited by your computing capabilities and lifetime of the universe. Remember, the SD of a simulated estimate scales only as $N^{-1/2}$. That will make it prohibitively difficult to attain more than about six significant figures; your example only gets about four sig figs. (Every once in a while even six sig figs is not good enough as a check, so this is not entirely nit-picking.) $\endgroup$
    – whuber
    May 31, 2012 at 18:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.