MATLAB: Strange behaviour in integral function in MATLAB

Question

I create two functions fx and fy as follows:

 fx = @(x) ncx2pdf(x, 10, 1);
 fy = @(y) ncx2pdf(y + 25, 10, 10);

Then, I define fs function as follows:

 shift = 0
 fs = @(s) fx(s) .* fy(shift - s)

Note that fs is always positive (product of two probability density function). If I compute:

 integral(fs, -Inf, 100)

I obtain the true value 0.0413, but If I compute

 integral(fs, -Inf, 1000)

I obtain 0. Why this strange behavior happens using integral function? Note that if I compute

 integral(fs, -Inf, Inf)

I obtain the true value 0.0413.

Answer 1

Numeric quadrature works by iteratively splitting the interval apart and approximating the integral of each partition. The approximation is done by discretizing the subintervals and using some form of integral approximation (Riemann sums etc.). Since there is some discretization, there will be some error and this is dependent on the partitioning, therefore we are best to keep our limits of integration closer to the points which will contribute most to the integral.

Since you know this function is approximately zero almost everywhere and has its largest weight around the origin (relatively speaking), the integral which samples points closer to the origin will perform better. For this reason it is a good idea to avoid using Inf and -Inf in this example. If however you really want to use these extremes, then it is best to use the Waypoints to let MATLAB know where it should be splitting the integral. Either of these will produce the results you are looking for:

integral(fs,-10000,10000) % Even this far away from the largest function values.
integral(fs,-Inf,1000,'Waypoints',-50:50) % MATLAB will split integral at these locations