Integral of product of Gaussian pdf and cdf $\phi(x) = \frac{1}{\sqrt{2\pi}} e^{-\frac{x^2}{2}}$ is the pdf of a standard normal distribution.
$\Phi(x) = \int_{- \infty}^x \phi(t) dt$ is the cdf of a standard normal distribution.
How does one calculate the following: 
$$
T = \int\limits_0^{\infty} \Phi^2(bx) \phi(x) dx 
$$
In fact, I can find the result in wiki,
https://en.wikipedia.org/wiki/List_of_integrals_of_Gaussian_functions#CITEREFOwen1980
But I would like to know how to do it.
 A: Maybe this proof is too large, but you can deduce the result you want, sorry for that.
First of all:


*

*$\phi(x)= \frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}}$

*$\Phi(x)=\int_{-\infty}^{x}\phi(t)dt=\frac{1}{2}(1+erf(\frac{x}{\sqrt{2}}))$


To prepare the proof, start with this facts:


*

*$\int_{0}^{\infty}xe^{-x^2a}dx=\frac{1}{2a}$ (Proof by direct integral)

*$erf'(x)=\frac{2}{\sqrt{\pi}}e^{-x^2}$ (Proof by definition of $erf$ function)

*$$\int_0^{\infty}xe^{-\alpha^2x^2}erf(\beta x)dx=\frac{\beta}{2\alpha^2 \sqrt{\alpha^2+\beta^2}} $$
Hint for 3.
Integrating by parts, with 
\begin{equation}u=erf(\beta x) \rightarrow du=\frac{2\beta}{\sqrt{\pi}}e^{-(x\beta)^2}dx\end{equation}
$$dv=xe^{-\alpha^2x^2} \rightarrow v=-\frac{e^{-x^2\alpha^2}}{2\alpha^2}$$
the result follows really easy.


Now, you have all the ingredients to calculate this monster:
$$T(b)=\int_{0}^{\infty}\Phi^2(bx)\phi(x)dx$$
Derivating $T(b)$ respect $b$, changing $\Phi(bx)$ by the definition in terms of $erf$ function, and using the fact 2, will lead you to:
$$T'(b)=...=\frac{1}{2\pi}\Big(\int_0^{\infty}{xe^{-x^2(b^2+1)/2}}dx + \int_0^{\infty}{x \cdot erf\big(\frac{bx}{\sqrt{2}}\big)e^{-x^2(b^2+1)/2}}dx \Big) =\frac{1}{2\pi}(A+B)$$
Compute A using the fact 1 to get:
$$A=\frac{1}{b^2+1}$$
and compute B using the fact 3 with $\alpha=\sqrt{\frac{b^2+1}{2}}$ and $\beta= \frac{b}{\sqrt{2}}$ to get:
$$B=\frac{b}{(b^2+1)\sqrt{2b^2+1}}.$$
It remains to integrate $T'(b)=\frac{A+B}{2\pi}$ respect $b$:
\begin{align}T(b)=&\frac{1}{2\pi}\Big(\int_{0}^{b}\frac{1}{t^2+1}dt+\int_{0}^{b}\frac{t}{(t^2+1)\sqrt{2t^2+1}}dt\Big)=\\
=&\frac{1}{2\pi}\Big(arctan(b)+arctan\big(\sqrt{2b^2+1}\big)\Big)+C. \hspace{0.5cm}\square
\end{align}
A: By a rather circuitous route I get $$T = \dfrac{\arctan\left(\sqrt{2b^2+1}\right) + \arctan(b)}{2\pi}$$
