Numerical problems when plotting with fplot
5 次查看(过去 30 天)
显示 更早的评论
The following are two different methods to plot an exponentially modified Gaussian distribution, https://en.wikipedia.org/wiki/Exponentially_modified_Gaussian_distribution
with parameters mu, sigma, lambda.
The first one is a point-by-point calculation of the integrals based on a for cycle; the second one is direct.
WIth parameters mu=0, sigma=1, lambda=1 all works fine, and the two plots are identical. Instead, if one sets mu = 0; sigma = 2; lambda = 3.5 we have an oscillation in the second plot. If we set mu = 0; sigma = 3; lambda = 3, the second plot is completely crazy.
I cannot understand why, since we have not pointS of singularity, that can affect the numerical convergence. Note, however, that the integral between -Inf and +Inf is in any case correct, i.e. =1, since this is a probability density function.
Here are the codes:
% FIRST METHOD
syms t
mu = 0;
sigma = 2;
lambda = 3.5;
x = -5:0.1:5;
F_lambda = zeros(1,length(x));
exp_t = @(t) exp(-t.^2);
for i=1:numel(x)
xmax = +Inf;
xmin = ((mu+lambda.*sigma.^2-x(i))./(sqrt(2).*sigma));
erfc = (2./sqrt(pi)).*int(exp_t,t,xmin,xmax);
f_lambda = (lambda./2).*exp((lambda./2).*(2.*mu+lambda.*sigma.^2-2.*x(i))).*erfc;
F_lambda(i) = f_lambda;
end
plot(x,F_lambda)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1398269/image.png)
% SECOND METHOD
syms x t
mu = 0;
sigma = 2;
lambda = 3.5;
exp_t = @(t) exp(-t.^2);
xmax = +Inf;
xmin = ((mu+lambda.*sigma.^2-x)./(sqrt(2).*sigma));
erfc = (2./sqrt(pi)).*int(exp_t,t,xmin,xmax);
f_lambda = (lambda./2).*exp((lambda./2).*(2.*mu+lambda.*sigma.^2-2.*x)).*erfc;
double(int(@(x) f_lambda,x,-Inf,+Inf))
fplot(f_lambda)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1398274/image.png)
0 个评论
采纳的回答
Paul
2023-5-30
Hi FRANCESCO,
I can't explain the inner workings of fplot, but I too have had occasion where it comes up with odd results.
Here's the second method, defining all constants as sym objects.
syms x t
mu = sym(0);
sigma = sym(2);
lambda = sym(3.5);
Define exp_t as a sym expression, didn't see a need to make it an anonymous function
exp_t = exp(-t.^2);
xmax = +Inf;
xmin = ((mu+lambda.*sigma.^2-x)./(sqrt(2).*sigma));
Define erfc expression. Matlab is smart enough to identify the integral and express erfc in term of erf
erfc = (2./sqrt(sym(pi))).*int(exp_t,t,xmin,xmax)
Here, I'll define f_lambda as a symfun.
f_lambda(x) = (lambda./2).*exp((lambda./2).*(2.*mu+lambda.*sigma.^2-2.*x)).*erfc
%double(int(@(x) f_lambda,x,-Inf,+Inf))
figure
fplot(f_lambda,'MeshDensity',200)
But the Symbolic Math Toolbox has its own erfc function (and there is a numeric version as well erfc), so maybe things will work better if we use that function.
clear erfc
f_lambda(x) = (lambda./2).*exp((lambda./2).*(2.*mu+lambda.*sigma.^2-2.*x))*erfc(xmin)
Same problem with fplot
figure
fplot(f_lambda,'MeshDensity',200)
But we can plot f_lambda at specified values of x
xval = -5:0.1:5;
figure
plot(xval,double(f_lambda(xval)))
2 个评论
Walter Roberson
2023-5-30
@Paul I notice the erfc comes out as erf(expression)+1 . I would have expected 1 - erf(some expression) ??
Paul
2023-5-30
I noticed that too and it threw me for a loop. But when I used the SMT version of erfc, it came up with an expression containing -(efrc() - 2), which is the same as
(2 - erfc()) = 1 + 1 - erfc() = 1 + erf()
so I didn't pursue any further. Probably has to do with erf being an odd function so that:
erfc(a-x) = 1 - erf(a-x) = 1 + erf(x-a)
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Applications 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!