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

箫剑→残念

function parafit
' B, v' C# w) m  ^  {- x%  k1->k-1，k2->k1，k3->k2,k4->k3,k5->k4
9 L2 ]5 i/ F! }- ]% k6->k6 k7->k7
% 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);
% dLadt = k(7)*C(Hmf);
0 M! ]' k; s5 ^%dHmfdt = k(3)*C(Fru)-（k6+k7)*C(Hmf);
8 P: n( Q* S  I" i  ]( lclear all
; W; R. b6 n' S/ O6 ~1 v4 |clc
format long
%        t/min   Glc    Fru        Fa   La   HMF/ mol/L
  Kinetics=[0    0.25    0           0    0       0
          15    0.2319    0.01257    0.0048    0    2.50E-04
# K, k0 A- p2 E8 J/ y- D% g          30    0.19345    0.027    0.00868    0    7.00E-04
          45    0.15105    0.06975    0.02473    0    0.0033
8 }! Q- Z8 ]1 J- Q          60    0.13763    0.07397    0.02615    0    0.00428
          90    0.08115    0.07877    0.07485    0    0.01405
! o4 J' h1 O: _8 f, ^) A2 c- E4 K3 v          120    0.0656    0.07397    0.07885    0.00573    0.02143
' _) M5 M* O0 ?1 z& o; {          180    0.04488    0.0682    0.07135    0.0091    0.03623
          240    0.03653    0.06488    0.08945    0.01828    0.05452
# g. ~* W( ?+ L1 S8 b' k          300    0.02738    0.05448    0.09098    0.0227    0.0597
          360    0.01855    0.04125    0.09363    0.0239    0.06495];
