seqcurvefit21.m

多元模型回归与分析-2 matlab 程序

function seqcurvefit21
clear all;
load seqdata;
[y_row,y_col]=size(c);

beta0=[0.005 0.001];
c0=[1 0];
lb=[0 0];ub=[inf inf];
[beta,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@seqfun,beta0,lb,ub,[],t,c,y_col,c0);

ci = nlparci(beta,residual,jacobian);

% print result
fprintf('\n Estimated Parameters by Lsqnonlin():\n')
fprintf('\t k1 = %.4f ± %.4f\n',beta(1),ci(1,2)-beta(1))
fprintf('\t k2 = %.4f ± %.4f\n',beta(2),ci(2,2)-beta(2))
fprintf('  The sum of the residual squares is: %.1e\n\n',sum(residual.^2))

% plot of fit results
tspan = [0  max(t)];
[tt yc] = ode45(@modeleqs,tspan,c0,[],beta);
tc=linspace(0,max(t),200);
yca = spline(tt,yc(:,1),tc);
plot(t,c(:,1),'ro',tc,yca,'r-');
hold on
ycb = spline(tt,yc(:,2),tc);
plot(t,c(:,2),'b+',tc,ycb,'b-');
xlabel('Time');
ylabel('Concentration');
hold off

function y = seqfun(beta,t,c,y_col,c0)
tspan = [0  max(t)];
[tt yy] = ode45(@modeleqs,tspan,c0,[],beta); 
for col = 1:y_col
    yc(:,col) = spline(tt,yy(:,col),t);
end
y=[c(:,1)-yc(:,1); c(:,2)-yc(:,2)];


function dydt = modeleqs(t,y,beta)                 % Model equations
dydt = [-beta(1)*y(1);
   beta(1)*y(1)-beta(2)*y(2)];