Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
Have you seen the Method option? Maybe it can help, but I'm not sure it's enough for your case. –  FJRA Mar 14 '12 at 13:56
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. –  Sjoerd C. de Vries Mar 14 '12 at 14:42
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. –  faleichik Mar 14 '12 at 16:01
@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. –  Szabolcs Mar 14 '12 at 16:36
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. –  Daniel Lichtblau Mar 14 '12 at 16:56
add comment

1 Answer

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.

share|improve this answer
add comment

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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