solve() returning empty, does a solution exist?

Hey all. For a project I'm doing, I need to generate a random number from a custom probability density function. To do that, I am generating a random number using rand, and finding its corresponding y-value on the related cumulative distribution function. However, solve() and vpasolve() only seem to return answers for some of these random values. Looking at the graph of this function, however, it would seem that these values exist. The relevant code, as well as the cumulative distribution function Fx, are below. Is it possible that my cumulative distribution function is too complicated, or that a solution doesn't exist after all? Any help would be greatly appreciated! Also, I am new to Matlab, so apologies if my code contains rookie mistakes!

clear all;
%Define the mean and standard deviation of the target normal distribution
%curve.
mu = 130;
stdev = 5/2;

%Define the function whose input is being tested (fy), as well as the target
%normal distribution PDF (normdist).
syms x p t;
fy(x) =  (-0.003867*x^3 + 0.01413*x^2 + 3.742*x + 101.9);
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) == 0, x);
localextrema_yvals = fy(localextrema_xvals);
localmin = [localextrema_xvals(1) localextrema_yvals(1)];
localmax = [localextrema_xvals(2) localextrema_yvals(2)];
localmin_xints = vpasolve(fy == localmin(2), x, [-50 50]);
localmax_xints = vpasolve(fy == localmax(2), x, [-50, 50]);
interval1ub = localmax_xints(1);
interval3lb = localmin_xints;

%Set up function (piecewise if necessary)
interval1 = int(fxtemp, t, -Inf, interval1ub);
interval2 = int(fxtemp, t, localmin(1), localmax(1));

Fx(x) = piecewise(x < interval1ub, int(fxtemp, t, -Inf, x), ...
        (x >= localmin(1)) & (x <= localmax(1)), interval1 + int(fxtemp, t, localmin(1), x),...
        x > interval3lb, interval1 + interval2 + int(fxtemp, t, interval3lb, x));

%Example of solve that cannot be done
xval = solve(Fx == 0.0284, x);

% %Generate 10000 random x values via inverse CDF method, not working!
% randvals = rand(10000, 1);
% xvals = zeros(10000, 1, "sym");
% for i = 1:10000
%     xvals(i) = solve(Fx == randvals(i), x);
% end
% yvals = fy(xvals);
% yvals = double(yvals);

Answer 1

[John D'Errico]

Just because solve returns empty does not mean a solution does not exist. In fact, most relations you can write down have absolutely no analytical solution.

When solve does fail to find a nice pleasing analytical solution however, it then tries to pass the problem to vpasolve. And USUALLY, vpasolve can find a solution, if a solution exists for a nice, well-behaved real valued function. HOWEVER...

Again, however, vpasolve need not always succeed. It is just a numerical solver, and that means it uses a starting value. The default starting value as I recall, is a random number between 0 and 1. And sometimes I have seen a starting value in that interval will diverge. Let me see if I can come up with a simple example.

syms x real
f(x) = exp(x) + 5 - 70*exp(-(x-4)^2*750)

f(x) = 

fplot(f,[-2,4.5])
grid on

That is a pretty simple function, one that has no analytical solution.

solve(f)

Warning: Unable to find explicit solution. For options, see help.

 
ans =
 
Empty sym: 0-by-1
 

vpasolve(f)

 
ans =
 
Empty sym: 0-by-1
 

As you can see, vpasolve also fails. I put a double solution in there, very near x==4. If I do start the solver in the vicinty of the known solution, it does succeed. The problem is, that solution is a difficult one to resolve, because it comes from a virtual spike in the function shape.

vpasolve(f,4)

ans = 

4.0140472717157835578354696409968

In your case of course, your function is purely monotomic. And that means vpasolve SHOULD pretty much always suceed. So, looking more carefully at your problem, what happens?

mu = 130;
stdev = 5/2;
%Define the function whose input is being tested (fy), as well as the target
%normal distribution PDF (normdist).
syms x p t;
fy(x) =  (-0.003867*x^3 + 0.01413*x^2 + 3.742*x + 101.9);
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) == 0, x);
localextrema_yvals = fy(localextrema_xvals);
localmin = [localextrema_xvals(1) localextrema_yvals(1)];
localmax = [localextrema_xvals(2) localextrema_yvals(2)];
localmin_xints = vpasolve(fy == localmin(2), x, [-50 50]);
localmax_xints = vpasolve(fy == localmax(2), x, [-50, 50]);
interval1ub = localmax_xints(1);
interval3lb = localmin_xints;
%Set up function (piecewise if necessary)
interval1 = int(fxtemp, t, -Inf, interval1ub);
interval2 = int(fxtemp, t, localmin(1), localmax(1));
Fx(x) = piecewise(x < interval1ub, int(fxtemp, t, -Inf, x), ...
    (x >= localmin(1)) & (x <= localmax(1)), interval1 + int(fxtemp, t, localmin(1), x),...
    x > interval3lb, interval1 + interval2 + int(fxtemp, t, interval3lb, x));

