Got Questions? Get Answers.
Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Bug with dblquad

Subject: Bug with dblquad

From: Guillaume

Date: 30 Jun, 2009 09:58:33

Message: 1 of 6

Hi everyone,

I got to compute a highly complex double integration, which gives erroneous results. That's why I decided to verify if dblquad was working well. To proof this, I decided to compute the integral of the gaussian probability density :

fufunc = @(x,y)((1/sqrt(2*pi)).*exp(-0.5*(x).^2));
result = quad2d(fufunc,-1000,1000,-1000,1000)/2000

result = 1.0 as expected

But then if I try to move the distribution by changing the mean:

fufunc = @(x,y)((1/sqrt(2*pi)).*exp(-0.5*(x-10).^2));
result = quad2d(fufunc,-1000,1000,-1000,1000)/2000

result = 9.3739e-21.

instead of 1 ...

Does anyone has an idea, where the problem comes from?

Subject: Bug by integrating with dblquad AND QUAD

From: Guillaume

Date: 30 Jun, 2009 10:32:01

Message: 2 of 6

I just notticed that the quad function produce the same error, while a standard trapz integration works fine....

 fufunc2 = @(x)((1/sqrt(2*pi)).*exp(-0.5*(x-10).^2));

axe=-100:0.01:100;
result=trapz(axe, ((1/sqrt(2*pi)).*exp(-0.5*(axe-10).^2)))
resultat2 = quad(fufunc2,-1000,1000)

result =
     1

resultat2 =
   3.5175e-20

"Guillaume " <fufubynight@free.fr> wrote in message <h2cnk9$a3v$1@fred.mathworks.com>...
> Hi everyone,
>
> I got to compute a highly complex double integration, which gives erroneous results. That's why I decided to verify if dblquad was working well. To proof this, I decided to compute the integral of the gaussian probability density :
>
> fufunc = @(x,y)((1/sqrt(2*pi)).*exp(-0.5*(x).^2));
> result = quad2d(fufunc,-1000,1000,-1000,1000)/2000
>
> result = 1.0 as expected
>
> But then if I try to move the distribution by changing the mean:
>
> fufunc = @(x,y)((1/sqrt(2*pi)).*exp(-0.5*(x-10).^2));
> result = quad2d(fufunc,-1000,1000,-1000,1000)/2000
>
> result = 9.3739e-21.
>
> instead of 1 ...
>
> Does anyone has an idea, where the problem comes from?

Subject: Bug by integrating with dblquad AND QUAD

From: John D'Errico

Date: 30 Jun, 2009 11:18:01

Message: 3 of 6

"Guillaume " <fufubynight@free.fr> wrote in message <h2cpj1$6vl$1@fred.mathworks.com>...
> I just notticed that the quad function produce the same error, while a standard trapz integration works fine....
>
> fufunc2 = @(x)((1/sqrt(2*pi)).*exp(-0.5*(x-10).^2));
>
> axe=-100:0.01:100;
> result=trapz(axe, ((1/sqrt(2*pi)).*exp(-0.5*(axe-10).^2)))
> resultat2 = quad(fufunc2,-1000,1000)
>
> result =
> 1
>
> resultat2 =
> 3.5175e-20

Give me a single penny for every time that someone
has found a "bug" in a numerical integration routine
for exactly this, and I would be wealthy beyond belief.

I.e., they have tried to do a numerical integration over
a function that is essentially zero almost everywhere
relative to the width of the interval, and the numerical
integration returns the wrong result.

When you integrate a function like the PDF of a normal
distribution over a very wide domain, quad finds that
it is essentially zero at every point that it samples the
function at. So it returns zero. It does not know that
you have given it a black box that is an essential delta
function, so a spike in some position.

When you used a trapezoidal rule, you used a fine
enough interval to evaluate the function 10000 times.
You manage to successfully catch the region where
the PDF has a non-zero value. But this many function
evaluations is far too many for quad to make every
time you call it. And don't forget that dblquad would
then be forced to make 10000^2 = 1e8 function
evaluations to have equal probability of success.

