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

3 vues (au cours des 30 derniers jours)
Matthieu le 6 Mar 2023
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.

Réponse acceptée

Torsten le 6 Mar 2023
Modifié(e) : Torsten le 6 Mar 2023
It's not possible for integral2 to capture the small peak around 50 to return 1 as result.
result = mvncdf([-Inf -Inf],[Inf Inf],[50 50],[sqrt(6),0;0, sqrt(6)])
result = 1
  2 commentaires
Matthieu le 6 Mar 2023
Thanks. I was beginning to wonder about the algorithm integral2 uses.
Torsten le 6 Mar 2023
Modifié(e) : Torsten le 6 Mar 2023
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.

Connectez-vous pour commenter.

Plus de réponses (2)

Askic V
Askic V le 6 Mar 2023
Modifié(e) : Askic V le 6 Mar 2023
@Torsten already answered, please keep in mind that integral2 still performs numerical integration.
have a look at this:
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
  1 commentaire
Matthieu le 6 Mar 2023
Yes, but I was surprised the loss was brutal around 34 from 1.0000 to 1e-10. Thanks for your addition to the answer.

Connectez-vous pour commenter.

John D'Errico
John D'Errico le 6 Mar 2023
Modifié(e) : John D'Errico le 6 Mar 2023
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.
ans = 0.3989
ans = 6.0759e-09
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?


En savoir plus sur Characters and Strings dans Help Center et File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by