why PSI values are diffrent when nlmefit using ode45 solver compared to when it uses analytical solution?
2 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Sepi Sharif
el 8 de Sept. de 2021
Comentada: Sepi Sharif
el 20 de Sept. de 2021
nlmefit estimated parameters differently when it uses analytical solution compared to when it uses ode45 solver for the same dataset. PSI values are completely underestimated.
[beta,PSI,stats,b] = nlmefit(timeNew,DVnew,IDnew,dosePerID,model,logparameter0,...
'REParamsSelect',[1 2 3 ],'ErrorModel','Constant','ApproximationType','FOCE','CovPattern',P,'RefineBeta0','on',...
'ParamTransform',repelem(1,3),'OptimFun','fminunc');
function concentration=Getsolutionode(parameter,t,doseNew)
%% initial values
V=parameter(:,1);
ka=parameter(:,2);
CL=parameter(:,3);
predictor(:,1)=t;
minTime=min(t);
maxTime=max(t);
y0=[doseNew 0];
%% ode45
options=odeset('RelTol',10-9,'AbsTol',10-9);
solode=ode45(@odefunction,[minTime maxTime], [y0],options,V,ka,CL);
rspl=deval(solode,t);
concentration=rspl(2,:)'./V;
end
function dxdy = odefunction(~,y,V,ka,CL)
kel=CL/V;
dxdy= [-ka*y(1)
ka*y(1)-kel*y(2)];
end
11 comentarios
Respuesta aceptada
Star Strider
el 8 de Sept. de 2021
From the description, it appears that the code is attempting to fit data to the solution of a differential equation, or a system of differential equations.
To do this with ode45 (or more appropriately ode15s for kinetic equations), see Coefficient estimation for a system of coupled ODEs for a useful illustration of a similar problem.
.
7 comentarios
Star Strider
el 8 de Sept. de 2021
As always, my pleasure!
I see that Arthur Goldsipe has responded, so I will leave you in his capable hands.
.
Más respuestas (2)
John D'Errico
el 8 de Sept. de 2021
Editada: John D'Errico
el 8 de Sept. de 2021
There can sometimes be problems when you try to apply a nonlinear estimation scheme on top of an ODE solver. The nonlinear fitting tool presumes that the objective function is a differentiable function of the parameters. But is it? Remember that an ODE solver is a tool designed to work with a tolerance. Change the parameters slightly, and you will get a subtly different result, but only to within a tolerance. Effectively, your objective function is not a differentiable function of the parameters, because there is an arbitrary and unpredictable fuzz on top of your objective.
How is this different when you use an analytical solution to the problem? Now the objective will generally be a smooth, differentiable function of the parameters.
Effectively, this CAN be a problem anytime you layer one numerical solver on top of another. One thing to try is to crank down on the tolerance of the internal solver. But will that always be sufficient? Possibly not, since if the problem is a stiff one, then the solution can vary significantly as you cary the parameters.
Another trick may be to use a stiff solver. ODE45 is a tool that is not designed to solve stiff problems. So you may gain by simply moving to a stiff solver. That is, solvers where the name ends in s, like ODE15s, or ODE23s.
Can I know that this is the reason? Not without a concerted effort to dig into exactly what you did. Is it possible that your ODE may be a little stiff? Possibly. That would certainly exaggerate any issues in the solve. And of course, there may be a bug in the code you wrote. :) Surely not! That would never happen. ;)
1 comentario
Arthur Goldsipe
el 8 de Sept. de 2021
The reason the results from nlmefit differ significantly is because your ODE solution is not equivalent to your analyitical solution. I think you need to update file Getsolutionode.m to set minTime to 0 rather than min(t). Otherwise, your initial dose condition gets applied at min(t) instead of at time 0, and the ODE solution essentially has a time-delay when compared to your analytical solution.
In addition, as I mention in my comments on the question, I recommend changing 10-9 to 1e-9 when setting the tolerances. Otherwise, you will be using extremely loose tolerances when solving the ODEs.
Finally, I see from your example code that you are looking at a standard problem from pharmacokinetics (the "warfarin" dataset). This is exactly the sort of problem that SimBiology is designed to solve. You can build your pharmacokinetic model and run nlmefit and other analyses using SimBiology's graphical user interface or using command line functions and scripts. Here are some advantages I see to using SimBiology to solve this sort of problem:
- You can build your models graphically.
- You can select many standard PK models (including this one) from a built-in library.
- You can represent your model in a more natural way, using compartments, species, and reactions.
- You don't have to worry about the "bookkeeping" required to write the underlying differential equations.
- You can use built-in programs/tasks in the graphical user interface to perform mixed effects modeling.
- Your simulations almost always run faster. SimBiology can solve linear ODE models like this one analytically, using highly optimzed C++ code. And models that can't be solved anlytically can be solved much faster by taking advatange of its acceleration feature.
If you have a SimBiology license, I can show you how I would solve this problem in SimBiology. Just let me know whether you prefer to use the graphical user interface, the command line, or both.
0 comentarios
Ver también
Categorías
Más información sobre Scan Parameter Ranges en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!