John

Subject: Bug by integrating with dblquad AND QUAD

From: Guillaume

Date: 30 Jun, 2009 12:25:02

Message: 4 of 6

"John D'Errico" <woodchips@rochester.rr.com> wrote in message <h2cs99$sau$1@fred.mathworks.com>...
> Give me a single penny for every time that someone
> has found a "bug" in a numerical integration routine
> for exactly this, and I would be wealthy beyond belief.
>
> I.e., they have tried to do a numerical integration over
> a function that is essentially zero almost everywhere
> relative to the width of the interval, and the numerical
> integration returns the wrong result.
>
> When you integrate a function like the PDF of a normal
> distribution over a very wide domain, quad finds that
> it is essentially zero at every point that it samples the
> function at. So it returns zero. It does not know that
> you have given it a black box that is an essential delta
> function, so a spike in some position.
>
> When you used a trapezoidal rule, you used a fine
> enough interval to evaluate the function 10000 times.
> You manage to successfully catch the region where
> the PDF has a non-zero value. But this many function
> evaluations is far too many for quad to make every
> time you call it. And don't forget that dblquad would
> then be forced to make 10000^2 = 1e8 function
> evaluations to have equal probability of success.
>
> John

Thanks for this answer. I reduced the integration domain and it works. I was wrong by thinking that the bigger the domain is, more accurate is the result...

The thing is that if my integration domain is not well selected, the computation... What a chance ! Furthermore, the double integration is within a loop and the gaussian distribution moves along the axes...

I can't either find where the spike is and compute around, or use the trapz method with a double integration...

Subject: Bug by integrating with dblquad AND QUAD

From: John D'Errico

Date: 30 Jun, 2009 12:57:01

Message: 5 of 6

"Guillaume " <fufubynight@free.fr> wrote in message <h2d06u$kk0$1@fred.mathworks.com>...

> Thanks for this answer. I reduced the integration domain and it works. I was wrong by thinking that the bigger the domain is, more accurate is the result...
>
> The thing is that if my integration domain is not well selected, the computation... What a chance ! Furthermore, the double integration is within a loop and the gaussian distribution moves along the axes...
>
> I can't either find where the spike is and compute around, or use the trapz method with a double integration...

What I forgot to follow up with is that you can easily do
normal PDF computations with a call to erf possibly using
a transformation. erf will be correct, as it is not an adaptive
integration rule.

John

Subject: Bug by integrating with dblquad AND QUAD

From: Steven Lord

Date: 30 Jun, 2009 13:32:46

Message: 6 of 6


"Guillaume " <fufubynight@free.fr> wrote in message
news:h2d06u$kk0$1@fred.mathworks.com...
> "John D'Errico" <woodchips@rochester.rr.com> wrote in message
> <h2cs99$sau$1@fred.mathworks.com>...

*snip*

> Thanks for this answer. I reduced the integration domain and it works. I
> was wrong by thinking that the bigger the domain is, more accurate is the
> result...

Nope. The bigger the region over which you integrate, the easier it is to
miss the interesting details. Fly over the town/city where your house is
located at 30,000 feet and try to identify your house. Now fly over the
town/city at 1,000 feet and try to identify your house -- you're going to be
able to identify your house much more precisely, aren't you?

> The thing is that if my integration domain is not well selected, the
> computation...

fails to give the answer you expect. Yes. You need to solve the correct
problem the correct way to get the results you expect, and this isn't
limited just to numeric computations.

>What a chance ! Furthermore, the double integration is within a loop and
>the gaussian distribution moves along the axes...
>
> I can't either find where the spike is and compute around, or use the
> trapz method with a double integration...

Looking back at your original posting, the function you're integrating has
no dependence on y at all, so instead of using QUAD2D or DBLQUAD, you should
be using QUADGK with appropriate limits of integration. Figuring out the
right limits of integration is the key step here -- perhaps evaluating the
integrand at each integer until you get a result that's small enough that
you can safely neglect the area under the function in that 'tail' region
past that point.

--
Steve Lord
slord@mathworks.com

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us