You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
solving equation with if statements
2 views (last 30 days)
Show older comments
Hello
I am trying to solve the following equation given input time and values (y) vectors.
I tried the following
syms theta kappa alpha
if time < alpha
eqn = exp(-kappa*time)==y;
else
eqn = exp(-kappa)*exp(-theta*(time-alpha))==y;
end
vars = [theta kappa alpha];
S=solve(eqn);
but it raises an error...
Conversion to logical from sym is not possible.
I tried to convert alpha into double (i.e. double(alpha) ) in the if statement but it did not work...
Any suggestions?
the data look like :
Accepted Answer
Walter Roberson
on 17 Jan 2020
syms theta kappa alpha time y real
eqn = y == piecewise(time < alpha, exp(-kappa*time), exp(-kappa)*exp(-theta*(time-alpha)));
vars = [theta kappa alpha];
S = solve(eqn, kappa, 'returnconditions', true); %solve for which variable ??
21 Comments
antoine
on 17 Jan 2020
thank you so much Walter! I want to solve all vars
Here 's what I did (mat files are vectors of proportion of survival overtime)
%%%equations
load('values_bupr1.mat')
load('days_bupr1.mat')
y=values_bupr_nodup;
time=days_bupr_nodup;
syms theta kappa alpha time y real
eqn = y == piecewise(time < alpha, exp(-kappa*time), exp(-kappa)*exp(-theta*(time-alpha)));
vars = [theta kappa alpha];S
S = solve(eqn, vars, 'returnconditions', true);
and I got empty output (I expect a value for each var)
S =
struct with fields:
theta: [2×1 sym]
kappa: [2×1 sym]
alpha: [2×1 sym]
parameters: [1×4 sym]
conditions: [2×1 sym]
Does that make sense?
antoine
on 17 Jan 2020
to complete my previous comment. The output are
>> S.kappa
ans =
-log(y)/time
z
>> S.alpha
ans =
u
z1
>> S.theta
ans =
x
log(y*exp(z))/(z1 - time)
but I don't know how to get a value for each of my var...I am probably missing an obvious step
Walter Roberson
on 18 Jan 2020
Edited: Walter Roberson
on 18 Jan 2020
You only have one equation, you can only solve for one variable.
Or is your y a vector?
Walter Roberson
on 18 Jan 2020
syms theta kappa alpha Time real
eqn = piecewise(Time < alpha, exp(-kappa*Time), exp(-kappa)*exp(-theta*(Time-alpha)));
residue = sum((subs(eqn, Time, time) - y).^2);
residueF = matlabFunction(residue, 'vars', {[theta kappa alpha]}, 'file', 'residue.m'); %must write to file!
Now residueF is a function handle to minimimize, such as
kappa_guess = mean(-log(y)./time);
theta_guess = randn;
alpha_guess = randn;
p0 = [theta_guess, kappa_guess, alpha_guess];
[best_parameters, best_residue] = fminsearch(residueF, p0);
It would be common that any particular initial guess will not get you to the best fit; you will probably need to try several values.
You must use the 'file' option with matlabFunction in this situation: because the value of alpha is not known, the piecewise() cannot be resolved even though the times are known, so the expression will have have a piecewise() expression in it, and matlabFunction can only handle piecewise() when writing to a file.
antoine
on 18 Jan 2020
My progress so far:
the code above is not working if I try to solve all 3 variables... Since I expect alpha to be between 60 and 300 , I tried to simplify the code and assign an alpha like:
alpha=90; % my first guess
syms theta kappa Time real
eqn = piecewise(Time < alpha, exp(-kappa*Time), exp(-kappa)*exp(-theta*(Time-alpha)));
residue = sum((subs(eqn, Time, time) - y).^2);
residueF = matlabFunction(residue, 'vars', {[theta kappa ]}, 'file', 'residue.m'); %must write to file!
kappa_guess = mean(-log(y)./time);
theta_guess = randn;
p0 = [theta_guess, kappa_guess];
[best_parameters, best_residue] = fminsearch(residueF, p0);
Here your code is working ONLY and I got the following theta and kappa
kappa_guess = 0.0023
theta_guess = 0.5377
Any suggestion? any ways alpha can also be an output?
thank you!
(I can share the data if you want)
antoine
on 18 Jan 2020
using strictly the code you suggested;
load('values_bupr1.mat')
load('days_bupr1.mat')
y=values_bupr_nodup;
time=days_bupr_nodup;
syms theta kappa alpha Time real
eqn = piecewise(Time < alpha, exp(-kappa*Time), exp(-kappa)*exp(-theta*(Time-alpha)));
residue = sum((subs(eqn, Time, time) - y).^2);
residueF = matlabFunction(residue, 'vars', {[theta kappa alpha]}, 'file', 'residue.m'); %must write to file!
Now residueF is a function handle to minimimize, such as
kappa_guess = mean(-log(y)./time);
theta_guess = randn;
alpha_guess = randn;
p0 = [theta_guess, kappa_guess, alpha_guess];
[best_parameters, best_residue] = fminsearch(residueF, p0);
I have the following error
Error using symengine
Unable to evaluate to Boolean.
Error in sym/mupadmexnout (line 1057)
out = mupadmex(fcn,args{:});
Error in sym/matlabFunction>optimize (line 468)
[tvalues,f,tnames] = mupadmexnout('symobj::optimizeWithIntermediates',f{:});
Error in sym/matlabFunction>writeMATLAB (line 443)
[f,tvalues,tnames] = optimize(f,optim);
Error in sym/matlabFunction (line 183)
g = writeMATLAB(funs,file,varnames,outputs,body, opts.Optimize, opts.Sparse, opts.Comments);
Error in SOLVER_3VARS (line 14)
residueF = matlabFunction(residue, 'vars', {[theta kappa alpha]}, 'file', 'residue.m'); %must write to file!
Walter Roberson
on 18 Jan 2020
Please attach your values.
Walter Roberson
on 19 Jan 2020
It looks like it is probably easiest to construct
residue = @(tka) sum((project_y(tka, time) - y).^2);
function projected = project_y(tka, time)
theta = tka(1); kappa = tka(2); alpha = tka(3);
projected = exp(-kappa.*time);
mask = time >= alpha;
projected(mask) = exp(-kappa).*exp(-theta.*(time(mask)-alpha));
end
Walter Roberson
on 20 Jan 2020
load('values_bupr1.mat')
load('days_bupr1.mat')
y = values_bupr_nodup;
time = days_bupr_nodup;
residueF = @(tka) sum((project_y(tka, time) - y).^2);
kappa_guess = mean(-log(y)./time);
theta_guess = randn;
alpha_guess = 150;
p0 = [theta_guess, kappa_guess, alpha_guess];
opts = optimset('fminsearch');
opts = optimset(opts, 'MaxIter', 2000, 'MaxFunEvals', 10000, 'FunValCheck', 'on');
[best_parameters, best_residue] = fminsearch(residueF, p0, opts)
function projected = project_y(tka, time)
theta = tka(1); kappa = tka(2); alpha = tka(3);
projected = exp(-kappa.*time);
mask = time >= alpha;
projected(mask) = exp(-kappa).*exp(-theta.*(time(mask)-alpha));
end
My tests so far suggest that the best fit is near
theta = 0.000844900288029776
kappa = -556.133604457634
alpha = -659149.32154226
That is, with alpha so negative that time < alpha is false for all entries, and the exp(-theta*(time-alpha)) term heads towards 0.
Walter Roberson
on 20 Jan 2020
When you restrict to positive alpha, then fitting pretty much cannot tell the difference between large-ish negative theta and large-ish positive theta. The only consistency you get is that kappa narrows down to close to 0.00292734523205026836 and otherwise it becomes difficult to tell apart locations with theta anywhere above 788 or so, and the larger you make alpha the more decimal places you can grind out on essentially the same fit. At the moment, the best fit I have is at
theta = 5087.59871387777093
kappa = 0.00292734523205026836
alpha = 803477.437304047635
but the confidence intervals for theta and alpha are terrible.
Walter Roberson
on 20 Jan 2020
All of which to say that the model lacks explanatory power for the data.
Walter Roberson
on 20 Jan 2020
Yup, but that gives residue in the range of 62, whereas 5087.59871387777093, 0.00292734523205026836, 803477.437304047635 gives you residue in the range of 39.02 which is significantly better. Given that model, alpha should be large, over 100,000.
More Answers (2)
antoine
on 18 Jan 2020
time and y are vectors (see plot above for data input). all suggestions are welcome.
I can share the data if it helps
thank you ++
1 Comment
Walter Roberson
on 18 Jan 2020
Unless time and y are vectors of length 3 exactly, it seems unlikely that there would be any solutions.
It looks to me as if what you have should be a modeling task to estimate parameters rather than a set of simultaneous equations.
antoine
on 18 Jan 2020
good not know...Not sure I have the skills to do that.
I really appreciated your input and advices.
1 Comment
antoine
on 22 Feb 2020
dear Walter
Just realized I forgot to thank you for your help
Problem is fixed! I got my equations working.
(I am know stuck with other timer/loop function : see 'time-dependant iteration through a loop?' if you have time...)
See Also
Categories
Find more on Calculus in Help Center and File Exchange
Tags
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)