slwcon and eincr in integrating multivariate Gaussian random variables

Question

Why is the following integration involving multivariate Gaussian distribution so slow and generating an error? Is there a better integration strategy? All that I'm doing is considering $(X_1,X_2,X_3,X_4) \sim N(\mu, \Sigma)$ and computing the (trivial expectation) $E[1]$, which obviously equals to $1$.

muvec = ConstantArray[0.1, 4];
sigmat = IdentityMatrix[4];
npdf[x_] := PDF[MultinormalDistribution[muvec, sigmat], x]; 
NIntegrate[ 1 * npdf[{x1, x2, x3, x4}], 
           {x1, -Infinity, Infinity}, 
           {x2, -Infinity, Infinity}, 
           {x3, -Infinity, Infinity}, 
           {x4, -Infinity, Infinity} ] // AbsoluteTiming

This results in the errors

NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small. >>

NIntegrate::eincr: The global error of the strategy GlobalAdaptive has increased more than 2000 times. The global error is expected to decrease monotonically after a number of integrand evaluations. Suspect one of the following: the working precision is insufficient for the specified precision goal; the integrand is highly oscillatory or it is not a (piecewise) smooth function; or the true value of the integral is 0. Increasing the value of the GlobalAdaptive option MaxErrorIncreases might lead to a convergent numerical integration. NIntegrate obtained 1.000000567592556` and 0.000021505642234874004` for the integral and error estimates. >>

and the result

{12.521, 1.}

The result is correct, but it has taken far too long (in my opinion), and moreover, I'm a little bit puzzled by the error message (despite reading http://reference.wolfram.com/language/tutorial/NIntegrateIntegrationStrategies.html), and I've tried other integration methods and it's still relatively slow.

In the actual application I have in mind, I want to consider $h(y) := E[ g(y;X_1,X_2,X_3,X_4)] $, where $g$ is somewhat complex and I need to use the expression $h(y)$ over and over again, and hence if I can't even compute $E[1]$ quickly and accurately, I have little hopes for computing the more difficult expectation.

Edit I emphasize that a symbolic integration of the solution is not of interest. In particular, since the question at hand is a "toy example", of which the real application is computing a more complex expectation $E[g(X_1,\ldots, X_n)]$ of which the function $g$ is sufficiently complex that there's no hope of a symbolic integral solution.

Answer 1

Because the integrand is highly localized at the origin; e.g.,

Log[10, npdf[{6, 0, 0, 0}]] // N
(* -9.16177 *)

limiting the range of integration reduces run time by a factor of three and moderately improves accuracy.

muvec = ConstantArray[1/10, 4];
sigmat = IdentityMatrix[4];
npdf[x_] := PDF[MultinormalDistribution[muvec, sigmat], x];
NIntegrate[1*npdf[{x1, x2, x3, x4}], {x1, -6, 6}, {x2, -6, 6}, {x3, -6, 6}, 
    {x4, -6, 6}, PrecisionGoal -> 5] // AbsoluteTiming

Not nearly as good as the answer by ciao but perhaps helpful. (Note that only the error message NIntegrate::slwcon persists.)

By the way, I first tried several NIntegrate methods but with little improvement. The issue, it seems to me, is that the integrand becomes ever more localized as the dimension increases. 3D works fine, for instance. Also, if a Precision of, say, 3 is acceptable, big savings can be achieved: {0.275684, 1.00009}. Decreasing the range of integration in this case helps further, because it concentrates the integration more on the central peak: {0.241199, 1.00001}. Taking advantages of symmetries, if any, also could be helpful. Finally, a coordinate transformation to spread the central peak relative to the wings could be helpful.

Answer 2

While the unknown g might or might not be amenable to this approach, it's really fast on g == 1:

Expectation[
  1, {x1, x2, x3, x4} \[Distributed] 
   MultinormalDistribution[muvec, sigmat]] // RepeatedTiming

(23000 times faster than NIntegrate.)

Then there's also NExpectation to try, too.

---

Ever since I came across Guess who it is's use of the Chebyshev polynomial proxy this morning in How to find numerically all roots of a function in a given range?, I've been itching for something to try it out on. It's really cool. It approximates a univariate function, but chebfun has a bivariate method; perhaps it is extensible. What I want to show here is how well this works on the following.

We approximate g over the interval from -10 to 10; outside that the Gaussian is so small that practically there is no contribution to the integral.

r = 10;
g = BesselJ[2, #] &;
n = 32;
cnodes = Rescale[N[Cos[Pi Range[0, n]/n], 30], {-1, 1}, {-r, r}];
cc = Sqrt[2/n] FourierDCT[g /@ cnodes, 1];
cc[[{1, -1}]] /= 2;

gcheb = Expand[cc.Table[ChebyshevT[n - 1, x1/r], {n, Length@cc}]] /. 
   c_ /; c == 0 :> 0;
Expectation[gcheb,
   {x1, x2, x3, x4} \[Distributed] 
   MultinormalDistribution[muvec, sigmat]] // AbsoluteTiming

(*  {0.081965, 0.0988631125184331759104}  *)

There are coefficients of the form 0``23.443340842020458, and it saves time with the integration to remove them; hence the replacement c_ /; c == 0 :> 0. The use of n = 32 was inherited from @Guess; the precision of cnodes was adjusted from 20 to 30 after examining gcheb, but it made little difference in the result.

The first NIntegrate misses the mark widely. Increasing MinRecursion helps, but it takes much, much longer, even though it recognizes a convergence problem.

NIntegrate[
  BesselJ[2, x1]*npdf[{x1, x2, x3, x4}], {x1, -Infinity, 
   Infinity}, {x2, -Infinity, Infinity}, {x3, -Infinity, 
   Infinity}, {x4, -Infinity, Infinity}] // AbsoluteTiming

(*  {105.429, 0.00890915}  *)

NIntegrate[
  BesselJ[2, x1]*npdf[{x1, x2, x3, x4}], {x1, -Infinity, 
   Infinity}, {x2, -Infinity, Infinity}, {x3, -Infinity, 
   Infinity}, {x4, -Infinity, Infinity}, 
  MinRecursion -> 3] // AbsoluteTiming

NIntegrate::eincr:...

(*  {1563.69, 0.0988631}  *)

Answer 3

About 30X faster integration, no errors...

Clear[x, x1, x2, x3, x4]
muvec = ConstantArray[1/10, 4];
sigmat = IdentityMatrix[4];
pdf = PDF[MultinormalDistribution[muvec, sigmat], {x1, x2, x3, x4}] //
   FullSimplify
Integrate[
  1*pdf, {x1, -Infinity, Infinity}, {x2, -Infinity, Infinity}, 
         {x3, -Infinity, Infinity}, {x4, -Infinity, Infinity}] // AbsoluteTiming