Minimal numerical integration problem demonstration.sagews

%md
I have a problem with numerical integration in Sagemath.

The following exemple tries to use a (ridiculously oversimplified) example of a Uniform(-1,1) density, and to integrate it. (And yes, I know that the sought function is $\min(0,\max(\frac{x+1}{2},1))$, thank you very much...).

What bothers me is not the "interesting" part of the integral ; it's what happens in the "uninteresting" (constant) parts...

I have a problem with numerical integration in Sagemath.

The following exemple tries to use a (ridiculously oversimplified) example of a Uniform(-1,1) density, and to integrate it. (And yes, I know that the sought function is $\min(0,\max(\frac{x+1}{2},1))$, thank you very much...).

What bothers me is not the "interesting" part of the integral ; it's what happens in the "uninteresting" (constant) parts...

def dunif(x,a,b):
    if(x<a):return(0)
    elif(x>b):return(0)
    else:return(1/(b-a))

The function itsel plots (and therefore probably computes) fine.

plot(lambda x:dunif(x,-1,1),(x,-10,10),figsize=3)

However, its numerical integral has some serious problems :

plot(lambda x:numerical_integral(lambda t:dunif(t,-1,1),-oo,x)[0],(x,-10,10), figsize=3)

This is probably not an artifact of plot : "zooming" on a problematic region, justs zooms on the same problem (i. e. not dependent of choice of plotting points) :

plot(lambda x:numerical_integral(lambda t:dunif(t,-1,1),-oo,x)[0],(x,7,8), figsize=3)

This might originate in Maxima, which exhibits the same behaviour :

maxima("define(dunif(x,a,b),if x<a then 0 else if x>b then 0 else 1/(b-a));")
maxima("define(punif(x,a,b),quad_qagi(dunif(t,a,b),t,minf,x)[1]);")
def punif(x,a,b): return(maxima("punif("+str(x)+","+str(a)+","+str(b)+");"))
plot(lambda x:punif(x,-1,1),(x,-10,10),figsize=3)

dunif(x,a,b):=ifx<athen0else(ifx>bthen0else1/(b-a)) punif(x,a,b):=quad_qagi(ift<athen0else(ift>bthen0else1/(b-a)),t,minf,x,epsrel=1.0e-8,epsabs=0.0,limit=200)[1]

Let's try quad (from mpmath), a nice low-dimension quadrature function, used by sympy :

from mpmath import quad
plot(lambda x:RR(quad(lambda t:dunif(t,-1,1), (-oo, x))),(x,-10,10), figsize=3)

This becomes ridiculous. Let's zoom in the first problematic region...

plot(lambda x:RR(quad(lambda t:dunif(t,-1,1), (-oo, x))),(x,1,3), figsize=3)

This is becoming preposterous. Let's zoom again...

plot(lambda x:RR(quad(lambda t:dunif(t,-1,1), (-oo, x))),(x,2.85, 2.95), figsize=3)

Theorem: Numerical analysis in hard.

Corollary: Numerical analysis is harder in apparently simple cases...

Proofs: Plotting a constant function is not (not hard). an erroneous plot is not a not-hard plot. See above...

More seriously : What can be the source of the problem ? I am aware that I am using very general algorithms in a special case susceptible of many simplifications (starting with using the closed form of the integral...). But this is a test case, prelude to more complex cases with no obvious simplifications.

Any ideas ?

Well, let's try to force dunif to accept and return floating point values : that should take care of "surprising" conversions :

%cython
cdef double typed_dunif(double x, double a, double b):
    if(x<a):
        return(0.0)
    elif(x>b):
        return(0.0)
    else:
        return(1.0/(b-a))

No functions defined.

Auto-generated code...

typed_dunif(2,-1,-1)

Error in lines 1-1 Traceback (most recent call last): File "/projects/sage/sage-7.3/local/lib/python2.7/site-packages/smc_sagews/sage_server.py", line 957, in execute exec compile(block+'\n', '', 'single') in namespace, locals File "", line 1, in <module> NameError: name 'typed_dunif' is not defined

plot(lambda x:numerical_integral(lambda t:typed_dunif(t,-1,1),-oo,x)[0],(x,-10,10), figsize=3)

plot(lambda x:RR(quad(lambda t:typed_dunif(t,-1.0,1.0), (-oo, x))),(x,-10.0,10.0), figsize=3)

Same difference, again... The problem is not here.

Let's try R's integration routine :

%r
my_punif<-Vectorize(function(x,a,b)integrate(Vectorize(function(t)dunif(t,a,b)),-Inf,x)$value)
## Sanity check
my_punif(seq(-2,2,by=0.5),-1,1)

[1] 0.00 0.00 0.00 0.25 0.50 0.75 1.00 1.00 1.00

%r
curve(my_punif(x,-1,1),from=-10,to=10)

** Same problem ?!?!?! **

So far, we tried three implementations (Sage-Maxima, mpmath, R) of numerical integration routines on the same trivial example, which gave variations of the same error.

Either these routines implement variations of the same flawed algorithm, my proble is ill-thought or these routines use the same flawed piece of system software in some critical part of their algorithms...