帮忙做下统计显著性检验和K值的误差以及灵敏度分析

箫剑→残念

function parafit
%  k1->k-1，k2->k1，k3->k2,k4->k3,k5->k4
% k6->k6 k7->k7
2 F1 G+ b, a, `9 ^: D! D$ |% dGlcdt = k-1*C(Fru)-(k1+k2)*C(Glc);
% dFrudt = k1*C(Glc)-(k-1+k3+k4)C(Fru);
% dFadt = k(2)*C(Glc)+k4*C(Fru)+（k6+k7)*C(Hmf);
/ r4 d* [- V4 S! E% m" G% @% dLadt = k(7)*C(Hmf);
. r/ g/ `+ U$ K2 Y+ M%dHmfdt = k(3)*C(Fru)-（k6+k7)*C(Hmf);
! }) _; W+ Y6 D1 L% Xclear all
clc
5 {- C4 D; o* Pformat long
4 A. u6 D9 t: K+ @! u9 d%        t/min   Glc    Fru        Fa   La   HMF/ mol/L
  Kinetics=[0    0.25    0           0    0       0
  s/ ~( k3 c4 T" c          15    0.2319    0.01257    0.0048    0    2.50E-04
          30    0.19345    0.027    0.00868    0    7.00E-04
, _7 F/ a: y2 Z          45    0.15105    0.06975    0.02473    0    0.0033
          60    0.13763    0.07397    0.02615    0    0.00428
          90    0.08115    0.07877    0.07485    0    0.01405
, v7 T+ s5 k! ?          120    0.0656    0.07397    0.07885    0.00573    0.02143
3 k% X. `6 h; C& I          180    0.04488    0.0682    0.07135    0.0091    0.03623
: @% \5 Y& v4 T. b          240    0.03653    0.06488    0.08945    0.01828    0.05452
          300    0.02738    0.05448    0.09098    0.0227    0.0597
          360    0.01855    0.04125    0.09363    0.0239    0.06495];
9 w8 g' |. G0 N* Wk0 = [0.0000000005  0.0000000005  0.0000000005  0.00000000005  0.00005  0.0134  0.00564  0.00001  0.00001  0.00001];        % 参数初值
lb = [0  0  0  0  0  0  0  0  0  0];                  % 参数下限
ub = [1  1  1  1  1  1  1  1  1  1];    % 参数上限
5 M# N- _  `: a& f$ B* z% r* @x0 = [0.25  0  0  0  0];
3 D. s$ c1 o: v) ^' l5 P( Oyexp = Kinetics;                 % yexp: 实验数据[x1        x4        x5        x6]
; m, R& ]2 {9 L# n/ \% warning off
/ X; h1 }% k$ Z$ E% 使用函数 ()进行参数估计
- T$ }& g, k+ K0 {1 N+ J[k,fval,flag] = fmincon(@ObjFunc7Fmincon,k0,[],[],[],[],lb,ub,[],[],x0,yexp);
: K, h$ b2 j  Jfprintf('\n使用函数fmincon()估计得到的参数值为:\n')
fprintf('\tk1 = %.11f\n',k(1))
% f/ Z8 Y4 a/ Q* ifprintf('\tk2 = %.11f\n',k(2))
fprintf('\tk3 = %.11f\n',k(3))
fprintf('\tk4 = %.11f\n',k(4))
fprintf('\tk5 = %.11f\n',k(5))
fprintf('\tk6 = %.11f\n',k(6))
fprintf('\tk7 = %.11f\n',k(7))
fprintf('\tk8 = %.11f\n',k(8))
fprintf('\tk9 = %.11f\n',k(9))
fprintf('\tk10 = %.11f\n',k(10))
5 m: n1 v9 Q- V  E, Bfprintf('  The sum of the squares is: %.1e\n\n',fval)
& m9 ]$ T  r! K1 A2 K; Z7 C' ?7 Ck_fm= k;
3 O# P" r4 o3 Z% warning off
% 使用函数lsqnonlin()进行参数估计
; U5 T+ V4 x7 Y/ @* u/ N[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
; h2 V. X; m) z; E) Z# `' U; H0 ^    lsqnonlin(@ObjFunc7LNL,k0,lb,ub,[],x0,yexp);      
+ H; X7 b: v. x: o- \3 tci = nlparci(k,residual,jacobian);
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\tk1 = %.11f\n',k(1))
* E- i( q6 o5 b1 B0 I  dfprintf('\tk2 = %.11f\n',k(2))
- H9 n5 G6 H, k3 n/ Dfprintf('\tk3 = %.11f\n',k(3))
fprintf('\tk4 = %.11f\n',k(4))
fprintf('\tk5 = %.11f\n',k(5))
fprintf('\tk6 = %.11f\n',k(6))
fprintf('\tk7 = %.11f\n',k(7))
fprintf('\tk8 = %.11f\n',k(8))
fprintf('\tk9 = %.11f\n',k(9))
fprintf('\tk10 = %.11f\n',k(10))
9 E+ i8 |5 ^1 ~) T. `0 G7 A, `. f# Q" Qfprintf('  The sum of the squares is: %.1e\n\n',resnorm)
8 E- ~; C5 V9 ak_ls = k;
output
* u, l) _, G% x9 ?6 n0 _9 xwarning off
% 以函数fmincon()估计得到的结果为初值，使用函数lsqnonlin()进行参数估计
k0 = k_fm;
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
# S; ~$ k1 S5 W: q) T" j: N    lsqnonlin(@ObjFunc7LNL,k0,lb,ub,[],x0,yexp);      
3 L/ E0 d1 q/ O& K1 r2 R9 Ici = nlparci(k,residual,jacobian);
fprintf('\n\n以fmincon()的结果为初值，使用函数lsqnonlin()估计得到的参数值为:\n')
2 Y# A) k! p) W9 N' ifprintf('\tk1 = %.11f\n',k(1))
9 O) s+ s; D$ R* `& x, lfprintf('\tk2 = %.11f\n',k(2))
fprintf('\tk3 = %.11f\n',k(3))
fprintf('\tk4 = %.11f\n',k(4))
fprintf('\tk5 = %.11f\n',k(5))
- q; z$ n5 c! C( c& nfprintf('\tk6 = %.11f\n',k(6))
fprintf('\tk7 = %.11f\n',k(7))
fprintf('\tk8 = %.11f\n',k(8))
6 M' d( a% y+ h- H6 t* I8 mfprintf('\tk9 = %.11f\n',k(9))
: S' }( Z6 K6 k5 ?( Jfprintf('\tk10 = %.11f\n',k(10))
fprintf('  The sum of the squares is: %.1e\n\n',resnorm)
k_fmls = k;
0 T( J. |5 l5 p# E4 E% i, F# ?2 Qoutput
& `. z7 s+ ~& _tspan = [0 15 30 45 60 90 120 180 240 300 360];
[t x] = ode45(@KineticEqs,tspan,x0,[],k_fmls);
) K3 m# F; U5 J( Zfigure;
plot(t,x(:,1),t,yexp(:,2),'*');legend('Glc-pr','Glc-real')
figure;plot(t,x(:,2:5));
p=x(:,1:5)
" }+ Q: B* p5 [: X8 @hold on
) z6 d; D3 g) e8 L; }" ?: j; qplot(t,yexp(:,3:6),'o');legend('Fru-pr','Fa-pr','La-pr','HMF-pr','Fru-real','Fa-real','La-real','HMF-real')
! P. Z+ E  `- K' \0 ~
& O; I9 M2 }  O$ Q
: s- f: G/ @) O8 p
function f = ObjFunc7LNL(k,x0,yexp)
% h; q% R" A- ?5 b/ Xtspan = [0 15 30 45 60 90 120 180 240 300 360];
9 p& `9 Z3 t% p$ g; Y$ ?[t, x] = ode45(@KineticEqs,tspan,x0,[],k);   
y(:,2) = x(:,1);
y(:,3:6) = x(:,2:5);
+ c, |$ @" s! z! J+ Nf1 = y(:,2) - yexp(:,2);
f2 = y(:,3) - yexp(:,3);
% u" X9 H7 h+ h5 jf3 = y(:,4) - yexp(:,4);
% F( `1 R" e8 C$ Lf4 = y(:,5) - yexp(:,5);
) t" `0 R2 a5 w4 r* ]# of5 = y(:,6) - yexp(:,6);
# e( U9 i. w- z7 r  W$ e$ @f = [f1; f2; f3; f4; f5];

