Weird results of Calculating Integration of Normal/Gaussian Distribution

7 views (last 30 days)
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.

Accepted Answer

Star Strider
Star Strider on 3 Aug 2018
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.
  5 Comments
David Goodmanson
David Goodmanson on 3 Aug 2018
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
leon t on 4 Aug 2018
Thank you very much, Walter Roberson and David Goodmanson! Your suggestion of using the 'Waypoints' option do solve my problem.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!