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

2 visualizzazioni (ultimi 30 giorni)
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.

Risposta accettata

Torsten
Torsten il 6 Mar 2023
Modificato: Torsten il 6 Mar 2023
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.
  2 Commenti
Matthieu
Matthieu il 6 Mar 2023
Thanks. I was beginning to wonder about the algorithm integral2 uses.
Torsten
Torsten il 6 Mar 2023
Modificato: Torsten il 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.

Accedi per commentare.

Più risposte (2)

Askic V
Askic V il 6 Mar 2023
Modificato: Askic V il 6 Mar 2023
@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
  1 Commento
Matthieu
Matthieu il 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.

Accedi per commentare.


John D'Errico
John D'Errico il 6 Mar 2023
Modificato: John D'Errico il 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.
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?

Community Treasure Hunt

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

Start Hunting!

Translated by