, g& }+ {+ t5 C4 J! y3 [* B

function f = ObjFunc7Fmincon(k,x0,yexp)
tspan = [0 15 30 45 60 90 120 180 240 300 360];
9 E1 n$ Z) q, C* X0 F[t x] = ode45(@KineticEqs,tspan,x0,[],k);   
* u: f) m/ N8 n3 m. Dy(:,2) = x(:,1);
y(:,3:6) = x(:,2:5);
f =  sum((y(:,2)-yexp(:,2)).^2) + sum((y(:,3)-yexp(:,3)).^2)   ...
    + sum((y(:,4)-yexp(:,4)).^2) + sum((y(:,5)-yexp(:,5)).^2)   ...
0 M. G: b" Q/ b+ c0 ^( D    + sum((y(:,6)-yexp(:,6)).^2) ;




/ s3 G5 C: B' g) D$ lfunction dxdt = KineticEqs(t,x,k)
' ?6 X# j6 V' ]dGldt = k(1)*x(2)-(k(2)+k(3)+k(8))*x(1);
! z6 d! O  ]# b8 j+ c; ]7 z2 EdFrdt = k(2)*x(1)-(k(1)+k(4)+k(5)+k(9))*x(2);
, s. f1 Q2 B8 @7 J( D, G2 ~, Z1 pdFadt = k(3)*x(1)+k(5)*x(2)+(k(6)+k(7))*x(5);
dLadt = k(7)*x(5);
# w& c7 R% z4 M' RdHmdt = k(4)*x(2)-(k(6)+k(7)+k(10))*x(5);
- U9 l% o/ ]1 B. kdxdt = [dGldt; dFrdt; dFadt; dLadt; dHmdt];
; R* G! N# v+ D: Z% {
6 ^; L5 c% h( O# M! G! S
3 G' G, z' F- t3 l$ M5 M