Strange behaviour in integral function in MATLAB

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

[Brendan Hamm]

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

[Abdal-bassir Abou El-ailah]

Thanks for your response,

Normally, I do not know the support of the function to use Waypoints on the support (Above, I just gave a simple example). Is there an efficient function (way) in Matlab to estimate the support of a function (at precision 99.99 % of the energy for example).

[Brendan Hamm]

I cannot think of a specific way to ensure that the precision is 99.99% as this would require knowing the analytic solution of the integral. You could always perform a "sanity" check on a reasonable interval by using the fplot function and adjusting the bounds. For your example you could do something like:

figure
fplot(fs,-1000,1000)
figure
fplot(fs,-100,100)

[Abdal-bassir Abou El-ailah]

In fact, I am looking for an automatic way to compute this integral. I apply this procedure a lot of times. So, it is time-consuming to handle this by hand. What exactly I am doing is the convolution between two probability density functions (fx and fy), so if I know the support of fx and fy, I can estimte the support of the convolution between them. Normally, I use this function in Matlab for fx (fx = @(x) ncx2pdf(x, 10, 1)), but I did not find in Matlab a function that gives me the support of the PDF?

Thanks a lot for your efforts

[Brendan Hamm]

This is providing a non-central Chi-squared distribution which has support [0,+Inf). This function will return the value of 0 if you were to pass in a negative value as technically speaking the pdf is defined piece wise on the entire real line, whereby the pdf at negative values is 0.

[Abdal-bassir Abou El-ailah]

It is true but it is very large ([0, +Inf)), normally the most of the energy is somewhere depending on the parameters (degrees of freedom and non-central parameter). For example, technically speaking the support of a Gaussian distribution is ]-Inf, +Inf[, but 99.99 % of the energy is within [mu - 3*sigma, mu + 3*sigma]. I am trying to obtain a compact support in order to reduce the computational time of the convolution and the integral.

Thank you

[Brendan Hamm]

There is not going to be an easy way about this as the result of the convolution is not itself a pdf. One thing you may consider to decide the waypoints is the use of the icdf evaluated at the 99.99% level for each of the pdfs.

>> fzero(@(x) ncx2cdf(x,1,1e1)-0.9999,1e1)
ans =
   47.3522
>> fzero(@(x) ncx2cdf(x,1,1e2)-0.9999,1e2)
ans =
  188.2114
>> fzero(@(x) ncx2cdf(x,1,1e3)-0.9999,1e3)
ans =
   1.2490e+03
>> fzero(@(x) ncx2cdf(x,1,1e6)-0.9999,1e6)
ans =
   1.0075e+06

The effect of the dof parameter is shown below:

>> fzero(@(x) ncx2cdf(x,1e2,1e3)-0.9999,1e3)
ans =
   1.3536e+03
>> fzero(@(x) ncx2cdf(x,1e3,1e3)-0.9999,1e3)
ans =
   2.2995e+03
>> fzero(@(x) ncx2cdf(x,1e5,1e3)-0.9999,1e3)
ans =
   1.0269e+05

As you can see the dof parameter has a similar effect, so we could reasonable just add the two parameters together to get a reasonable upper bound on each of the pdfs separately then choose an appropriate interval. You should be careful as the shift you have on the y-value in fy will alter this region, but it is easy to take care of this:

dofx = 10;
nonx = 1;
fx = @(x) ncx2pdf(x,dofx,nonx);

% Note: fx(x) has support [0,Inf);

dofy = 10;
nony = 10;
yshift = 25;
fy = @(y) ncx2pdf(y + yshift,dofy,nony);

% Note fy(y) has support [-25,Inf) == [-yshift,Inf)

multFact = 10;
upperLimit = multFact * max(dofx+nonx,dofy+nony);

shift = 0 % If you change the shift the problem will change
fs = @(s) fx(s) .* fy(shift - s)
% fy(-y) has support (-Inf,25] so technically fs now has the support:
% [0,25] and we can approximate the integral with these bounds.

integral(fs,0,upperLimit,'Waypoints',linspace(0,upperLimit))

i am not sure how you plan on using the shift parameter, does this change?

Will you always have fs = @(s) fx(s) .* fy(shift-s)? If so, you should always be able to find the exact support of the function.

[Abdal-bassir Abou El-ailah]

It is a nice idea to use ncx2cdf to find the support at precision X. I implemented a simple approach based on the parameters to estimate the support of non-central Chi distribution, which gives sophisticated results. But now, for example, I can directly use fzero(@(x) ncx2cdf(x,1,1e1)- 0.9999,1e1) and fzero(@(x) ncx2cdf(x,1,1e1)- 0.001,1e1) to find the support at the precision that I want. Thanks for that.

Yes, I am always interested in finding the convolution of fx and fy. Now, I am trying to sample fx and fy and use directly 'conv' function in Matlab to estimate the convolution. Normally, the support of the convolution is [minFx + minFy, maxFx + maxFy], where [minFx, maxFx] and [minFy, maxFy] are the supports of fx and fy respectively (I do that to avoid the 'integral' function because it takes a lot of time, especially I must call it for each 'shift' value in the support).

By the way, in integral(fs,0,upperLimit,'Waypoints',linspace(0,upperLimit)), we can only use the support of fx (instead of [0, upperLimit]), but shift can vary in [0, upperLimit].

Thank you a lot.

[Brendan Hamm]

No Problem. Best of luck.