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.

  • $\begingroup$ Have you seen the Method option? Maybe it can help, but I'm not sure it's enough for your case. $\endgroup$
    – FJRA
    Mar 14, 2012 at 13:56
  • $\begingroup$ Using StepMonitor you could try to find out what steps are used in the integration of G, but I don't have a clue how that would help in the integration of L. I haven't found a method in the plethora of methods that NIntegrate has that lets you specify a sampling grid. $\endgroup$ Mar 14, 2012 at 14:42
  • $\begingroup$ Is it possible to reduce you integration to root finding?.. And maybe it is a good idea to cross-post this question, perhaps in slightly modified form, to scicomp.stackexchange.com. $\endgroup$
    – faleichik
    Mar 14, 2012 at 16:01
  • 3
    $\begingroup$ @faleichik Cross posting is not allowed on SE, or at least heavily discouraged. One of the two would be closed. Let's leave it here for now, this is a very good on-topic question here. $\endgroup$
    – Szabolcs
    Mar 14, 2012 at 16:36
  • 2
    $\begingroup$ Chances are sampling is relatively heavy near the max value. Maaybe cut out a box around that, integrate L[x] elsewhere in your region, and integrate G[x] in the cut out box? Another possibility of course is to do an Interpolation of L and numerically integrate that instead of L. $\endgroup$ Mar 14, 2012 at 16:56

1 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]

Mathematica graphics

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]

Mathematica graphics

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.


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.