Absurd result using 'integral2' with non-"centered" functions

Hi,

I have been trying to compute an expression for a probability which involves four normal density functions ; those are by definition positive on all real numbers. That is, when integrating over any bounds the expression of a density function, we shouldn't get a smaller result by integrating over larger bounds.

Everything comes up to this problem : this code snippet given should yield 1 !

sigma = sqrt(6) ;
Theta = 50 ;                                        % Replace Theta by something less than 34 and it gives 1 !

f = @(x) normpdf(x,Theta,sigma) ;

integral2(@(x,y) f(x).*f(y),-inf,inf,-inf,inf)      % Integrating product of two normal probability density functions over R^2

ans = 1.9949e-30

This is not a problem with normpdf as I still observe the problem by replacing f by the litteral expression of the normal probability density function.

Am I missing something ?

Thanks for your time.

Answer 1

[Torsten]

It's not possible for integral2 to capture the small peak around 50 to return 1 as result.

Use

result = mvncdf([-Inf -Inf],[Inf Inf],[50 50],[sqrt(6),0;0, sqrt(6)])

result = 1

instead.

[Matthieu]

Thanks. I was beginning to wonder about the algorithm integral2 uses.

[Torsten]

But imagine you search for a small peak in the large infinite plane. I wonder which algorithm should help here. It shows integrating special functions with general integrators often causes trouble and special tools have to be supplied that are adequate for the special function in question.

The method used in mvncdf is especially suited for the multivariate normal pdf.

Answer 2

[Askic V]

@Torsten already answered, please keep in mind that integral2 still performs numerical integration.

have a look at this:

clear
clc
close all
sigma = sqrt(6) ;
Theta = 50 ;                                       
g = @(x,y) (normpdf(x,Theta,sigma) .* normpdf(y,Theta,sigma));
Q2 = integral2(g,-inf,inf,-inf,inf, 'Method','iterated','AbsTol',0)

Q2 = 1.0000

[Matthieu]

Yes, but I was surprised the loss was brutal around 34 from 1.0000 to 1e-10. Thanks for your addition to the answer.

Answer 3

[John D'Errico]

Why should you be surprised?

Integral and integral2 are tools that look at your function as a complete black box. A black box in this context means the code passes some numbers in, and then looks at the output. They cannot look at your code and "understand" it, parsing it, and then somehow know this is a Normal distribution involved, and to look near the mode of that normal? The code would then need to go back, and understand how the distribution is parameterized, and from that, infer where to look. The tool does not know the function has a very localised peak at one point in the domain, and where that may be.

A numerical integration tool looks at the function values it sees. If they are all zero to within a very high degree of precision EVERY where it has looked, what could you expect? A normal PDF has the property that in a moderately small region, so within plus or minus 6 or 10 sigma, it is ZERO to a fairly high tolerance.

normpdf(0)

ans = 0.3989

normpdf(6)

ans = 6.0759e-09

normpdf(30)

ans = 1.4736e-196

What did you do? You passed the integration tool limits of -inf and +inf. Is there a wide range between -inf and +inf? What are the odds the integration tool should somehow know to look in a VERY relatively narrow band for something to happen?