solving equation with if statements

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 :

Answer 1

[Walter Roberson]

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

[antoine]

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]

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]

You only have one equation, you can only solve for one variable.

Or is your y a vector?

[Walter Roberson]

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]

thank you+++. Really appreciated! I will try and get back to you asap.

[antoine]

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]

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]

Please attach your values.

[antoine]

see attached

thank you!

[Walter Roberson]

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

[antoine]

thank you Walter!

I will give it a try but sounds promising.

[antoine]

Really feel bad but I can't make it work... can you share the full code for my input data? sorry but my lack of matlab experience is strickingly apparent.

thank you ++

[antoine]

Hi Walter. Can you let me know if it works on your machine ? Would you mind sharing the full code for my input? thank you so much!

[Walter Roberson]

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.

[antoine]

mmhhh.. that does help but alpha cannot be negative (it is a duration)

can we inforce alpha as positive? thank you!

[Walter Roberson]

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.

[antoine]

bummer... I will try to figure this out but this is definitely VERY helpful. I will keep you informed

Thank you so much for your patience and support. I learnt a lot

[Walter Roberson]

All of which to say that the model lacks explanatory power for the data.

[antoine]

If I run the code above once, I got the following

If 266 is alpha, this is really what I would expect... I am probably missing something (I used alpha guess=150)

best_parameters =

    0.0065    0.0060  266.0000

[Walter Roberson]

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.

[antoine]

ahaha.. I get it. thank you!

Answer 2

[antoine]

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

[Walter Roberson]

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.

Answer 3

[antoine]

good not know...Not sure I have the skills to do that.

I really appreciated your input and advices.

[antoine]

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