What is the probability that a Gaussian random variable is the largest of three Gaussian random variables? Consider three independent,  normally distributed RVs: $YA \sim N(a,\sigma ^{2}),$ $%
YB\sim N(b,\sigma ^{2})$ and $YC\sim N(c,\sigma ^{2})$.
What is the probability that $YA$ is the maximum?: $$\Pr (YA>YB \;\;\mbox{and} \;\;YA>YC)=IA=\int_{-\infty }^{\infty }\Phi (\frac{Y-b}{%
\sigma })\Phi (\frac{Y-c}{\sigma })\phi (\frac{Y-a}{\sigma })dY$$
My first thought was to differentiate $IA$ with respect to $b$  and $c$, complete the square and exploit that $$\int_{-\infty }^{\infty }\frac{\sqrt{3}}{\sigma}\phi \left(\frac{Y-\frac{a+b+c}{2}}{%
\sigma /\sqrt{3} }\right)dY=1$$ to remove the integral over $Y$. 
Therefore, $$\frac{d^2IA}{dbdc}\propto \exp \left(-\frac{a^2-2 a y+b^2-2 b y+c^2-2 cy+3 y^2}{2 \sigma ^2}\right)$$.
Unfortunately, I cannot figure how to integrate with respect to $b$ and $c$ to get back to $IA$. 
Any suggestions would be much appreciated. 
 A: Since Google brought me here, I may as well share where I got solving this exact problem.
Note that:
$\Pr \left( YA>YB\text{ and }YA>YC\right) =\Pr \left( YB-YA<0\text{ and } YC-YA<0\right) $
and the distribution of $YB-YA\mid YA$ is simply $\phi _{\mu _{B}-YA,\sigma_{B}^{2}}$
The probability density of $YA=x$ and $YB-x=y$ and $YC-x=z$ is
$\Pr \left( YA=x\cap YB-x=y\cap YC-x=z\right) =\phi _{\mu _{B}-x,\sigma_{B}^{2}}\left( y\right) \phi _{\mu _{C}-x,\sigma _{C}^{2}}\left( z\right)\phi _{\mu _{A},\sigma _{A}^{2}}\left( x\right) $
Hence
$\Pr \left( YA>YB\text{ and }YA>YC\right) =\int_{-\infty }^{\infty }\left(\int_{-\infty }^{0}\phi _{\mu _{B}-x,\sigma _{B}^{2}}\left( y\right)dy\int_{-\infty }^{0}\phi _{\mu _{C}-x,\sigma _{C}^{2}}\left( z\right)dz\right) \phi _{\mu _{A},\sigma _{A}^{2}}\left( x\right) dx$
This expression can be simplified by convolution to get rid of the outer integral. But, like the cumulative normal distribution, there is no analytical solution. 
If anyone wants the final answer, you can email me. But, also, if anyone knows a good algorithm for estimating this expression numerically, please email me.
A: Since
$\begin{pmatrix}YA-YB\\YA-YC\end{pmatrix}\stackrel{\text{dist}}{=}
N\left(\begin{pmatrix}a-b\\a-c\end{pmatrix},\begin{pmatrix}2\sigma^2&\sigma^2\\ \sigma^2&2\sigma^2\end{pmatrix}\right)$, then
\begin{align} \mathsf{P}((YA>YB)\cap(YA>YC))&=\mathsf{P}(YA-YB>0,YA-YC>0)\\
&=L\Bigl(\frac{b-a}{\sqrt2\sigma},\frac{c-a}{\sqrt2\sigma},\frac12\Bigr),
\end{align}
where
$$L(h,k,\rho)=\frac{1}{2\pi\sqrt{1-\rho^2}}\int_h^\infty dx\int_k^\infty\exp\Bigl[
-\frac{x^2-2\rho xy+y^2}{2(1-\rho^2)}\Bigr]\,dy.
$$
In the following book: S. Kotz et al., Continuous Multivariate Distributions,  Vol. 1, 2nd Ed(2000), The Ch.46, Sec.4, p.264--. It gave a complete discussion about above bivariate normal integral--table and approximation. 
