Given two means and deviations, how can I compute the probability that x < y? I've experimentally achieved my goal by running many random trials, generating two points according to scaled and translated gaussian distributions and counting how many times x is less than y.
However this is becoming a problem from a performance point of view for my application.
Is there a straightforward way to compute P(X < Y) where X and Y are variables with a given mean and deviation?
 A: If your assumptions regarding the normality of $X$ and $Y$ are correct, then you have two normally distributed random variables $X\sim\mathcal{N}(\mu_X, \sigma_X^2)$ and $Y\sim\mathcal{N}(\mu_Y, \sigma_Y^2)$, and you are searching for $\mathrm{Pr}(X < Y)$, aka $\mathrm{Pr}(X-Y < 0)$.
Assuming $X$ and $Y$ are independent, this is easy to analyze, because in that case $X-Y\sim\mathcal{N}(\mu_X-\mu_Y, \sigma^2_X+\sigma^2_Y)$. Therefore your desired probability can be obtained from the standard normal cdf: $\mathrm{Pr}(X-Y < 0) = \Phi(\frac{\mu_Y-\mu_X}{\sqrt{\sigma^2_X+\sigma^2_Y}})$.
To link this to your approach of simulating many samples from the random variables, consider the case where $X\sim\mathcal{N}(1, 3^2)$ and $Y\sim\mathcal{N}(2, 4^2)$. Then we can calculate that $\mathrm{Pr}(X < Y) = \Phi(\frac{2-1}{\sqrt{3^2+4^2}})\approx 0.57926$. We can obtain a similar result via simulation in R:
set.seed(144)
x.samples <- rnorm(1e6, 1, 3)
y.samples <- rnorm(1e6, 2, 4)
mean(x.samples < y.samples)
# [1] 0.579655

If $X$ and $Y$ are not independent, then we need more information to calculate $\mathrm{Pr}(X < Y)$. For instance, if they are jointly normally distributed with correlation coefficient $\rho$, then $X-Y\sim\mathcal{N}(\mu_X-\mu_Y, \sigma_X^2+\sigma_Y^2+2\rho\sigma_X\sigma_Y)$, meaning $\mathrm{Pr}(X < Y) = \Phi(\frac{\mu_Y-\mu_X}{\sqrt{\sigma_X^2+\sigma_Y^2+2\rho\sigma_X\sigma_Y}})$. If they are not jointly normally distributed, then $X-Y$ may have some other distribution, and that will impact the calculation of $\mathrm{Pr}(X < Y)$.
