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

箫剑→残念

function parafit
2 L) c* |  V, A" W9 ~1 r' |. ?%  k1->k-1，k2->k1，k3->k2,k4->k3,k5->k4
% k6->k6 k7->k7
. w# t: O3 B+ f$ D9 f7 J% dGlcdt = k-1*C(Fru)-(k1+k2)*C(Glc);
% dFrudt = k1*C(Glc)-(k-1+k3+k4)C(Fru);
" r" A1 I  \6 E  G* ]% dFadt = k(2)*C(Glc)+k4*C(Fru)+（k6+k7)*C(Hmf);
% dLadt = k(7)*C(Hmf);
%dHmfdt = k(3)*C(Fru)-（k6+k7)*C(Hmf);
% R7 h( A9 o; u9 I$ M. y- |5 `clear all
" d6 o/ R9 P, i5 [5 O# Rclc
format long
! v1 h. Q% Q" f# l( ]# |%        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
  [" b0 c! u5 B- g/ Y- d7 c          30    0.19345    0.027    0.00868    0    7.00E-04
* M. t9 f% M$ d% X' H0 h$ q8 ^% z' Y2 @          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
$ y. K0 o5 o$ U1 s5 M, |$ K1 I          120    0.0656    0.07397    0.07885    0.00573    0.02143
          180    0.04488    0.0682    0.07135    0.0091    0.03623
- |6 b$ a* t  z' \  c* a          240    0.03653    0.06488    0.08945    0.01828    0.05452
          300    0.02738    0.05448    0.09098    0.0227    0.0597
8 \6 ?4 f3 i& j3 Z          360    0.01855    0.04125    0.09363    0.0239    0.06495];
) \4 H) @2 |& n# m  y  s) gk0 = [0.0000000005  0.0000000005  0.0000000005  0.00000000005  0.00005  0.0134  0.00564  0.00001  0.00001  0.00001];        % 参数初值
8 w6 Z; j7 f' qlb = [0  0  0  0  0  0  0  0  0  0];                  % 参数下限
ub = [1  1  1  1  1  1  1  1  1  1];    % 参数上限
& i6 V9 o, _/ H$ c4 a* g0 y  c; E/ \x0 = [0.25  0  0  0  0];
5 ]5 Q" g  i  R% t; H; R' ^! Byexp = Kinetics;                 % yexp: 实验数据[x1        x4        x5        x6]
) C2 h, Q( y/ a# m, U7 e% warning off
% 使用函数 ()进行参数估计
[k,fval,flag] = fmincon(@ObjFunc7Fmincon,k0,[],[],[],[],lb,ub,[],[],x0,yexp);
; K) i( b+ t1 y5 g1 W2 p# Ofprintf('\n使用函数fmincon()估计得到的参数值为:\n')
8 w+ @2 _, P& v+ S; Xfprintf('\tk1 = %.11f\n',k(1))
fprintf('\tk2 = %.11f\n',k(2))
6 a9 m& Q. ]) N- Ifprintf('\tk3 = %.11f\n',k(3))
! e+ r) p  ]5 p7 m9 ]; Mfprintf('\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))
fprintf('  The sum of the squares is: %.1e\n\n',fval)
( x& ?/ [- l9 R+ v% ak_fm= k;
% warning off
$ L" f; V+ B4 o6 }# w' a! h- i( a% 使用函数lsqnonlin()进行参数估计
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
6 g/ \& n% Z9 d( I" ~1 m3 g    lsqnonlin(@ObjFunc7LNL,k0,lb,ub,[],x0,yexp);      
$ E  R1 ?7 K1 a7 L& s: b3 x1 b8 xci = nlparci(k,residual,jacobian);
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
' [- R1 v$ V, U+ ?fprintf('\tk1 = %.11f\n',k(1))
2 N4 |% T" z; k! \7 j. mfprintf('\tk2 = %.11f\n',k(2))
) T! }8 j- }9 J( D9 Kfprintf('\tk3 = %.11f\n',k(3))
- }* t& e: U' ufprintf('\tk4 = %.11f\n',k(4))
fprintf('\tk5 = %.11f\n',k(5))
fprintf('\tk6 = %.11f\n',k(6))
' e, X# d1 ~) j& K& [+ Y# Efprintf('\tk7 = %.11f\n',k(7))
+ q5 X3 [5 ^+ r8 n3 H: U9 [7 nfprintf('\tk8 = %.11f\n',k(8))
" |* P  i* G+ g$ efprintf('\tk9 = %.11f\n',k(9))
fprintf('\tk10 = %.11f\n',k(10))
9 U3 T0 h7 \; e: pfprintf('  The sum of the squares is: %.1e\n\n',resnorm)
k_ls = k;
output
warning off
% 以函数fmincon()估计得到的结果为初值，使用函数lsqnonlin()进行参数估计
8 @5 P) y. u  m2 b0 Y3 Zk0 = k_fm;
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc7LNL,k0,lb,ub,[],x0,yexp);      
ci = nlparci(k,residual,jacobian);
  Q# K9 L) J" \" Z5 r0 ffprintf('\n\n以fmincon()的结果为初值，使用函数lsqnonlin()估计得到的参数值为:\n')
! `8 G8 s8 Y/ n% K5 U) [( xfprintf('\tk1 = %.11f\n',k(1))
6 z9 i9 N: x- ^* _fprintf('\tk2 = %.11f\n',k(2))
; K! t+ x* a5 p( j5 ^) ifprintf('\tk3 = %.11f\n',k(3))
7 g3 D3 a+ I  x/ t# |9 Efprintf('\tk4 = %.11f\n',k(4))
fprintf('\tk5 = %.11f\n',k(5))
fprintf('\tk6 = %.11f\n',k(6))
# {$ l# P/ \. l0 Q5 g1 cfprintf('\tk7 = %.11f\n',k(7))
' S7 L4 E$ v% Mfprintf('\tk8 = %.11f\n',k(8))
- ~$ J5 G" O( ^" U) wfprintf('\tk9 = %.11f\n',k(9))
fprintf('\tk10 = %.11f\n',k(10))
fprintf('  The sum of the squares is: %.1e\n\n',resnorm)
k_fmls = k;
output
tspan = [0 15 30 45 60 90 120 180 240 300 360];
' [( t& _, L  K4 x[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 g+ E  u1 G1 W. Q5 whold on
plot(t,yexp(:,3:6),'o');legend('Fru-pr','Fa-pr','La-pr','HMF-pr','Fru-real','Fa-real','La-real','HMF-real')



function f = ObjFunc7LNL(k,x0,yexp)
% k+ T" z0 p8 F% otspan = [0 15 30 45 60 90 120 180 240 300 360];
+ x, h! ]9 h9 S! X# ]4 c# R9 c[t, x] = ode45(@KineticEqs,tspan,x0,[],k);   
y(:,2) = x(:,1);
y(:,3:6) = x(:,2:5);
f1 = y(:,2) - yexp(:,2);
  O; H' _, @% v* i8 H4 r. O/ Df2 = y(:,3) - yexp(:,3);
