Symbolic expression output changes when turned into function?

Hey all! I'm having a weird issue with my code where defining an expression as a function changes the output. For instance, I defined a symbolic function fx with the following code:

clear all;
%Define the mean and standard deviation of the target normal distribution
%curve.
mu = 7.5;
stdev = 2.5/2;
%Define the function whose input is being tested (fy), as well as the target
%normal distribution PDF (normdist).
syms x p t real
fy(x) =  (4.763*x^3 - 7.939*x^2 - 10.62*x + 20.78);
normdist(p) = (1/(sqrt(2*pi*(stdev)^2))*exp(-(p - mu)^2/(2*(stdev)^2)));

%Create PDF of x (fx)

fx(x) = (normdist(fy) * abs(diff(fy)));
%Temp fx function in terms of t for integrating
fxtemp(t) = fx(t);

%Create CDF of x (Fx)
%Determine integral bounds (if fy is not monotonic)
localextrema_xvals = vpasolve(diff(fy), x);
localextrema_yvals = fy(localextrema_xvals);
localmin = [localextrema_xvals(2) localextrema_yvals(2)];
localmax = [localextrema_xvals(1) localextrema_yvals(1)];
localmin_xints = vpasolve(fy - localmin(2), x, [-50 50]);
localmax_xints = vpasolve(fy - localmax(2), x, [-50, 50]);
interval1ub = localmin_xints;
interval3lb = localmax_xints;

Following that, I wished to define a function whose variable was the upper bound of an integral of fx. However, the function was not behaving the way I expected it to. After some testing, I realized that the behavior of the integral expression changes once it is defined as a function, as shown by this code:

int(fxtemp, t, localmax(1), 1)

ans = 

0.66012342427575267765952529754189

f2 = int(fxtemp, t, localmax(1), x);
subs(f2, x, 1)

ans = 

1.6601234242757526705599930293611

In this case, the function displays an answer that does not make sense. Strangely, this issue disappears when localmax(1) is replaced with its equivalent value:

double(int(fxtemp, t, -0.47003, 1))

ans = 0.6601

f1 = int(fxtemp, t, -0.47003, x);
double(subs(f1, x, 1))

ans = 0.6601

If anyone has any ideas as to what is going on, it would be greatly appreciated!

Answer 1

[Walter Roberson]

%Define the mean and standard deviation of the target normal distribution
%curve.
mu = 7.5;
stdev = 2.5/2;
%Define the function whose input is being tested (fy), as well as the target
%normal distribution PDF (normdist).
syms x p t real
fy(x) =  (4.763*x^3 - 7.939*x^2 - 10.62*x + 20.78);
normdist(p) = (1/(sqrt(2*pi*(stdev)^2))*exp(-(p - mu)^2/(2*(stdev)^2)));

%Create PDF of x (fx)

fx(x) = (normdist(fy) * abs(diff(fy)));
%Temp fx function in terms of t for integrating
fxtemp(t) = fx(t);

%Create CDF of x (Fx)
%Determine integral bounds (if fy is not monotonic)
localextrema_xvals = vpasolve(diff(fy), x);
localextrema_yvals = fy(localextrema_xvals);
localmin = [localextrema_xvals(2) localextrema_yvals(2)];
localmax = [localextrema_xvals(1) localextrema_yvals(1)];
localmin_xints = vpasolve(fy - localmin(2), x, [-50 50]);
localmax_xints = vpasolve(fy - localmax(2), x, [-50, 50]);
interval1ub = localmin_xints;
interval3lb = localmax_xints;

syms A B real
assumeAlso(A<=B)
II = int(fxtemp, t, A, B)

II = 

II1 = subs(II, B, 1)

II1 = 

subs(II1, A, localmax(1))

ans = 

vpa(ans)

ans = 

1.1601234242757526741097591634515

vpa(II1)

ans = 

subs(ans, A, localmax(1))

ans = 

1.6601234242757526705599930293611

Your localmax(1) is right beside the boundary, and that is throwing everything off.

Your localmax(1) is close enough to the boundary that the fact that you derived it with vpasolve() might be significant.

[Walter Roberson]

format long g
sympref('FloatingPointOutput', false)

