Numerical problems when plotting with fplot

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)

% 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)

Answer 1

[Paul]

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)

erfc = 

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

f_lambda(x) = 

%double(int(@(x) f_lambda,x,-Inf,+Inf))

Sometimes, fplot works better by increasing the MeshDensity. Didn't help here.

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)

f_lambda(x) = 

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)))

[Walter Roberson]

@Paul I notice the erfc comes out as erf(expression)+1 . I would have expected 1 - erf(some expression) ??

[Paul]

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)