lsqcurvefit Fitting kinetic model to estimate kinetic parameter.
Show older comments
I found this code in the community and used it for my own data(Thanks to the provider of the code), but no matter how I ran it, the data didn't fit well. How can I find the right initial value? Or what other function should I use to estimate my parameters?
These are the parameters I got from running.
Theta( 1) = 0.00221
Theta( 2) = 0.00427
Theta( 3) = 0.00042
Theta( 4) = 0.00816
Theta( 5) = 0.00076
Theta( 6) = 0.00000
Theta( 7) = 0.00000
Theta( 8) = 0.00000
Theta( 9) = 0.00002
Theta(10) = 2.94054
Theta(11) = 0.02476
The code is as follows:
clc;
clear all;
t=[ 0
10
20
30
40
50
60];
c=[ 100 100 0 0 0 0 0
94.77333333 93.73333333 0.400145927 0.27701074 0.107176434 0.158130584 0.000505722
84.48 88.87333333 0.604820791 0.42039699 0.148464431 0.359954181 0.001077632
73.91333333 83.23333333 0.589530533 0.461026019 0.17071451 0.508774341 0.001205872
68.80666667 76.43333333 0.529979 0.48559649 0.180350987 0.635601375 0.001720902
63.56666667 73.92 0.502282172 0.52315421 0.198918834 0.700114548 0.002064255
59.76666667 70.59333333 0.478273783 0.574401193 0.207380132 0.845269187 0.002895749
];
theta0=rand(11,1)*1E-3;
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c,zeros(size(theta0)));
fprintf(1,'\tRate Constants:\n')
for k1 = 1:numel(theta)
fprintf(1, '\t\tTheta(%2d) = %8.5f\n', k1, theta(k1))
end
tv = linspace(min(t), max(t));
Cfit = kinetics(theta, tv);
figure(1)
hd = plot(t, c, 'p');
for k1 = 1:size(c,2)
CV(k1,:) = hd(k1).Color;
hd(k1).MarkerFaceColor = CV(k1,:);
end
hold on
hlp = plot(tv, Cfit);
for k1 = 1:size(c,2)
hlp(k1).Color = CV(k1,:);
end
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, 'c(1)', 'c(2)', 'c(3)', 'c(4)', 'c(5)','c(6)','c(7)', 'Location','best');
function C=kinetics(theta,t)
c0=[100; 100; 0; 0; 0; 0; 0];
[T,Cv]=ode45(@DifEq,t,c0);
function dC=DifEq(t,c)
dcdt=zeros(7,1);
dcdt(1)=-theta(1)*c(1)-theta(2)*c(1)-theta(3)*c(1)-theta(9)*c(1)*c(2);
dcdt(2)=-theta(4)*c(2)*c(3)-theta(5)*c(2)*c(4)-theta(6)*c(2)*c(3)-theta(7)*c(2)*c(4)-theta(8)*c(2)-theta(9)*c(2)*c(1);
dcdt(3)= theta(2)*c(1)-theta(4)*c(2)*c(3)-theta(6)*c(2)*c(3);
dcdt(4)= theta(3)*c(1)-theta(5)*c(2)*c(4)-theta(7)*c(2)*c(4);
dcdt(5)= theta(4)*c(2)*c(3)-theta(10)*c(5)+theta(9)*c(2)*c(1);
dcdt(6)= theta(5)*c(2)*c(4)-theta(11)*c(6);
dcdt(7)= theta(6)*c(2)*c(3)+theta(7)*c(2)*c(4)+theta(8)*c(2);
end
C=Cv;
end
Accepted Answer
That does not appear to me to be the same data.
I get decent results except for 
and 
, leading me to believe that there is something wrong with their differential equations. Please check them.
and , leading me to believe that there is something wrong with their differential equations. Please check them. % clc;
% clear all;
t=[ 0
10
20
30
40
50
60];
c=[ 100 100 0 0 0 0 0
94.77333333 93.73333333 0.400145927 0.27701074 0.107176434 0.158130584 0.000505722
84.48 88.87333333 0.604820791 0.42039699 0.148464431 0.359954181 0.001077632
73.91333333 83.23333333 0.589530533 0.461026019 0.17071451 0.508774341 0.001205872
68.80666667 76.43333333 0.529979 0.48559649 0.180350987 0.635601375 0.001720902
63.56666667 73.92 0.502282172 0.52315421 0.198918834 0.700114548 0.002064255
59.76666667 70.59333333 0.478273783 0.574401193 0.207380132 0.845269187 0.002895749
];
theta0 = [rand(11,1)*1E-3; c(1,:).'];
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c,zeros(size(theta0)));
fprintf(1,'\tRate Constants:\n')
for k1 = 1:numel(theta)
fprintf(1, '\t\tTheta(%2d) = %8.5f\n', k1, theta(k1))
end
tv = linspace(min(t), max(t));
Cfit = kinetics(theta, tv);
figure(1)
hd = plot(t, c, 'p');
for k1 = 1:size(c,2)
CV(k1,:) = hd(k1).Color;
hd(k1).MarkerFaceColor = CV(k1,:);
end
hold on
hlp = plot(tv, Cfit);
for k1 = 1:size(c,2)
hlp(k1).Color = CV(k1,:);
end
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, 'c(1)', 'c(2)', 'c(3)', 'c(4)', 'c(5)','c(6)','c(7)', 'Location','best');
figure
for k = 1:size(c,2)
subplot(4,2,k)
plot(t,c(:,k),'bp', 'MarkerFaceColor','b')
hold on
plot(tv, Cfit(:,k))
hold off
grid
title("C_"+k)
end
function C=kinetics(theta,t)
% c0=[100; 100; 0; 0; 0; 0; 0];
c0 = theta(end-6:end);
[T,Cv]=ode45(@DifEq,t,c0);
function dcdt=DifEq(t,c)
dcdt=zeros(7,1);
dcdt(1)=-theta(1)*c(1)-theta(2)*c(1)-theta(3)*c(1)-theta(9)*c(1)*c(2);
dcdt(2)=-theta(4)*c(2)*c(3)-theta(5)*c(2)*c(4)-theta(6)*c(2)*c(3)-theta(7)*c(2)*c(4)-theta(8)*c(2)-theta(9)*c(2)*c(1);
dcdt(3)= theta(2)*c(1)-theta(4)*c(2)*c(3)-theta(6)*c(2)*c(3);
dcdt(4)= theta(3)*c(1)-theta(5)*c(2)*c(4)-theta(7)*c(2)*c(4);
dcdt(5)= theta(4)*c(2)*c(3)-theta(10)*c(5)+theta(9)*c(2)*c(1);
dcdt(6)= theta(5)*c(2)*c(4)-theta(11)*c(6);
dcdt(7)= theta(6)*c(2)*c(3)+theta(7)*c(2)*c(4)+theta(8)*c(2);
end
C=Cv;
end
I changed the code slightly to add the initial conditions as parameters to be estimated, since in my experience, that generally produces better results, even if the initial conditions themselves do not change much from their original values. (I added the subplot array because otherwise the last 5 compartments are squashed together at the lower end of the first plot.)
.
14 Comments
Kaihua Liu
on 30 Jun 2023
Star Strider
on 30 Jun 2023
As always, my pleasure!
I will help as much as I can if you have other problems with it (since I wrote the original version of the code you are using, and feel responsible for being certain that it works correctly). I do not know what you are calculating, and since I do not have access to those differential equations, I cannot check it myself.
Kaihua Liu
on 3 Aug 2023
I was in the process of responding to this when Micro$oft crashed my computer again with more of this azure and autopilot s**t. I am never buying a Windows machine again, and I am going to seriously consider archiving my files, wiping Win 11 from this machine, and installing Linux (that I have to learn about first, however that errort will be more than worth the effort).
Anyway, my intense hatred for Micro$oft aside, the only option fot the confidence intervals is nlpredci. (Use nlparci to get confidence intervals on the parameters.) I can implement that here by tweaking calls to it. I am not certain how to calculate the R² or other goodness-of-fit statistics for a regression such as this. (Occasonally an svd error message appears in the nlpredci call. If that occurs, just run the code again.)
Try this —
% clc;
% clear all;
t=[ 0
10
20
30
40
50
60];
c=[ 100 100 0 0 0 0 0
94.77333333 93.73333333 0.400145927 0.27701074 0.107176434 0.158130584 0.000505722
84.48 88.87333333 0.604820791 0.42039699 0.148464431 0.359954181 0.001077632
73.91333333 83.23333333 0.589530533 0.461026019 0.17071451 0.508774341 0.001205872
68.80666667 76.43333333 0.529979 0.48559649 0.180350987 0.635601375 0.001720902
63.56666667 73.92 0.502282172 0.52315421 0.198918834 0.700114548 0.002064255
59.76666667 70.59333333 0.478273783 0.574401193 0.207380132 0.845269187 0.002895749
];
theta0 = [rand(11,1)*1E-3; c(1,:).'];
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c,zeros(size(theta0)));
fprintf(1,'\tRate Constants:\n')
for k1 = 1:numel(theta)
fprintf(1, '\t\tTheta(%2d) = %8.5f\n', k1, theta(k1))
end
tv = linspace(min(t), max(t));
Cfit = kinetics(theta, tv);
for k = 1:size(c,2)
[Ypred,delta] = nlpredci(@(theta,t)kinetics1(theta,t,k),tv,theta,Rsd,'Jacobian',Jmat);
ypred(:,k) = Ypred;
CI(:,k) = delta;
end
figure(1)
hd = plot(t, c, 'p');
for k1 = 1:size(c,2)
CV(k1,:) = hd(k1).Color;
hd(k1).MarkerFaceColor = CV(k1,:);
end
hold on
hlp = plot(tv, Cfit);
for k1 = 1:size(c,2)
hlp(k1).Color = CV(k1,:);
end
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, 'c(1)', 'c(2)', 'c(3)', 'c(4)', 'c(5)','c(6)','c(7)', 'Location','best');
figure
for k = 1:size(c,2)
subplot(4,2,k)
hp1 = plot(t,c(:,k),'bp', 'MarkerFaceColor','b', 'DisplayName','Data');
hold on
hp2 = plot(tv, ypred(:,k), '-r', 'DisplayName','Regression');
hp3 = plot(tv, ypred(:,k)+CI(:,k)*[-1 1], ':r', 'DisplayName','±95% Conofidence Intervals');
hold off
grid
title("C_"+k)
end
hl = legend([hp1 hp2 hp3(1)]);
subplot(4,2,8)
Ax = gca;
Ax.Visible = 0;
hl.Position = Ax.Position;
function C=kinetics(theta,t)
% c0=[100; 100; 0; 0; 0; 0; 0];
c0 = theta(end-6:end);
[T,Cv]=ode45(@DifEq,t,c0);
function dcdt=DifEq(t,c)
dcdt=zeros(7,1);
dcdt(1)=-theta(1)*c(1)-theta(2)*c(1)-theta(3)*c(1)-theta(9)*c(1)*c(2);
dcdt(2)=-theta(4)*c(2)*c(3)-theta(5)*c(2)*c(4)-theta(6)*c(2)*c(3)-theta(7)*c(2)*c(4)-theta(8)*c(2)-theta(9)*c(2)*c(1);
dcdt(3)= theta(2)*c(1)-theta(4)*c(2)*c(3)-theta(6)*c(2)*c(3);
dcdt(4)= theta(3)*c(1)-theta(5)*c(2)*c(4)-theta(7)*c(2)*c(4);
dcdt(5)= theta(4)*c(2)*c(3)-theta(10)*c(5)+theta(9)*c(2)*c(1);
dcdt(6)= theta(5)*c(2)*c(4)-theta(11)*c(6);
dcdt(7)= theta(6)*c(2)*c(3)+theta(7)*c(2)*c(4)+theta(8)*c(2);
end
C=Cv;
end
function C=kinetics1(theta,t,k)
% c0=[100; 100; 0; 0; 0; 0; 0];
c0 = theta(end-6:end);
[T,Cv]=ode45(@DifEq,t,c0);
function dcdt=DifEq(t,c)
dcdt=zeros(7,1);
dcdt(1)=-theta(1)*c(1)-theta(2)*c(1)-theta(3)*c(1)-theta(9)*c(1)*c(2);
dcdt(2)=-theta(4)*c(2)*c(3)-theta(5)*c(2)*c(4)-theta(6)*c(2)*c(3)-theta(7)*c(2)*c(4)-theta(8)*c(2)-theta(9)*c(2)*c(1);
dcdt(3)= theta(2)*c(1)-theta(4)*c(2)*c(3)-theta(6)*c(2)*c(3);
dcdt(4)= theta(3)*c(1)-theta(5)*c(2)*c(4)-theta(7)*c(2)*c(4);
dcdt(5)= theta(4)*c(2)*c(3)-theta(10)*c(5)+theta(9)*c(2)*c(1);
dcdt(6)= theta(5)*c(2)*c(4)-theta(11)*c(6);
dcdt(7)= theta(6)*c(2)*c(3)+theta(7)*c(2)*c(4)+theta(8)*c(2);
end
C=Cv(:,k);
end
NOTE — I added the ‘kinetics1’ function since that is necessary to run the nlpredci function with individual vectors, since it does not accept matrix arguments. That also explains the for loop that calls it.
EDIT — Corrected typographical errors.
.
Kaihua Liu
on 4 Aug 2023
Star Strider
on 4 Aug 2023
That call would be:
ci = nlparci(theta,Rsd,'jacobian',Jmat)
The two-sided confidence intervals will be in the order of the parameter vector ‘theta’.
The parameters are significant if both confidence intervals have the same sign. If they are opposite signs, that means the 95% confidence intervals include zero, and that parameter is not necessary for the model.
.
Kaihua Liu
on 4 Aug 2023
Star Strider
on 4 Aug 2023
Edited: Star Strider
on 5 Aug 2023
The way I have them set (‘lb’ is zero, ‘ub’ not specified, so essentially infinite), all the parameter values are positive. Check to be certain that all the arguments are in the correct locations in the argument list. Otherwise, copy my lsqcurvefit call and it should work.
EDIT — (5 Aug 2023 at 16:59)
Thinking about this further, I believe that you are referring to the 95% confidence limits of the parameters, not the parameters themselves, none of which should be negative in this particular model.
- The parameters are significant if both confidence intervals have the same sign. If they are opposite signs, that means the 95% confidence intervals include zero, and that parameter is not necessary for the model.
So wide confidence intervals are occasionally to be expected, and the parameter value should be midway between them.
If the confidence intervals are of the same sign (either positive or negative), the parameter is significantly different from zero with 95% probability, and so contributes to the regression. If the confidence intervals are of opposite signs, that means that they include zero (that parameter is not statistically significantly different from zero with 95% probability), so that specific parameter is actually not necessary in the regression.
.
Kaihua Liu
on 20 Aug 2023
Star Strider
on 20 Aug 2023
Getting the standard errors for the parameters is likely going to require writing your own function (to either calculate them directly or back-calculate from the confidence intervals), and probably is not necessary anyway. The nlparci function (discussed earlier here and here) will give the 95% confidence intervals for the parameters, and that is actually the most significant information.
I agree that there is a problem with the code, particularly with 
and 
. Look at them specifically.
and . Look at them specifically.
Kaihua Liu
on 20 Aug 2023
Star Strider
on 20 Aug 2023
As always, my pleasure!
dinesh kumar s
on 10 Sep 2023
@Star Strider Hi, I was working on similar problem like this and ofcourse I was using the same code written by you for the kinetic modeling. I hope the parameters are good. I just want to thank you @Star Strider for the code. It works well and good. Thanks.
Star Strider
on 10 Sep 2023
@dinesh kumar s — My pleasure!
A Vote would be appreciated!
More Answers (0)
Categories
Find more on MATLAB in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
