Join Probability of three exponential random variables? I'm trying to derive the probability $Pr[Y<a+bZ,Y<X,Z<X]$, where $X$, $Y$, and $Z$ are independent exponential distributed random variables with parameters $1/\lambda_x$, $1/\lambda_y$, and $1/\lambda_z$, respectively.
What I did is as follows:
$\Pr \left[ {Y < a + bZ,Y < X,Z < X} \right] = \int\limits_0^\infty  {\int\limits_0^x {\int\limits_0^{a + bz} {{f_{X,Y,Z}}\left( {x,y,z} \right)dydzdx} } } \\
 = \int\limits_0^\infty  {\int\limits_0^x {\int\limits_0^{a + bz} {{f_X}\left( x \right){f_Y}\left( y \right){f_Z}\left( z \right)dydzdx} } } \\
 = 1 - \frac{{{e^{\frac{{ - a}}{{{\lambda _y}}}}}}}{{{\lambda _z}}}{\left( {\frac{b}{{{\lambda _y}}} + \frac{1}{{{\lambda _z}}}} \right)^{ - 1}} - \frac{1}{{{\lambda _x}}}{\left( {\frac{1}{{{\lambda _x}}} + \frac{1}{{{\lambda _z}}}} \right)^{ - 1}} + \frac{{{e^{\frac{{ - a}}{{{\lambda _y}}}}}}}{{{\lambda _x}{\lambda _z}}}{\left( {\frac{b}{{{\lambda _y}}} + \frac{1}{{{\lambda _z}}}} \right)^{ - 1}}\left( {\frac{1}{{{\lambda _x}}} + \frac{b}{{{\lambda _y}}} + \frac{1}{{{\lambda _z}}}} \right)$
I did simulation to verify the analytic result. What I got is
.
Do you have any idea where the gap between simulated and analytic curves come from? Thank you very much.
 A: Your limits of integration don't include the restriction $Y\lt X$. To fix it, the integral needs to be split into three parts and the results added:
