function k1k2k32
format long
clear all
clc
tspan = [0 6 24 44 68 72 74 92 104 116]';%%这是时间
yexp= [3.111,3.639,3.887,4.289,4.658,5.531,6.218,6.979,7.111,7.114]';%%%这是菌落总数
x0 = [3.111];
k0 = [0.1345 4 37 8.1];
lb = [0000 00];
ub = [100100100 100 100 100];
%
% % opts = statset('nlinfit');
% % opts.RobustWgtFun = 'bisquare';
% % mdll = fitnlm(tspan,yexp,ObjFunc,b0,'Options')
% mdll = fitnlm(tspan,yexp,ObjFunc,b0,'Options',opts)
options = optimoptions(@lsqnonlin,'Algorithm','trust-region-reflective');
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
lsqnonlin(@ObjFunc,k0,lb,ub,[],tspan,x0,yexp);
ci = nlparci(k,residual,jacobian);
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\tk1 = %.9f ± %.9f\n',k(1),ci(1,2)-k(1))
fprintf('\tk2 = %.9f ± %.9f\n',k(2),ci(2,2)-k(2))
fprintf('\tk3 = %.9f ± %.9f\n',k(3),ci(3,2)-k(3))
fprintf('\tk4 = %.9f ± %.9f\n',k(4),ci(4,2)-k(4))
fprintf('\tk5 = %.9f ± %.9f\n',k(5),ci(5,2)-k(5))
fprintf('\tk6 = %.9f ± %.9f\n',k(6),ci(6,2)-k(6))
fprintf('The sum of the squares is: %.9e\n\n',resnorm)
tsa=0:0.01:max(tspan);
[tsa ysa]=ode45(@KineticsEqs,tsa,x0,[],k);
figure(1),
plot(tsa,ysa(:,1),'b',tspan,yexp,'or'),legend('计算值','实验值','Location','best');
function f = ObjFunc(k,tspan,x0,yexp) % 目标函数
[tspany] = ode45(@KineticsEqs,tspan,x0,[],k);
f = (y-yexp);
function dCdt = KineticsEqs(tspan ,C,k)
T=4*(tspan >=0 & tspan <=30)+10*(tspan >31 & tspan <=90)+4*(tspan >90 ); % 假设这是我写的温度关于时间的分段函数
umax=((k(2)*(T-k(3))*(T-k(4))*(T-k(4)))/(((k(5)-k(4))*(T-k(5))-(k(5)-k(3))*(k(5)+k(4)-2*T))*(k(5)-k(4))));
dCdt=(1/(1+exp(-4*(tspan -k(1)))))*umax*(1-exp(C-k(6)));
这是根据他人的代码更改的,得到的结果是置信区间过大?想有什么办法解决?以及想用fitnlm函数进行拟合,达到一些参数评价t值p值等,但用fitnlm运行不起来。 @beefly