Ever since I came across Guess who it is's use of the Chebyshev polyonomial 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 coefficient 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} *)