I do admit that Fx is virtually constant on the interval [0,1], which is again where the starting values for vpasolve will be. And that can potentially cause a problem, but not here, as long as you use vpasolve properly.

Fx(0) - 0.0284

ans = 

Fx(1) - 0.0284

ans = 

However, your problem is a subtly different one. Rather than trying to use vpasolve to find an exact solution to that equality statement, you want to SUBTRACT OFF 0.0284.

vpasolve(Fx(x) == 0.0284)

 
ans =
 
Empty sym: 0-by-1
 

vpasolve(Fx(x) - 0.0284) 

ans = 

6.3490049686606501343309322051234

Indeed, that worked nicely. Why did it work, when the first call to vpasolve failed? Look at what MATLAB sees for these two expressions:

x0 = 6;
Fx(x0) == 0.0284

ans = 

Fx(x0) - 0.0284

ans = 

Warning: Graphics acceleration hardware is unavailable. Graphics quality and performance might be diminished. See MATLAB System Requirements.

Do you see the equality statement will alsys return false for any number in the vicinity of the root? At the same time, the difference of Fx and a number returns a floating point number. Now vpasolve is a happy camper.

[Aidan]

Thanks for your help! Your recommendation worked, but I'm still a bit curious as to why it did! According to the documentation for vpasolve, if a symbolic expression is the input rather than a symbolic equation, the solver assumes that the right side of the equation is 0, and solves eqn == 0. However, these two code segments produce different outputs:

vpasolve(Fx(x) - 0.0284)

ans = 6.3490049686606501343309322051234

vpasolve(Fx(x) - 0.0284 == 0)

ans =

Empty sym: 0-by-1

Why is that? Does the first one have a tolerance that the second one does not have?

[Walter Roberson]

vpasolve(Fx(x) - 0.0284 == 0) tries to solve for exact equality.

vpasolve(Fx(x) - 0.0284) looks for a numeric solution and is willing to accept a numeric approximation of zero.

Chances are that up to the maximum number of digits, that Fx(x) - 0.0284 produces an expression that is slightly below zero for candidate x, and is slightly greater than zero for the numerically adjacent x. vpasolve(Fx(x) - 0.0284) says that situation is good enough. vpasolve(Fx(x) - 0.0284 == 0) compares to exact 0 and says that the comparisons fail on both sides, so says there is no solution.

[John D'Errico]

Again, look at the result returned by Fx == 0.0284.

The result is either true or false. And so vpasolve will only ever be happy if the result is EXACTLY true.

The same problem arises when you write

Fx(x) - 0.0284 == 0

Again, this is either true or false. Remember this case?

x0 = 6;
Fx(6) == 0.0284
ans =
0.0084279305055892282445002708293545 == 0.0284

Look carefully at ans. What do you see? It is a statement of potential equality. Is the left hand side IDENTICALLY equal to the right hand side? Of course not. In fact, the two sides will be non-identical, and so the expression is false for essentially all but one real value of x.

But vpasolve is a rootfinder. It wants to see a function that varies, and is smoothly differentiable. It is different for this:

Fx(x0) - 0.0284
ans =
-0.019972069494410771755499729170645

When you write the latter form, the result is a NUMBER. And when you do something like

vpasolve(fun)

then vpasolve looks for a zero of fun, whatever fun happens to be. This is not a question of a tolerance, but a question of what vpasolve sees, what you pass to vpasolve.

All of this is different when you use solve, since solve can look at an equality statement and understand it. That is, if we do this:

syms y
solve(y - 3 == 0)
ans =
3

solve can handle the equailty. But when solve is unable to find a solution, then it calls vpasolve. And vpasolve does not like the equality.