Is it possible to calculate a probability by plugging a random variable into the cumulative distribution function? I was calculating $p = {\rm Pr}(X < Y)$ where both $X$ and $Y$ are independent exponentially distributed with rate parameters $\lambda_X$ and $\lambda_Y$ respectively. I had the idea that, since the $p$ is just a constant, $p = \mathbb{E}[p]$, and since ${\rm Pr}(k < Y) = 1 - F_Y(k) = e^{-k\lambda_Y}$, then maybe \begin{align}p = \mathbb{E}[1-F(X)] = \mathbb{E}[e^{-X\lambda_Y}] = \int_0^\infty e^{-x\lambda_Y}\lambda_Xe^{-x\lambda_X} dx  = \frac{\lambda_X}{\lambda_X+\lambda_Y}\end{align} which is the correct answer as $$\int_0^\infty (\lambda_X+\lambda_Y)e^{-(\lambda_X+\lambda_Y)x}dx = \frac{\lambda_X}{\lambda_X+\lambda_Y}.$$ However, when I tried using the same logic to calculate $${\rm Pr}(Z < X + Y) = 1 - \frac{\lambda_X\lambda_Y}{(\lambda_X + \lambda_Z)(\lambda_Y + \lambda_Z)} $$
using the fact that $X$ and $Y$ are independent to separate the expected values however, some preliminary simulation results suggest this is not true. Is there something I am missing?
 A: If I were you, I would double check your simulation. I think the result you are claiming is true. You can formalize the argument with conditional expectations:
$$\mathsf{P}\left(Z<X+Y\right)=\mathsf{E}\left[\mathbf{1}_{Z<X+Y}\right]=\mathsf{E}\left[\mathsf{E} \left[\mathbf{1}_{Z<X+Y} \mid(X,Y)\right]\right]$$
using the law of total expectation. Now as you noticed,
$$\mathsf{E} \left[\mathbf{1}_{Z<X+Y} \mid(X,Y)\right]=\mathsf{E} \left[1-e^{-\lambda_z(X+Y)} \mid(X,Y)\right]=1-\mathsf{E} \left[e^{-\lambda_z X} \mid X\right]\mathsf{E} \left[e^{-\lambda_z Y} \mid Y\right].$$
Note that the last equality is always true, even if $X$ and $Y$ were not independent. It follows that:
$$\mathsf{P}\left(Z<X+Y\right)=1-\mathsf{E} \left[\mathsf{E} \left[e^{-\lambda_z X} \mid X\right]\mathsf{E} \left[e^{-\lambda_z Y} \mid Y\right]\right]=\mathsf{E} \left[\mathsf{E} \left[e^{-\lambda_z X} \mid X\right]\right]\mathsf{E}\left[\mathsf{E} \left[e^{-\lambda_z Y} \mid Y\right]\right]$$
where this time we used that $X$ and $Y$ are independent. Therefore
$$\mathsf{P}\left(Z<X+Y\right)=1-\frac{\lambda_X\lambda_Y}{(\lambda_X+\lambda_Z)(\lambda_Y+\lambda_Z)}.$$

Here is a "simulation" in Python:
from numpy.random import exponential

N = 50000
lx, ly, lz = 4.434, 2.343, 3.789
s = 0
for i in range(N):
    x, y, z = exponential(1 / lx), exponential(1 / ly), exponential(1 / lz)
    if x < y + z:
        s += 1
print(s / N, 1 - ly * lz / (lx + ly) / (lx + lz))