$$\int_{x=0}^a \int_{z=0}^{x} \int_{y=0}^{x} {f_X(x) f_Y(y) f_Z(z)\;dy\;dz\;dx} + \int_{x=a}^\infty \int_{z=0}^{\frac{x-a}{b}} \int_{y=0}^{a+bz} {f_X(x) f_Y(y) f_Z(z)\;dy\;dz\;dx} + \int_{x=a}^\infty \int_{z=\frac{x-a}{b}}^{x} \int_{y=0}^{x} {f_X(x) f_Y(y) f_Z(z)\;dy\;dz\;dx}.$$
The upper limit $\frac{x-a}{b}$ comes from finding that value of $z$ such that $y$ can range between $0$ and $x$ without exceeding $a+bz$. It is found by solving $x=a+bz$.
For $x\lt a$ we have $\frac{x-a}{b}\lt 0$ so this range for $x$ needs separating.
Note that these integral limits rely on $a\geq 0$ and $b\geq 1$. If this is not the case, the integrals need reworking.
We now calculate the three definite integrals in turn. For easier notation, we let:
$$\alpha=1/\lambda_x \\ \beta=1/\lambda_y \\ \gamma=1/\lambda_z.$$
$$$$
\begin{eqnarray*}
I_1 &=& \int_{x=0}^a \int_{z=0}^{x} \int_{y=0}^{x} {\alpha\beta\gamma\; e^{-\alpha x-\beta y-\gamma z}\;dy\;dz\;dx} \\
&=& \int_{x=0}^a \int_{z=0}^{x} \left[ -\alpha\gamma\; e^{-\alpha x-\beta y-\gamma z}\right]_{y=0}^{x} \;dz\;dx \\
&=& \int_{x=0}^a \int_{z=0}^{x} \left( -\alpha\gamma\; e^{-\alpha x-\beta x-\gamma z} + \alpha\gamma\; e^{-\alpha x-\gamma z}\right) \;dz\;dx \\
&=& \int_{x=0}^a \left[ \alpha\; e^{-\alpha x-\beta x-\gamma z} - \alpha\; e^{-\alpha x-\gamma z}\right]_{z=0}^{x} \;dx \\
&=& \int_{x=0}^a \left( \alpha\; e^{-\alpha x-\beta x-\gamma x} - \alpha\; e^{-\alpha x-\gamma x} - \alpha\; e^{-\alpha x-\beta x} + \alpha\; e^{-\alpha x}\right) \;dx \\
&=& \left[ \frac{-\alpha}{\alpha+\beta+\gamma }\; e^{-\alpha x-\beta x-\gamma x} + \frac{\alpha}{\alpha +\gamma}\; e^{-\alpha x-\gamma x} + \frac{\alpha}{\alpha +\beta}\; e^{-\alpha x-\beta x} - e^{-\alpha x}\right]_{x=0}^a \\
&& \\
&=& \frac{-\alpha}{\alpha+\beta+\gamma }\; e^{-a(\alpha+\beta+\gamma)} + \frac{\alpha}{\alpha +\gamma}\; e^{-a(\alpha+\gamma)} + \frac{\alpha}{\alpha +\beta}\; e^{-a(\alpha +\beta)} - e^{-a\alpha} \\
&&\qquad + \frac{\alpha}{\alpha+\beta+\gamma } - \frac{\alpha}{\alpha +\gamma} - \frac{\alpha}{\alpha +\beta} + 1. \\
\end{eqnarray*}
$$\\$$
\begin{eqnarray*}
I_2 &=& \int_{x=a}^{\infty} \int_{z=0}^{(x-a)/b} \int_{y=0}^{a+bz} {\alpha\beta\gamma\; e^{-\alpha x-\beta y-\gamma z}\;dy\;dz\;dx} \\
&=& \int_{x=a}^{\infty} \int_{z=0}^{(x-a)/b} \left[ -\alpha\gamma\; e^{-\alpha x-\beta y-\gamma z}\right]_{y=0}^{a+bz} \;dz\;dx \\
&=& \int_{x=a}^{\infty} \int_{z=0}^{(x-a)/b} \left( -\alpha\gamma\; e^{-\alpha x-\beta (a+bz)-\gamma z} + \alpha\gamma\; e^{-\alpha x-\gamma z}\right) \;dz\;dx \\
&=& \int_{x=a}^{\infty} \left[ \frac{\alpha\gamma}{b\beta+\gamma}\; e^{-\alpha x-\beta (a+bz)-\gamma z} - \alpha\; e^{-\alpha x-\gamma z}\right]_{z=0}^{(x-a)/b} \;dx \\
&=& \int_{x=a}^{\infty} \left( \frac{\alpha\gamma}{b\beta+\gamma}\; e^{-\alpha x-\beta x-\gamma (x-a)/b} - \alpha\; e^{-\alpha x-\gamma (x-a)/b} - \frac{\alpha\gamma}{b\beta+\gamma}\; e^{-\alpha x-a\beta} + \alpha\; e^{-\alpha x}\right) \;dx \\
&=& \left[ \frac{-\alpha\gamma}{(b\beta+\gamma)(\alpha+\beta+\gamma/b) }\; e^{-\alpha x-\beta x-\gamma (x-a)/b} + \frac{\alpha}{\alpha +\gamma/b}\; e^{-\alpha x-\gamma (x-a)/b} + \frac{\gamma}{b\beta + \gamma}\; e^{-\alpha x-a\beta} - e^{-\alpha x}\right]_{x=a}^\infty \\
&& \\
&=& \frac{\alpha\gamma}{(b\beta+\gamma)(\alpha+\beta+\gamma/b) }\; e^{-a(\alpha +\beta)} - \frac{\alpha}{\alpha +\gamma/b}\; e^{-a\alpha} - \frac{\gamma}{b\beta + \gamma}\; e^{-a(\alpha +\beta)} + e^{-a\alpha}. \\
\end{eqnarray*}
$$\\$$
\begin{eqnarray*}
I_3 &=& \int_{x=a}^{\infty} \int_{z=(x-a)/b}^{x} \int_{y=0}^{x} {\alpha\beta\gamma\; e^{-\alpha x-\beta y-\gamma z}\;dy\;dz\;dx} \\
&=& \int_{x=a}^{\infty} \int_{z=(x-a)/b}^{x} \left[ -\alpha\gamma\; e^{-\alpha x-\beta y-\gamma z}\right]_{y=0}^{x} \;dz\;dx \\
&=& \int_{x=a}^{\infty} \int_{z=(x-a)/b}^{x} \left( -\alpha\gamma\; e^{-\alpha x-\beta x-\gamma z} + \alpha\gamma\; e^{-\alpha x-\gamma z}\right) \;dz\;dx \\
&=& \int_{x=a}^{\infty} \left[ \alpha\; e^{-\alpha x-\beta x-\gamma z} - \alpha\; e^{-\alpha x-\gamma z}\right]_{z=(x-a)/b}^{x} \;dx \\
&=& \int_{x=a}^{\infty} \left( \alpha\; e^{-\alpha x-\beta x-\gamma x} - \alpha\; e^{-\alpha x-\gamma x} - \alpha\; e^{-\alpha x-\beta x-\gamma(x-a)/b} + \alpha\; e^{-\alpha x-\gamma(x-a)/b}\right) \;dx \\
&=& \left[ \frac{-\alpha}{\alpha+\beta+\gamma}\; e^{-\alpha x-\beta x-\gamma x} + \frac{\alpha}{\alpha +\gamma}\; e^{-\alpha x-\gamma x} + \frac{\alpha}{\alpha+\beta+\gamma/b}\; e^{-\alpha x-\beta x-\gamma(x-a)/b} - \frac{\alpha}{\alpha+\gamma/b}\; e^{-\alpha x - \gamma(x-a)/b}\right]_{x=a}^\infty \\
&& \\
&=& \frac{\alpha}{\alpha+\beta+\gamma}\; e^{-a(\alpha +\beta+\gamma)} - \frac{\alpha}{\alpha +\gamma}\; e^{-a(\alpha+\gamma)} - \frac{\alpha}{\alpha+\beta + \gamma/b}\; e^{-a(\alpha +\beta)} + \frac{\alpha}{\alpha+\gamma/b}\; e^{-a\alpha}. \\
\end{eqnarray*}