ans = logical
   1

%Define the mean and standard deviation of the target normal distribution
%curve.
mu = 7.5;
stdev = 2.5/2;
%Define the function whose input is being tested (fy), as well as the target
%normal distribution PDF (normdist).
syms x p t real
fy(x) =  (4.763*x^3 - 7.939*x^2 - 10.62*x + 20.78);
normdist(p) = (1/(sqrt(2*pi*(stdev)^2))*exp(-(p - mu)^2/(2*(stdev)^2)));

%Create PDF of x (fx)

fx(x) = (normdist(fy) * abs(diff(fy)));
%Temp fx function in terms of t for integrating
fxtemp(t) = fx(t);

%Create CDF of x (Fx)
%Determine integral bounds (if fy is not monotonic)
localextrema_xvals = solve(diff(fy), x)

localextrema_xvals = 

localextrema_yvals = fy(localextrema_xvals);
localmin = [localextrema_xvals(2) localextrema_yvals(2)];
localmax = [localextrema_xvals(1) localextrema_yvals(1)];
localmin_xints = vpasolve(fy - localmin(2), x, [-50 50]);
localmax_xints = vpasolve(fy - localmax(2), x, [-50, 50]);
interval1ub = localmin_xints;
interval3lb = localmax_xints;

syms A B real
assumeAlso(A<=B)
II = int(fxtemp, t, A, B)

II = 

II1 = subs(II, B, 1)

II1 = 

●

subs(II1, A, localmax(1))

ans = 

vpa(ans)

ans = 

0.1601234242757526812092914316323

vpa(II1)

ans = 

subs(ans, A, localmax(1))

ans = 

vpa(ans)

ans = 

0.66012342427575267765952529754189

So localmax(1) is so close to the boundary that evaluting at full precision and taking vpa() of the result gets you about 0.16, but vpa() first and then evaluating and vpa() the result gets you about 0.66

[Aidan]

Sorry, I'm a little bit confused as to what you mean. Close to what boundary? Also, why does only localmax(1) cause the issue and not the value it stores? If it helps at all I'm fairly certain the function is off from the expression by a flat value of 1, since evaluating the function at localmax(1), which should be 0, just yields 1.

[Walter Roberson]

The indefinite integral over symbolic bounds involves several piecewise() tests. When 1 is substituted for the upper bound, there are several piecewise conditions remaining that test the lower bound of the integral.

The exact bounds determined by the piecewise() are either at or are immediately adjacent numerically to the value held in localmax(1) . If you use solve() to derive localmax(1) then you get something of the form c/d - sqrt(e)/f and if you compare that to the symbolic bounds of the piecewise() components, you will find that there is an exact match . The indefinitely precise form of localmax(1) is exactly the same as the indefinitely precise form of the piecewise() bounds.

That being the case, exactly when you take the vpa() of localmax(1) and exactly when you take the vpa() of the piecewise bounds matters, because the comparison matters to the very last decimal place.

Your localmax(1) does not hold the exact value -0.47003

xyz = sym('-0.47003070272839625365530652424196')

xyz = 

format short
double(xyz)

ans = -0.4700

format long g
double(xyz)

ans = 
        -0.470030702728396

You have confused the display of the value with the value stored. The value stored is in symbolic floating point form, and has at least 32 decimal places (if you probe more closely you will find that the value stored has over 70 decimal places.) You are asking for double() of it while under the influence of "format short". "format short" does not change the value stored for it, only how much of it is displayed. format long g displays more of the internal number (but at that, format long g misses between 1 and 4 bits worth of the actual number.)

[Aidan]

I see! I had actually tested out the numeric replacement initially by copying down the value associated with localmax(1) in the workspace. So essentially, changing localmax(1) to be even one decimal place less precise fixes the issue since the issue only arises because localmax(1) and the boundary of the piecewise function generated by Matlab to represent the integral are exactly the same? I have made a similar function before seemingly without issue, which is why this confused me. Was this by chance, or is there a more systematic way to avoid this issue? Will the issue always be so obvious (like in this case where there was just an offset of 1 from the expected value) or could the effects be more insidious? Thanks for your help!