Weird results of Calculating Integration of Normal/Gaussian Distribution

5 次查看(过去 30 天)
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.

采纳的回答

Star Strider
Star Strider 2018-8-3
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 个评论
David Goodmanson
David Goodmanson 2018-8-3
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 2018-8-4
Thank you very much, Walter Roberson and David Goodmanson! Your suggestion of using the 'Waypoints' option do solve my problem.

请先登录,再进行评论。

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Calculus 的更多信息

Community Treasure Hunt

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

Start Hunting!

Translated by