k0 = [0.0000000005  0.0000000005  0.0000000005  0.00000000005  0.00005  0.0134  0.00564  0.00001  0.00001  0.00001];        % 参数初值
4 ^% T4 t* \' z  E( f3 |; ulb = [0  0  0  0  0  0  0  0  0  0];                  % 参数下限
ub = [1  1  1  1  1  1  1  1  1  1];    % 参数上限
6 h5 m; A! Q1 A, Z' }7 {) y) k/ ax0 = [0.25  0  0  0  0];
yexp = Kinetics;                 % yexp: 实验数据[x1        x4        x5        x6]
+ D" w3 y1 i( L+ R% warning off
% 使用函数 ()进行参数估计
: g) t& Q/ c7 t& I  p; L[k,fval,flag] = fmincon(@ObjFunc7Fmincon,k0,[],[],[],[],lb,ub,[],[],x0,yexp);
fprintf('\n使用函数fmincon()估计得到的参数值为:\n')
fprintf('\tk1 = %.11f\n',k(1))
fprintf('\tk2 = %.11f\n',k(2))
fprintf('\tk3 = %.11f\n',k(3))
* u  O: J- G2 ~- nfprintf('\tk4 = %.11f\n',k(4))
fprintf('\tk5 = %.11f\n',k(5))
fprintf('\tk6 = %.11f\n',k(6))
% q/ j% [& {. B3 y% Z! Bfprintf('\tk7 = %.11f\n',k(7))
2 Q# w; ?. S' P- i* Z1 ^  ^* n3 g. `, Ufprintf('\tk8 = %.11f\n',k(8))
  I0 a7 s6 A+ a/ {5 |; c% f. Sfprintf('\tk9 = %.11f\n',k(9))
2 ~: P* _- W; M8 b: A. xfprintf('\tk10 = %.11f\n',k(10))
fprintf('  The sum of the squares is: %.1e\n\n',fval)
6 C1 m5 t) O% P9 mk_fm= k;
8 \+ K) A8 _' c- e) P, ]% warning off
0 @2 w8 w, \; `, U; `( q. }% 使用函数lsqnonlin()进行参数估计
- h+ m1 ?$ n  N: g: t# s[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
2 U! w+ x& M) M" B( S    lsqnonlin(@ObjFunc7LNL,k0,lb,ub,[],x0,yexp);      
ci = nlparci(k,residual,jacobian);
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
% n7 C" A: k9 ?- Dfprintf('\tk1 = %.11f\n',k(1))
5 p7 u% C! s3 P1 @" Mfprintf('\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))
0 |# \# y1 `* m  Nfprintf('\tk9 = %.11f\n',k(9))
: t+ [# W, m6 f. T7 c2 ?2 yfprintf('\tk10 = %.11f\n',k(10))
fprintf('  The sum of the squares is: %.1e\n\n',resnorm)
k_ls = k;
, ~6 w$ r8 K; r. H9 [& routput
4 _( j% ^# \- O1 O: ywarning off
% 以函数fmincon()估计得到的结果为初值，使用函数lsqnonlin()进行参数估计
# @( X; _5 ^" e* d$ H: O4 r  Qk0 = k_fm;
% w2 h& w1 z8 ]7 s) r[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
4 K- G" E9 z: B4 I6 M  ~* s6 J4 M: b    lsqnonlin(@ObjFunc7LNL,k0,lb,ub,[],x0,yexp);      
+ }# h3 F8 c3 vci = nlparci(k,residual,jacobian);
' _, q! q3 t  M% Xfprintf('\n\n以fmincon()的结果为初值，使用函数lsqnonlin()估计得到的参数值为:\n')
9 ?! d% G) p8 W; c) U% j. _fprintf('\tk1 = %.11f\n',k(1))
fprintf('\tk2 = %.11f\n',k(2))
* X$ T; Z* ?; C, J' A2 _fprintf('\tk3 = %.11f\n',k(3))
* _" W/ K  a9 e  g, gfprintf('\tk4 = %.11f\n',k(4))
9 l& H% u3 ?  G' j5 Qfprintf('\tk5 = %.11f\n',k(5))
0 \- A" n4 l' d) qfprintf('\tk6 = %.11f\n',k(6))
fprintf('\tk7 = %.11f\n',k(7))
9 e( b* e* w6 Q0 D# j/ ?fprintf('\tk8 = %.11f\n',k(8))
fprintf('\tk9 = %.11f\n',k(9))
& c- C) E( D! {9 k& qfprintf('\tk10 = %.11f\n',k(10))
1 E0 H2 ]! e9 [fprintf('  The sum of the squares is: %.1e\n\n',resnorm)
k_fmls = k;
output
8 Q% D& L  p& ]* n) ztspan = [0 15 30 45 60 90 120 180 240 300 360];
3 y2 [  n; u4 Z; z, m# @0 P. V[t x] = ode45(@KineticEqs,tspan,x0,[],k_fmls);
figure;
plot(t,x(:,1),t,yexp(:,2),'*');legend('Glc-pr','Glc-real')
figure;plot(t,x(:,2:5));
p=x(:,1:5)
2 Y. w! B* O2 p' j; o# b) Ghold on
plot(t,yexp(:,3:6),'o');legend('Fru-pr','Fa-pr','La-pr','HMF-pr','Fru-real','Fa-real','La-real','HMF-real')


3 c! E6 K/ I# ^( u( I: ~$ G% g0 N$ H, q
3 _& Z: A( N; J" `function f = ObjFunc7LNL(k,x0,yexp)
( W: V8 ]* p5 b- ztspan = [0 15 30 45 60 90 120 180 240 300 360];
[t, x] = ode45(@KineticEqs,tspan,x0,[],k);   
y(:,2) = x(:,1);
y(:,3:6) = x(:,2:5);
f1 = y(:,2) - yexp(:,2);
, y: o# |/ B. S& E- T- \0 O( Jf2 = y(:,3) - yexp(:,3);
  H# K6 ^3 G9 g2 lf3 = y(:,4) - yexp(:,4);
f4 = y(:,5) - yexp(:,5);
f5 = y(:,6) - yexp(:,6);
7 m  Q/ Q8 n8 a. u  W4 @3 G( xf = [f1; f2; f3; f4; f5];

8 _3 s9 W+ _1 n2 H1 w
" u7 K* k: H) |: O
function f = ObjFunc7Fmincon(k,x0,yexp)
tspan = [0 15 30 45 60 90 120 180 240 300 360];
* }8 H6 y) g# s) f% F) [[t x] = ode45(@KineticEqs,tspan,x0,[],k);   
+ R! U2 H8 Y, P& y/ B3 T5 s) h* F# d4 oy(:,2) = x(:,1);
y(:,3:6) = x(:,2:5);
- `1 c, ^  |8 zf =  sum((y(:,2)-yexp(:,2)).^2) + sum((y(:,3)-yexp(:,3)).^2)   ...
6 A9 v, K) \! w; M3 w    + sum((y(:,4)-yexp(:,4)).^2) + sum((y(:,5)-yexp(:,5)).^2)   ...
6 D9 Q& o, N/ j% h0 M' y& D/ h    + sum((y(:,6)-yexp(:,6)).^2) ;
( v* G  l: m: B0 O

" X0 M$ N+ k: t8 y

function dxdt = KineticEqs(t,x,k)
dGldt = k(1)*x(2)-(k(2)+k(3)+k(8))*x(1);
dFrdt = k(2)*x(1)-(k(1)+k(4)+k(5)+k(9))*x(2);
dFadt = k(3)*x(1)+k(5)*x(2)+(k(6)+k(7))*x(5);
dLadt = k(7)*x(5);
dHmdt = k(4)*x(2)-(k(6)+k(7)+k(10))*x(5);
2 }" H6 M, m/ B* R( y* G' \2 Ddxdt = [dGldt; dFrdt; dFadt; dLadt; dHmdt];
" h. G& |* {) ]0 P