An ode system solution as nlinfit function

Hi everyone,

I am trying to use one of the solutions of a differential equations system to fit a set of data using nlinfit and nlparci since I found an example online. In particular, I think I am doing the iteration: solve the system -> obtain the solution -> use this solution to try to fit the data -> get better parameters from nlinfit and nlparci -> use this parameter to solve the system

and so on.

I would like to know a couple of thing:

1. Does nlinfit solve the system every iteration, in order to obtain the function, or not?
2. If not, how can I use the optimized parameters obtained from nlparci to solve again the system and to go on with the iteration?

Thank you for any help.

Here is my code:

Gamma_0 = 1.2; 
gamma_0 = 0.3; 
gamma_c_0 = 5.1; 
c_0 = 6.5;
A_0 = 1; 
t0_0 = -0.6; 
y0_0 = 0; 

parguess = [Gamma_0, gamma_0, gamma_c_0, c_0, A_0, t0_0, y0_0];


options = statset();
options.MaxIter = 1000; 
[pars, resid, J] = nlinfit(tdata,ydata,@myfunc,parguess,options);

alpha = 0.05; 
pars_co = nlparci(pars, resid, 'jacobian', J, 'alpha', alpha);


Gamma = pars(1);
gamma = pars(2);
gamma_c = pars(3);
c = pars(4);
A = pars(5);
t0 = pars(6);
y0 = pars(7);



parst = transpose(pars);
for i=1:length(pars)
    errors(1,i) = max(abs(parst(i)-pars_co(i,1)),abs(pars(i)-pars_co(i,2)));
    errors(2,i) = abs(errors(1,i)/parst(i))*100;
end


tfit = transpose(linspace(tdata(1),tdata(end),1000));
fit = myfunc(pars, tfit);

figure();
plot(tfit,fit,'b',tdata,ydata,'r');
legend('fit','data','FontSize',12);
figure();
semilogy(tfit,fit,'b',tdata,ydata,'r');
legend('fit','data','FontSize',12);


%%%% 
function uscita = myfunc(pars, tdata)
    Gamma = pars(1);
    gamma = pars(2);
    gamma_c = pars(3);
    c = pars(4);
    A = pars(5);
    t0 = pars(6);
    y0 = pars(7);
    k = 760;    
        function dydt = ode(t,y)
        dydt = zeros(8,1);
        dydt(1)=-y(1)*gamma_c*2*(exp(-gamma_c*(t+t0)))+y(4)*gamma+y(8)*k; %p00
        dydt(2)=y(1)*gamma_c*(exp(-gamma_c*(t+t0)))-y(2)*c*gamma_c*(exp(-gamma_c*(t+t0)))+y(5)*gamma; %p10
        dydt(3)=y(1)*gamma_c*(exp(-gamma_c*(t+t0)))-y(3)*c*gamma_c*(exp(-gamma_c*(t+t0)))+y(6)*gamma; %p01
        dydt(4)=-y(4)*(gamma+Gamma+gamma_c+2*(exp(-gamma_c*(t+t0))))+(y(2)+y(3))*c*gamma_c*(exp(-gamma_c*(t+t0)))+y(7)*gamma+y(8)*Gamma; %p11
        dydt(5)=y(4)*gamma_c*(exp(-gamma_c*(t+t0)))-y(5)*(gamma+c*gamma_c*(exp(-gamma_c*(t+t0)))); %p21
        dydt(6)=y(4)*gamma_c*(exp(-gamma_c*(t+t0)))-y(6)*(gamma+c*gamma_c*(exp(-gamma_c*(t+t0)))); %p12
        dydt(7)=(y(5)+y(6))*c*gamma_c*(exp(-gamma_c*(t+t0)))-y(7)*gamma; %p22
        dtdy(8)=y(4)*Gamma-y(8)*(k+Gamma); %pcm
        end
    y00 = 1;
    y10 = 0;
    y01 = 0;
    y11 = 0;
    y21 = 0;
    y12 = 0;
    y22 = 0;
    ycm = 0;
    [time,y] = ode45(@ode,tdata,[y00,y10,y01,y11,y21,y12,y22,ycm]);
    uscita = A*y(:,4)+y0;
end

Answer 1

[Star Strider]

‘Does nlinfit solve the system every iteration, in order to obtain the function, or not?’

It solves it at every iteration. It updates the ‘pars’ parameters and continues iterating until it converges on an acceptable parameter set. The nlparci results will be the parameter confidence intervals for the converged parameter set.

[Andrea Ristori]

Thank you for your answer.

So my code is ok and it should work, is it?

Since I tried to compare the results of my fit with the one obtained by another guy and they do not match.

[Star Strider]

My pleasure.

‘So my code is ok and it should work, is it?’

I have no idea. I cannot run it without your data, and I do not know what system you are simulating.

‘Since I tried to compare the results of my fit with the one obtained by another guy and they do not match.’

I would be less concerned about their matching and more concerned with how well the model fits the data. If the model is coded and implemented correctly, the best fit likely has the most accurate parameter estimates.

Onr thing that concerns me is:

y0 = pars(7);

and later the result that is passed back to nlinfit is:

uscita = A*y(:,4)+y0;

with ‘y0’ not being re-defined (that I can see) in the interim.

However I have no idea what you are doing, so that may be correct.

.

[Andrea Ristori]

With "So my code is ok and it should work, is it?" I meant that the process is correct and it does the iteration I need to do.

Regarding y0, do you mean that is not re-defined in the sytem of equations? If it is so, yes I wanted only to intoduce a constant value to the fitting function (and an amplitute, A), to have an additional parameter to better fit the data.

Thank you again.

[Star Strider]

As I wrote previously, I have no idea.

I have no idea what you are doing. If using ‘y0’ there is correct, then use it.