Integrating pdf times cdf squared I'm trying to solve a distribution problem. As part of it, I need to find $\int\limits_{-\infty}^\infty\phi(z)\Phi^2(z+q)dz$ and $\int\limits_{-\infty}^\infty\phi(z)\Phi(z)\Phi(z+q)dz$
So far, I've worked out that $\int\phi(z)\Phi^2(z)dz=\frac{1}{3}$ and that $\int\limits_{-\infty}^\infty\phi(z)\Phi(z+q)dz = \Phi\left(\frac{q}{\sqrt2}\right)$, but I can't get any further.
Edit: to get the first identity:
if $u = \Phi(z)$ then $du = \phi(z)dz$
$\int\phi(z)\Phi^2(z)dz$ becomes $\int u^2du = \left[\frac{u^3}{3}\right]_{-\infty}^\infty= \frac{\Phi^3(\infty)-\Phi^3(-\infty)}{3}=\frac{1}{3}$
The second comes from here:
https://mathoverflow.net/questions/101469/integration-of-the-product-of-pdf-cdf-of-normal-distribution
They are standard normal distributions.
 A: $$\int_{-\infty}^\infty \phi(z)\Phi^2(z+q)\,\mathrm dz=\int_{-\infty}^\infty \phi(y-q)\Phi^2(y)\,\mathrm dy.$$ Integrals of this type have no known analytical solution in terms of other well-known functions, but see David Robertson's comment re special functions.  
The value of the first integral that you seek is the probability that a $N(q,1)$ random variable is larger than both of two other independent $N(0,1)$ random variables. See this answer of mine for more discussion of the case when $\Phi^2(y)$ is replaced by $\Phi^n(y), n \geq 2.$ Indeed,
the value of the first integral as a probability makes the interpretation
of your derivation of $\int\phi(z)\Phi^2(z)dz=\frac{1}{3}$ straightforward: the probability that one specific random variable
among $3$ iid random variables has the largest value is $\frac 13$
(since all three random variables are equally likely to be the largest.)
A: When we follow the way suggested by Dilip Sarwate, we see that
\begin{eqnarray}\int_{-\infty}^{\infty}\Phi(z+q)^{2}\phi(z)dz=\mathbb{E}\left(P(x\leq z+q,y\leq z+q)\right)=P(x-z\leq q,y-z\leq q)\\=
\mathcal{MVN}\left(x=\{q,q\},\mu=\{0,0\},\Sigma=\begin{bmatrix}\sqrt{2} & 1 &\\ 1& \sqrt{2}\end{bmatrix}\right)\end{eqnarray}
where $\mathcal{MVN}$ is multivariate normal distribution CDF.
A quick R implementation results as follows,

*

*We define a function g for to inside the integral and then numerically integrate

g = function(a,b,x) {pnorm(a+b*x)^2*dnorm(x)}
integrate(g, -Inf, Inf, a=1., b=1.)

[1]0.633702 with absolute error < 4.7e-05



*We use the multivariate normal CDF we found above,

pmnorm(x = c(1.,1.), mean = rep(0.,2), varcov = matrix(c(2,1.,1.,2), 
       ncol=2, byrow=T))

[1] 0.633702

We see that they coincide and it might be possible to generalize this for the power $n$ using $\mathcal{MVN}$ with $n\times n$ dimensional covariance matrix.
