Qué es un método para calcular precisamente $P(Y \geq X, Y\leq Z)$, dado tres variables al azar independientes $X, Y$ y $Z$

Question

Para tres independiente distribuidas normalmente variables aleatorias continuas X, Y y Z (cada uno con su propia media y desviación estándar), necesito una manera de calcular

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

Sé que puedo hacerlo por el siguiente:

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

Soy capaz de calcular porcentajes $P(Y \geq X)$utiliza la siguiente relación: $P(Y \geq X) = P(Y - X \geq 0)$

Sin embargo, tengo problemas para calcular precisamente $P(Y \leq Z | Y \geq X)$

Answer 1

Uno relativamente fácil enfoque es considerar $X$, $Y$, y $Z$ como tener una articulación distribución normal multivariante.

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

Vamos

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

Añadido: Aquí el código R para calcular la probabilidad.

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

Sólo podría hacer muchos sorteos de la distribución y calcular la tasa que se produzca el evento que te interesa. En 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

Es sólo una estimación así que tal vez hacerlo un par de veces. Rápida y sucia