Weird results of Calculating Integration of Normal/Gaussian Distribution

I wrote a simple function as same as the Normal/Gaussian Distribution with mean of 0 and sigma as the standard deviation.

The integration results with "integral" function seemed so whimsical in some cases with specific boundaries.

To my knowledge, the integration of Normal Distribution with boundaries of (-inf,inf) should be 1, while the integration with boundaries of (-n*sigma,n*sigma) should be close enough to 1, when n>=5 .

clear;close all;format long
sigma = 0.05; % the standard deviation
G1 = @(X) exp(-X.^2/(2*sigma^2) )/(sigma*sqrt(2*pi) ); % the normal/Gaussian distribution
clc
n = 18;                         % -n*sigma defines the lower bound
x1 = 1000;                      % x1 defines the upper bound 1
x2 = 5000;                      % x2 defines the upper bound 2
[integral(G1, -n*sigma, 0),1/2; % should be close to 1/2
integral(G1, 0, x1),1/2;        % should be close to 1/2
integral(G1, -n*sigma, x1),1;   % should be close to 1
integral(G1, x1, x2),0;         % should be close to 0
integral(G1, -n*sigma, x2),1]   % should be close to 1

However, the result I got was like below:

ans =

   0.500000000000000   0.500000000000000
   0.500000000000252   0.500000000000000
   0.000000000013511   1.000000000000000
                   0                   0
   1.000000000001456   1.000000000000000

The 3rd row of the results, which was expected to be close to 1, was calculated to be 0.000000000013511

Besides, if the value of x1 was changed as x1 = 988, the result would be 0.999999999999933, which was quite close to the correct answer of 1; However, if x1 = 989 was chosen, the result would immediately turned to 0.000000000099978, which made me surprise.

My question is, why the "integral" function is so sensitive to the boundaries that even a slightly changing would lead to larger difference between the results? And, how could the "weird result" be avoided or excluded?

Thank you.

Answer 1

[Star Strider]

See the documentation section on: Accuracy of Floating-Point Data (link).

You are seeing the result of floating-point approximation error. There are many discussions of it in MATLAB Answers, so I’ll not go into detail here.

[leon t]

Thank you Star Strider for your advice.

I would chew the link you provided over first.

By the way, is there anything I could do with the initial code to prevent the catastrophic floating-point accuracy problem?

I tried the "int" function and it worked perfectly, however, the "int" function would be insufficient in calculating the value of integration of a more complicated G(x) (in that case, it seems that I have to use the "integral" function for a numerical solution, which may be effected by the accuracy problem ).

Would you please send me some suggestions? I'd be very grateful.

Again, thank you for your help.

[Star Strider]

My pleasure.

Unfortunately, we’re all stuck with the constraints of floating-point calculations. (For even further enlightenment, see: Why is 0.3 - 0.2 - 0.1 (or similar) not equal to zero? (link).) I doubt that there is anything you (or anyone else) can do to optimise your code in that respect.

The int function is part of the Symbolic Math Toolbox. It has extended precision, so it will give a more precise result. However, if you want to use its results in the rest of your code, you will again encounter the floating-point approximation problem.

One option could be to test every value that is sufficiently close to zero and assign it as zero:

zeroforce = @(f) f.*(abs(f) > 1E-10);

x1 = 0.000000000013511;
x2 = pi;

y1 = zeroforce(x1)
y2 = zeroforce(x2)

y1 =
     0

y2 =
   3.141592653589793

That’s the best I can do.

[Walter Roberson]

You should consider using the 'waypoints' option of integral()

[David Goodmanson]

Hi leon,

I don't think this issue has as much to do with floating point accuracy as it has to do with the nature of the function. G1 is a very sharply peaked function. At X = =-1 it is has a value of 1.1e-86 and is dropping like a rock. Past X = +-2 you get G1 = 0.

When you set integration limits of 1000 or so, as far as 'integral' is concerned you are setting a scale for the problem. Will 'integral' stumble across the peak, or not? If it does you get basically 1, and if not you get basically 0.

For large X, integral(G1,0,X) = 1/2, since the lower limit is on top of the peak and the peak is found automatically. Here X could be 1000 but 'large' just means X>2.

integral(G1,-10,10) = 1 since 10 is comparable to the peak width and the peak will be found. integral(G1,-1000,1000) misses. And so forth.

For a well chosen waypoint, I would think that Walter's waypoint suggestion assures getting the right answer.

[leon t]

Thank you very much, Walter Roberson and David Goodmanson! Your suggestion of using the 'Waypoints' option do solve my problem.