Scipy integral comes out wrong

Question

I'm trying to integrate a function. The function is guaranteed to be non negative:

def function(x):
    something = ...
    something_else = ...
    return exp(something) / sqrt(something_else)

Now I'm integrating it:

def integrand(a, b):
    return scipy.integrate.quad(function, a, b)

The results I'm getting are not what I expect, so to check I did this:

for x in range(0, 10000):
    if integrand(0,x+1) < integrand(0,x):
        raise ValueError("Weird!")

Sure enough, I'm getting 'Weird' exceptions. How can that be?

Can you show an actual example that demonstrates the problem? — BrenBarn, Nov 09 '14 at 05:38
It's kind of hard, because the actual example involves a 10,000-point C14 dating calibration curve... — zmbq, Nov 09 '14 at 05:39
Well, it's also kind of hard to debug the problem without a small test case. Try to reduce the problem to a small test case. — BrenBarn, Nov 09 '14 at 05:52
Try doing a check of `function(0..10000)` for `< 0` to make 100% sure that it _never_ returns -ve numbers. — Steve Barnes, Nov 09 '14 at 06:25
it hopefully shouldn't - from the very vague function definition one can still gather that `exp(something) / sqrt(something_else) > 0` if something positive is fed into `something_else`. But this type of guessing game is usually a waste of time - more info please! — eickenberg, Nov 09 '14 at 09:06
I made sure that function is non-negative for range(0,10000,0.01). I'll try to find a small self-contained example for this. I thought this was a known issue with the way quad estimates the integral. — zmbq, Nov 09 '14 at 09:25

Warren Weckesser · Answer 1 · 2014-11-09T17:42:35.210

Without a working example, it is difficult to debug your code, but here are some suggestions based on what you have shown.

quad returns a tuple with two values: the estimate of the integral, and an estimate of the absolute error in the result. You don't appear to be using just the first value in your code. Try changing integrand to

def integrand(a, b):
    return scipy.integrate.quad(function, a, b)[0]

Note the [0] after the call to quad.

Next, you should print the values of integrand(0, x) and integrand(0, x+1) when you find a "weird" case. If they are approximately the same, then normal numerical imprecision could destroy the order that you would theoretically expect. Take a look at the second value returned by quad. If that value is bigger than the difference between integrand(0, x) and integrand(0, x+1), you can't really expect the estimates of the integrals to increase monotonically, because the expected increase is less than the error in the calculation. If that is the problem, you could try setting epsabs and/or epsrel to values smaller than the default values of 1.49e-8, but that will only get you so far. Eventually you run into the limits of floating point calculations.

Scipy integral comes out wrong

1 Answers1