How to compute normal integrals $\int_{-\infty}^\infty\Phi(x)N(x\mid\mu,\sigma^2)\,dx$ and $\int_{-\infty}^\infty\Phi(x)N(x\mid\mu,\sigma^2)x\,dx$ How to compute the following formula?
$$
\int_{-\infty}^{+\infty} \Phi(x) N(x\mid\mu,\sigma^2) \, dx
$$
$$
\int_{-\infty}^{+\infty} \Phi(x) N(x\mid\mu,\sigma^2) x\,dx
$$
where $\Phi(x)=\int_{-\infty}^{x}\frac{1}{\sqrt{2\pi}} \exp{(-t^2/2)} \,  dt$, namely, the cumulative distribution function of normal distribution $N(0,1)$. $N(x\mid\mu,\sigma^2) $ means the probability density function of Gaussian distribution with mean $\mu$ and variance $\sigma^2$
 A: The first integral is $$I_{\mu,\sigma^2}=\int_\mathbb R\Phi(x)\varphi_{\mu,\sigma^2}(x)dx=E(\Phi(\mu+\sigma X))$$ that is, $$I_{\mu,\sigma^2}=P(Y\leqslant \mu+\sigma X)=P(Z\leqslant\nu)=\Phi(\nu)$$ where $(X,Y)$ is i.i.d. standard normal, $$\nu=\frac{\mu}{\sqrt{\sigma^2+1}}$$ and $Z$ is also standard normal since $$Z=\frac{Y-\sigma X}{\sqrt{\sigma^2+1}}$$ Thus, 

$$I_{\mu,\sigma^2}=\int_\mathbb R\Phi(x)\varphi_{\mu,\sigma^2}(x)dx=\Phi\left(\frac{\mu}{\sqrt{\sigma^2+1}}\right)$$

The second integral is $$J_{\mu,\sigma^2}=\int_\mathbb R\Phi(x)\varphi_{\mu,\sigma^2}(x)xdx=\int_\mathbb R\Phi(\mu+\sigma x)\varphi(x)xdx$$ Since $x\varphi(x)=-\varphi'(x)$, an integration by parts yields $$J_{\mu,\sigma^2}=\int_\mathbb R\sigma\varphi(\mu+\sigma x)\varphi(x)dx$$ Now, there exists $(\tau,\kappa,\lambda)$ such that $$x^2+(\mu+\sigma x)^2=(\tau x+\kappa)^2+\lambda$$ hence $$J_{\mu,\sigma^2}=\frac\sigma{\sqrt{2\pi}}e^{-\lambda/2}\int_\mathbb R\varphi(\kappa+\tau x)dx=\frac\sigma{\sqrt{2\pi}}e^{-\lambda/2}\frac1{\tau}$$ By identification, $$\tau^2=1+\sigma^2\qquad\kappa=\frac{\mu\sigma}\tau\qquad\lambda=\mu^2-\kappa^2=\frac{\mu^2}{1+\sigma^2}$$ hence

$$J_{\mu,\sigma^2}=\frac{\sigma e^{-\mu^2/(2(1+\sigma^2))}}{\sqrt{2\pi(1+\sigma^2)}}$$

A: The first integral is solved in Distribution of the normal cdf.
So we know that
$$ \int_{-\infty}^{\infty} \Phi \left (x  \right ) N \left ( x|\mu,\sigma^2 \right ) dx = \Phi \left (\frac{\mu}{\sqrt{1+\sigma^2}}  \right ) $$
Now, if we differentiate both sides w.r.t. $\mu$,
$$\int_{-\infty}^{\infty} \frac{x-\mu}{\sigma^2} \Phi \left (x  \right ) N \left ( x|\mu,\sigma^2 \right ) dx = \frac{1}{\sqrt{1+\sigma^2}} N \left (\frac{\mu}{\sqrt{1+\sigma^2}}  \right ) $$
Notice the first term on the left hand side is what we want. Rearranging,
$$\int_{-\infty}^{\infty} \Phi \left (x  \right ) N \left ( x|\mu,\sigma^2 \right ) x dx = \mu \Phi \left ( \frac{\mu}{\sqrt{1+\sigma^2}} \right ) + \frac{\sigma^2}{\sqrt{1+\sigma^2}} N \left (\frac{\mu}{\sqrt{1+\sigma^2}}  \right ) $$
