이 질문을 팔로우합니다.
- 팔로우하는 게시물 피드에서 업데이트를 확인할 수 있습니다.
- 정보 수신 기본 설정에 따라 이메일을 받을 수 있습니다.
solving equation with if statements
조회 수: 4 (최근 30일)
이전 댓글 표시
antoine
2020년 1월 17일
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 :
채택된 답변
Walter Roberson
2020년 1월 17일
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
antoine
2020년 1월 17일
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
2020년 1월 17일
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
2020년 1월 18일
편집: Walter Roberson
2020년 1월 18일
You only have one equation, you can only solve for one variable.
Or is your y a vector?
Walter Roberson
2020년 1월 18일
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
2020년 1월 18일
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
2020년 1월 18일
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
2020년 1월 19일
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
2020년 1월 19일
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
2020년 1월 20일
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
2020년 1월 20일
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
2020년 1월 20일
mmhhh.. that does help but alpha cannot be negative (it is a duration)
can we inforce alpha as positive? thank you!
Walter Roberson
2020년 1월 20일
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
2020년 1월 20일
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
2020년 1월 20일
All of which to say that the model lacks explanatory power for the data.
antoine
2020년 1월 20일
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
2020년 1월 20일
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.
추가 답변 (2개)
antoine
2020년 1월 18일
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
Walter Roberson
2020년 1월 18일
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
2020년 1월 18일
good not know...Not sure I have the skills to do that.
I really appreciated your input and advices.
댓글 수: 1
antoine
2020년 2월 22일
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...)
참고 항목
카테고리
Help Center 및 File Exchange에서 Calculus에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!오류 발생
페이지가 변경되었기 때문에 동작을 완료할 수 없습니다. 업데이트된 상태를 보려면 페이지를 다시 불러오십시오.
웹사이트 선택
번역된 콘텐츠를 보고 지역별 이벤트와 혜택을 살펴보려면 웹사이트를 선택하십시오. 현재 계신 지역에 따라 다음 웹사이트를 권장합니다:
또한 다음 목록에서 웹사이트를 선택하실 수도 있습니다.
사이트 성능 최적화 방법
최고의 사이트 성능을 위해 중국 사이트(중국어 또는 영어)를 선택하십시오. 현재 계신 지역에서는 다른 국가의 MathWorks 사이트 방문이 최적화되지 않았습니다.
미주
- América Latina (Español)
- Canada (English)
- United States (English)
유럽
- 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)
아시아 태평양
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)