Boosting the performance of expensive NIntegrate by feeding in a cheap approximation of the integrand

Question

I need to integrate an expensive likelihood L[x] over its n-dimensional domain. I know that L[x] is decently approximated by a gaussian likelihood G[x], which is very cheap to evaluate. In particular, L and G share the same maximum. The approximation, however, is not good for the integral I need.

I've already got a good performance boost by restricting the integration to the 4-sigma Ellipsoid predicted by G, as discussed in NIntegrating within an Ellipsoid.

I would like to ask now if it is possible to feed NIntegrate with the information that G is a decent approximation for L so that the sampling in NIntegrate is done in a more efficient way.

Answer 1

Would it be a solution to evaluate $L$ at a set of points $x_i$ and morph $G$ with a function $G_m$ in such a way that $G_m[x]G[x]=L[x]$ for those points? That is, with $G_m[x]$ defined as the interpolated function through the set of points ${x_i,L[x_i]/G[x_i]}$.

For example, let your two functions be:

l[x_] := PDF[GammaDistribution[5, 2], x]
g[x_] := PDF[NormalDistribution[8, (3 E^4)/(16 Sqrt[2 \[Pi]])], x]

They share the same maximum, and are somewhat similar but not closely:

Plot[{l[x], g[x]}, {x, 1, 20}, Frame -> True, Axes -> False]

Then with:

gm = Interpolation[Table[l[x]/g[x], {x, 1, 21}]];

the plots are much more similar:

Plot[{l[x], gm[x] g[x]}, {x, 1, 20}, Frame -> True, Axes -> False]

NIntegrate[l[x], {x, 1, 20}]

(*
==> 0.9705751963
*)

NIntegrate[g[x], {x, 1, 20}]

(*
==> 0.9550846595
*)

NIntegrate[gm[x] g[x], {x, 1, 20}]

(*
==> 0.9703985212
*)

In this case, the result of integrating the morphed 'el cheapo' function is within 0.02% of the value of the 'expensive' function. Of course, for this to work the functions should be sufficiently smooth and G should never get too close to zero.

The set $x_i$ could be obtained with the EvaluationMonitor option in the integration of $G$:

 Last@ Reap[
   NIntegrate[g[x], {x, 1, 20}, EvaluationMonitor :> Sow[x]]]

    (*
    ==> {{1.151189079, 1.891291464, 3.335416098, 5.384541554, 7.843511075, 
  10.5, 13.15648893, 15.61545845, 17.6645839, 19.10870854, 
  19.84881092, 1.07559454, 1.445645732, 2.167708049, 3.192270777, 
  4.421755537, 5.75, 7.078244463, 8.307729223, 9.332291951, 
  10.05435427, 10.42440546, 10.57559454, 10.94564573, 11.66770805, 
  12.69227078, 13.92175554, 15.25, 16.57824446, 17.80772922, 
  18.83229195, 19.55435427, 19.92440546}}
    *)

This set is for one reason or another a meaningful set and could be used as anchorpoints for the interpolating function.

---

Actually, now that I come to think about it a bit more, this may not be so useful in many cases as the number of evaluations of the expensive function may not be drastically reduced. We need to evaluate the sample points and this set may be as large as the set NIntegrate might need. So, this saves time only if we can reduce the number of points used in the interpolation function. But if we could halve that set this would save a factor $2^n$ for $n$-dimensional functions.