10
$\begingroup$

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:

  1. Can this integral be further simplified?
  2. 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)?
Improve this question
$\endgroup$
6
  • 1
    $\begingroup$ Some questions and remarks: Note that there is a high degree of symmetry here; in particular, the implicit joint probability is invariant wrt to $n+1$ iid standard normal random variables. Are there any restrictions on $a_i$ and $b_i$, or can they vary arbitrarily? Over what range of error bounds are you interested in? Are you looking for hard bounds or are you potentially content with "softer" statistical ones, e.g., via some Monte Carlo approach? $\endgroup$
    – cardinal
    Apr 4, 2013 at 3:13
  • $\begingroup$ Cardinal, you're correct. In particular, I'm trying to estimate the probability that the value $(N+1)^{th}$ independent normal variable exceeds that of the other $N$ variables. Accordingly, $a_i$ is restricted to positive real numbers. Soft bounds, including Monte Carlo techniques are fine, though I'd prefer a deterministic solution. $\endgroup$
    – melchimm
    Apr 4, 2013 at 3:48
  • $\begingroup$ Why is it important to minimize the evaluations of $\Phi$? Isn't evaluating $\Phi$ only slightly slower than evaluating an elementary function? $\endgroup$
    – Douglas Zare
    Apr 4, 2013 at 4:57
  • $\begingroup$ Why calculating the integral and not just trying an acceptance-rejection type approach? - Simulating the $(n+1)$st variable and comparing it with other $n$ normal pdfs? This should be faster then integration. If you want hard bounds, you could try to some Quasi-Monte Carlo based sampling. [But this is maybe a bit trickier, as you have to keep the low discrepancy while transforming uniform distributions in normal ones...] $\endgroup$
    – Stephan Sturm
    Apr 4, 2013 at 5:13
  • $\begingroup$ @cardinal Could you be more explicit about the meaning of symmetry here? Why is it relevant for the computation? $\endgroup$
    – an12
    Apr 4, 2013 at 12:02

3 Answers 3

Reset to default
3
$\begingroup$

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.

Improve this answer
$\endgroup$
2
  • $\begingroup$ I should have mentioned this in my original question, but I am currently using a different quadrature rule that is inspired by Monte Carlo methods and that I believe (intuitively, I have no proof) should be more efficient than G-H. What I'm currently doing is choosing the abscissae to be centered on each of k quantile intervals for the normal pdf and giving the corresponding function evaluations equal weight in the resulting sum. Given that, as John points out below, my integrand is not well approximated by a polynomial, is there any reason to think that G-H would be more efficient? $\endgroup$
    – melchimm
    Apr 4, 2013 at 14:29
  • $\begingroup$ G-H integration with $n$ abscissae gives the exact integral of $e^{-x^2}p(x)$ when $p(x)$ is any polynomial of degree at most $2n-1$. So it gives a fair approximation when $p(x)$ can be fairly approximated by a polynomial of degree $2n-1$ over the region that matters, and $2n-1$ is a pretty high degree polynomial to play with. It will depend on the constants $a_i,b_i$. If they are all equal and $N$ is very large, the product approaches a step function and your integral is just a chunk of the guassian integral; G-H won't work in that case. $\endgroup$
    – Brendan McKay
    Apr 4, 2013 at 23:13
3
$\begingroup$

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.

Improve this answer
$\endgroup$
2
  • $\begingroup$ The important thing is whether the product term is well approximable by a polynomial in the central region of the guassian term. The tails don't matter as the guassian will kill their contribution. But I'm not at all confident it is true in this instance when $N$ is large. The product of $\Phi$s might look closer to a step function. Your suggestion of a transformation is a good one. $\endgroup$
    – Brendan McKay
    Apr 4, 2013 at 23:03
  • $\begingroup$ The higher order the Gauss-Hermite rule, the further out the tails matter. So one thing you can't do is just keep increasing the order of your rule hoping to get more accuracy. It's possible that increasing the order of the rule may not help and may even hurt. $\endgroup$
    – John D. Cook
    Apr 6, 2013 at 23:03
0
$\begingroup$

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.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.