0

Sometimes scipy.integrate.quad wrongly returns near-0 values. This has been addressed in this question, and seems to happen when the integration technique doesn't evaluate the function in the narrow range where it is significantly different than 0. In this and similar questions, the accepted solution was always to use the points parameter to tell scipy where to look. However, for me, this seems to actually make things worse.

Integration of exponential distribution pdf (answer should be just under 1):

from scipy import integrate
import numpy as np
t2=.01
#initial problem: fails when f is large
f=5000000
integrate.quad(lambda t:f*np.exp(-f*(t2-t)),0,t2)
#>>>(3.8816838175855493e-22, 7.717972744727115e-22)

Now, the "fix" makes it fail, even on smaller values of f where the original worked:

f=2000000
integrate.quad(lambda t:f*np.exp(-f*(t2-t)),0,t2)
#>>>(1.00000000000143, 1.6485317987792634e-14)
integrate.quad(lambda t:f*np.exp(-f*(t2-t)),0,t2,points=[t2])
#>>>(1.6117047218907458e-17, 3.2045611390981406e-17)
integrate.quad(lambda t:f*np.exp(-f*(t2-t)),0,t2,points=[t2,t2])
#>>>(1.6117047218907458e-17, 3.2045611390981406e-17)

What's going on here? How can I tell scipy what to do so that this will evaluate for arbitrary values of f?

Kalev Maricq
  • 617
  • 1
  • 7
  • 24
  • 1
    Try an array of points just under `t2`, e.g. `x = np.linspace(t2-.00001,t2,40)` – hpaulj May 03 '22 at 16:55
  • @hpaulj Any idea how to generalize that for different values of t2 and f? For example, if I try t2=1 and f=20000, then it stops working. – Kalev Maricq May 03 '22 at 17:44
  • All I did was to create a set of points where the function value was noticeably more than 0. I based that on a plot of the function, using different ranges. It has a limit of 50 points. – hpaulj May 03 '22 at 18:46
  • Ok, thanks. For my use case, I won't know ahead of time what t2 and f are. – Kalev Maricq May 03 '22 at 22:03
  • I don't know if there's a integrator that's more careful about identifying places that need greater sampling. Review the docs. Sometimes, even often, the user has to know something about their problem, and anticpate problems. Packages can only put so much effort into looking for problem cases. – hpaulj May 03 '22 at 22:10

1 Answers1

0

It's not a generic solution, but I've been able to fix this for the given function by using points=np.geomspace to guide the numerical algorithm towards heavier sampling in the interesting region. I'll leave this open for a little bit to see if anyone finds a generic solution.

Generate random values for t2 and f, then check min and max values for subset that should be very close to 1:

>>> t2s=np.exp((np.random.rand(20)-.5)*10)
>>> fs=np.exp((np.random.rand(20)-.1)*20)
>>> min(integrate.quad(lambda t:f*np.exp(-f*(t2-t)),0,t2,points=t2-np.geomspace(1/f,t2,40))[0] for f in fs for t2 in t2s if f>(1/t2)*10)
0.9999621825009719
>>> max(integrate.quad(lambda t:f*np.exp(-f*(t2-t)),0,t2,points=t2-np.geomspace(1/f,t2,40))[0] for f in fs for t2 in t2s if f>(1/t2)*10)
1.000000288722783
Kalev Maricq
  • 617
  • 1
  • 7
  • 24