Integrating pdf times cdf squared

Question

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.

Answer 1

$$\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.)

Answer 2

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,

1. 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

2. 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.