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)?