dblquad integral strange behaviour

Hello,

Following is my function (which approximately calculates overlap of 2 bivariate distributions with identity covariance matrices)

function [ res ] = bivariate_overlap_integral(mu_x1,mu_y1,mu_x2,mu_y2)
bpdf_vec1=@(x,y,mu_x,mu_y)(exp(-((x-mu_x).^2)./2.-((y-mu_y)^2)/2)./(2*pi));

overlap_point = @(x,y) min(bpdf_vec1(x,y,mu_x1,mu_y1),bpdf_vec1(x,y,mu_x2,mu_y2));

res=dblquad(overlap_point,-100,100,-100,100);

You can see that this involves taking a double integral (x: -100 to 100, y:-100 to 100, ideally -inf to inf but suffices at the moment) from the function overlap_point which is minimum of 2 pdf-s given by the function bpdf_vec1 of two distributions at the point x,y.

Now, PDF is never 0, so I would expect that the larger the area of the interval, the larger will the end result become, obviously with the negligible difference after a certain point. However, it appears that sometimes, when I decrease the size of the interval the result grows. For instance

>> mu_x1=0;mu_y1=0;mu_x2=5;mu_y2=0;
>> bpdf_vec1=@(x,y,mu_x,mu_y)(exp(-((x-mu_x).^2)./2.-((y-mu_y)^2)/2)./(2*pi));
>>  overlap_point = @(x,y) min(bpdf_vec1(x,y,mu_x1,mu_y1),bpdf_vec1(x,y,mu_x2,mu_y2));

>> dblquad(overlap_point,-10,10,-10,10)

ans =

    0.0124
>> dblquad(overlap_point,-100,100,-100,100)

ans =

1.4976e-005  -----> strange, as theoretically cannot be larger thn the first answer

>> dblquad(overlap_point,-3,3,-3,3)

ans =

    0.0110  -----> makes sense that the result is less than the first answer as the interval is decreased

Does this have to do with the implementation of dblquad, or am I making a mistake somewhere? I use MATLAB R2011a

thanks

Answer 1

[Mike Hosea]

Since adaptive quadrature relies on finite sampling, and the number of initial points tends to be fixed, regardless of interval size, it turns out that increasing the bounds with a decaying function is not a reliable way of approximating an improper integral. This is true of all the adaptive quadrature routines I have used.

DBLQUAD is deprecated now. You should be using INTEGRAL2 with infinite limits. If you want to use finite limits in a situation like this, break the integral up into smaller chunks to help ensure that the integrator "sees" the variation of the function.