* ?6 B7 W) N' E7 {% l3 ~f3 = y(:,4) - yexp(:,4);
f4 = y(:,5) - yexp(:,5);
1 R3 M' G5 t; h+ `/ `* @/ h8 wf5 = y(:,6) - yexp(:,6);
) P. f, m! ~) C. f: ]7 \5 Bf = [f1; f2; f3; f4; f5];



6 X( C' G$ G# |7 R% ?% S' ]function f = ObjFunc7Fmincon(k,x0,yexp)
* L" Q, R- Q( y) b% k- V* Ttspan = [0 15 30 45 60 90 120 180 240 300 360];
[t x] = ode45(@KineticEqs,tspan,x0,[],k);   
y(:,2) = x(:,1);
% Y+ }( n7 A+ t2 L& Qy(:,3:6) = x(:,2:5);
5 a7 c3 E' [+ ^7 V  af =  sum((y(:,2)-yexp(:,2)).^2) + sum((y(:,3)-yexp(:,3)).^2)   ...
    + sum((y(:,4)-yexp(:,4)).^2) + sum((y(:,5)-yexp(:,5)).^2)   ...
2 G+ X* V! \& Z9 \6 v6 |0 q/ m    + sum((y(:,6)-yexp(:,6)).^2) ;
" K$ g7 K- F: J0 `$ t

8 f( ^" l0 E, s5 \' @. w

function dxdt = KineticEqs(t,x,k)
0 X2 w* }  _2 Y- ?  b9 SdGldt = 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);
: j& g; o' y1 Y/ m+ ZdLadt = k(7)*x(5);
; _5 Y7 w* k) ~- jdHmdt = k(4)*x(2)-(k(6)+k(7)+k(10))*x(5);
2 Y% O7 }( n$ v+ q$ w* R' q2 ~6 Ddxdt = [dGldt; dFrdt; dFadt; dLadt; dHmdt];

" M  s% X& }  U* ~  _+ ]  G/ P