gaussian integral of power of cdf : $\int_{-\infty}^{+\infty} \Phi(x)^n \cdot \phi(a+bx) \cdot dx$

Question

Is there an analytic solution for the following Gaussian integral?

$$\int_{-\infty}^{+\infty} \Phi(x)^n \cdot \phi(a+bx) \cdot dx$$

with

• $n$, a positive integer (typically under 10)
• $a,b$, real numbers (typical values: $a$ between 1 and 30, and $b$ between 1 and 10)
• $\Phi(\cdot)$, the standard normal cumulative distribution function
• $\phi(\cdot)$, the standard normal density function

I found a solution for $n=1$ and $n=2$ (see http://en.wikipedia.org/wiki/List_of_integrals_of_Gaussian_functions)

However, I would need a general solution (for any $n$) if it exists. If not, is there a good approximation?

Thanks.

Answer 1

If we assume that $b$ is large, we can try the Laplace method see http://en.wikipedia.org/wiki/Laplace%27s_method look under "Other formulations". This assumes that the integrand is maximal at $-a/b$ and decays fast around this point. This gives $$ \int \Phi^n(x)\phi(bx+a) dx \approx \Phi^n(-a/b)\int \phi(bx+a) dx=$$ $$\Phi^n(-a/b) \frac{1}{b}$$ This is now the quick and dirty proposal. If it doesn't help, (please some information about the context), maybe one can try Laplace's method in some other way, i.e. with $n$ as large parameter.

Answer 2

Following the methods explained in the post integrating-pdf-times-cdf-squared (restates the integral in terms of probability arguments) and redefining the same integral we obtain, \begin{eqnarray} \int_{-\infty}^{\infty}\Phi(x)^{n}\phi(a+bx)dx&=&\frac{1}{b}\int_{-\infty}^{\infty}\Phi\left(\frac{y-a}{b}\right)^{n}\phi(y)dy\\&=&\frac{1}{b}P\left(z_{1}\leq \frac{y-a}{b},z_{2}\leq \frac{y-a}{b},z_{3}\leq \frac{y-a}{b}....z_{n}\leq \frac{y-a}{b}\right)\\&=&\frac{1}{b}P\left(bz_{1}-y\leq -a,bz_{2}-y\leq -a,bz_{3}-y\leq -a....bz_{n}-y\leq -a\right)\\&=& \frac{1}{b}\mathcal{MVN}\left(x=\{-a,-a,-a...\},\mu=\{0,0,0...\},\Sigma=\begin{bmatrix}b^{2}+1 & 1 & 1&...\\1&1&b^{2}+1&1&...\\1&1&1&b^{2}+1&1... \end{bmatrix} \right) \end{eqnarray} where $\mathcal{MVN}$ is multivariatenormal distribution CDF. An R-based quick application for a=1,b=1, n=3 and a=1,b=2,n=4 shows that indeed the integral in question can be computed using this solution, we can observe this from coinciding results of numerical integral using the definition and the multivariate normal CDF. The only needed effort is to define the covariance matrix.

 sg<-function(a,b,x){pnorm(x)^3*dnorm(a+b*x)}
 integrate(sg,-Inf,Inf,a=1.,b=1.)
0.06832488 with absolute error < 4.9e-07

pmnorm(x = c(-1.,-1.,-1),mean = rep(0.,3),varcov = matrix(c(2,1.,1.,1.,2,1.,1,1,2),ncol=3,byrow=T))
[1] 0.06832488

sg<-function(a,b,x){pnorm(x)^4*dnorm(a+b*x)}
> integrate(sg,-Inf,Inf,a=1.,b=2.)
0.01677928 with absolute error < 3.3e-05

pmnorm(x = c(-1.,-1.,-1,-1),mean = rep(0.,4),varcov = matrix(c(5,1.,1.,1.,1,5,1.,1,1,1,5,1,1,1,1,5),ncol=4,byrow=T))/2.
[1] 0.0167786