Rapid evaluation of multivariate normal integral I'm implementing a model that requires me to numerically evaluate a multivariate normal integral of the following form 
$$\int_{-\infty}^\infty \phi(z)\displaystyle\prod_{i=1}^N \Phi(a_iz+b_i) \, dz,$$
where $\phi(\cdot)$ and $\Phi(\cdot)$ represent the standard normal distribution function and integral, respectively. $N$ is a fairly large number and I need to be able to evaluate this integral rapidly. I have two questions:


*

*Can this integral be further simplified?

*What is the most efficient method for estimating this integral (e.g., requiring the fewest evaluations of $\Phi(\cdot)$ for an arbitrarily selected error bound)? 

 A: I suggest you try Gauss-Hermite integration. You can guess the precision by increasing the number of abscissae. Tables of abscissae and weights are here. 
A: Your integral gives the probability that a sample from a standard normal random variable is larger than the maximum of a set of normal random variables with means $b_i$ and standard deviations $a_i$. You could calculate the integral by simulation, drawing samples from these random variables and counting what proportion of the time the first sample is largest.
If N = 1, the integral can be evaluated in closed form. I don't know about larger N.
If you want to try Gauss-Hermite integration as suggested in another answer, do a change of variables first: Gauss-Hermite integration can work well with the right change of variables and terrible otherwise. Here's why: G-H assumes the integrand is of the form $\phi(x) P(x)$ where $\phi(x)$ is a normal PDF and $P(x)$ is a polynomial. The product term in your integrand is not even approximately like a polynomial because it is asymptotically 0 on the left and asymptotically 1 on the right. 
But if you find the mean and variance of the best normal approximation to your integrand, and do an affine change of variables so that your integrand is approximately a standard normal density, then G-H integration can work well.
A: I ended up using standard Gauss-Legendre quadrature for this problem, which provides a reasonably good approximation using relatively few abscissae. Brendan's solution (Gauss-Hermite) did a pretty horrible job in this case, for the reasons that John pointed out. I also tried using a transformation as suggested by John, but finding the appropriate change of variables (using the Laplace approximation) wound up being too computationally expensive.
