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

箫剑→残念

function parafit
%  k1->k-1，k2->k1，k3->k2,k4->k3,k5->k4
% k6->k6 k7->k7
) [8 w, M6 f& z. L) p# W) O% dGlcdt = k-1*C(Fru)-(k1+k2)*C(Glc);
% dFrudt = k1*C(Glc)-(k-1+k3+k4)C(Fru);
* r$ ~6 [2 X% y% P/ ?% dFadt = k(2)*C(Glc)+k4*C(Fru)+（k6+k7)*C(Hmf);
+ I2 C3 Q3 H  g+ T0 j! L; u/ r% dLadt = k(7)*C(Hmf);
# P9 u  N1 |% _% o6 N2 c. N%dHmfdt = k(3)*C(Fru)-（k6+k7)*C(Hmf);
clear all
. |# N; `' r. h2 O5 w; uclc
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
; g8 l: d+ H$ z9 _& X- d          30    0.19345    0.027    0.00868    0    7.00E-04
( ?% k1 J8 N3 ~" X* m          45    0.15105    0.06975    0.02473    0    0.0033
  s2 g. ]( h% H  E1 E( ]          60    0.13763    0.07397    0.02615    0    0.00428
- s1 i1 e. W! u; z          90    0.08115    0.07877    0.07485    0    0.01405
          120    0.0656    0.07397    0.07885    0.00573    0.02143
          180    0.04488    0.0682    0.07135    0.0091    0.03623
          240    0.03653    0.06488    0.08945    0.01828    0.05452
          300    0.02738    0.05448    0.09098    0.0227    0.0597
9 u/ y7 i5 B9 o% R  d6 U/ f: J7 P" F          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];        % 参数初值
lb = [0  0  0  0  0  0  0  0  0  0];                  % 参数下限
7 ~; n: C- u2 z+ j' E: C  p# Eub = [1  1  1  1  1  1  1  1  1  1];    % 参数上限
x0 = [0.25  0  0  0  0];
5 H( k# t3 t: d8 `; s/ F) g) X2 Gyexp = Kinetics;                 % yexp: 实验数据[x1        x4        x5        x6]
% warning off
% 使用函数 ()进行参数估计
) T0 q1 t% m, Z' J3 q% [* ~[k,fval,flag] = fmincon(@ObjFunc7Fmincon,k0,[],[],[],[],lb,ub,[],[],x0,yexp);
fprintf('\n使用函数fmincon()估计得到的参数值为:\n')
* y3 w8 p; k2 c$ E2 S' Kfprintf('\tk1 = %.11f\n',k(1))
) }# U" s9 M, I0 ?0 J& p3 L* ~8 _fprintf('\tk2 = %.11f\n',k(2))
6 B/ g$ X7 x) wfprintf('\tk3 = %.11f\n',k(3))
fprintf('\tk4 = %.11f\n',k(4))
3 v+ I, k* R) z% T: C- Q$ gfprintf('\tk5 = %.11f\n',k(5))
1 F% g8 X: c) y: [. L1 o! r- Kfprintf('\tk6 = %.11f\n',k(6))
fprintf('\tk7 = %.11f\n',k(7))
fprintf('\tk8 = %.11f\n',k(8))
0 _% r+ g% A& o; D. k1 F: ^fprintf('\tk9 = %.11f\n',k(9))
  Z; D) T4 F$ `fprintf('\tk10 = %.11f\n',k(10))
fprintf('  The sum of the squares is: %.1e\n\n',fval)
' E9 X2 Y7 W: \1 Ck_fm= k;
% warning off
% 使用函数lsqnonlin()进行参数估计
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc7LNL,k0,lb,ub,[],x0,yexp);      
# R8 K' i' ^8 M' G! r) p. J3 R, Ici = nlparci(k,residual,jacobian);
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\tk1 = %.11f\n',k(1))
1 f1 ^& I( p$ K9 m/ M: ufprintf('\tk2 = %.11f\n',k(2))
fprintf('\tk3 = %.11f\n',k(3))
; ~2 Y! N* X( f* l4 ^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))
! t- [. [0 d( r! X$ bfprintf('\tk9 = %.11f\n',k(9))
. P& m% r2 {: K7 w8 A  zfprintf('\tk10 = %.11f\n',k(10))
+ \; R" O0 X% Gfprintf('  The sum of the squares is: %.1e\n\n',resnorm)
k_ls = k;
output
! ~, n: i  r( N7 {1 fwarning off
! [6 V" u* E8 Y$ a* c# z/ _1 u. |% 以函数fmincon()估计得到的结果为初值，使用函数lsqnonlin()进行参数估计
- H6 Z8 [/ o) X" qk0 = k_fm;
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
+ s; |) ^3 w4 {* A8 o3 z    lsqnonlin(@ObjFunc7LNL,k0,lb,ub,[],x0,yexp);      
ci = nlparci(k,residual,jacobian);
. {* `, [9 ]4 z& q+ Q% O8 x* p- Ufprintf('\n\n以fmincon()的结果为初值，使用函数lsqnonlin()估计得到的参数值为:\n')
8 w: p' |7 u. ]5 F3 ^+ jfprintf('\tk1 = %.11f\n',k(1))
fprintf('\tk2 = %.11f\n',k(2))
fprintf('\tk3 = %.11f\n',k(3))
fprintf('\tk4 = %.11f\n',k(4))
+ W# j, G* y* l, Vfprintf('\tk5 = %.11f\n',k(5))
2 E1 e1 ^$ Q  g# t" jfprintf('\tk6 = %.11f\n',k(6))
fprintf('\tk7 = %.11f\n',k(7))
fprintf('\tk8 = %.11f\n',k(8))
fprintf('\tk9 = %.11f\n',k(9))
! z  q$ Y0 ]& |1 zfprintf('\tk10 = %.11f\n',k(10))
$ R3 [8 Q, |  X" X: S. Qfprintf('  The sum of the squares is: %.1e\n\n',resnorm)
k_fmls = k;
. s+ @0 f7 o7 R8 P8 I% Boutput
tspan = [0 15 30 45 60 90 120 180 240 300 360];
[t x] = ode45(@KineticEqs,tspan,x0,[],k_fmls);
figure;
8 r- |3 i& [: K( o( s$ ~plot(t,x(:,1),t,yexp(:,2),'*');legend('Glc-pr','Glc-real')
figure;plot(t,x(:,2:5));
p=x(:,1:5)
. O: Q- P$ _$ j4 A. P' [5 t, X7 m$ Lhold on
, {' F% a8 d% w: Pplot(t,yexp(:,3:6),'o');legend('Fru-pr','Fa-pr','La-pr','HMF-pr','Fru-real','Fa-real','La-real','HMF-real')

# ^( K( D) {) p# q1 I8 }5 g
) C& q  l9 e) Z7 |3 r3 R7 ^+ p1 O" d  s
) X4 ^: X8 C5 L- }6 s% {0 ?0 c: vfunction f = ObjFunc7LNL(k,x0,yexp)
tspan = [0 15 30 45 60 90 120 180 240 300 360];
[t, x] = ode45(@KineticEqs,tspan,x0,[],k);   
/ T' u2 y, j8 P! F4 c* vy(:,2) = x(:,1);
% W6 i1 n6 Y  `, _y(:,3:6) = x(:,2:5);
f1 = y(:,2) - yexp(:,2);
f2 = y(:,3) - yexp(:,3);
2 M$ ~" b; T8 E4 ]1 ~f3 = y(:,4) - yexp(:,4);
) Y! q. t, i! P) h# s! |# Hf4 = y(:,5) - yexp(:,5);
f5 = y(:,6) - yexp(:,6);
" d. x' r( M+ ~0 d$ nf = [f1; f2; f3; f4; f5];



# b0 V/ D4 W1 {8 X3 Kfunction f = ObjFunc7Fmincon(k,x0,yexp)
( H" D; B' n$ G3 ^tspan = [0 15 30 45 60 90 120 180 240 300 360];
, K& _5 {5 I0 q/ k[t x] = ode45(@KineticEqs,tspan,x0,[],k);   
, G  U0 U2 e% R1 m4 fy(:,2) = x(:,1);
- j1 `' _' C, _, ~. ky(:,3:6) = x(:,2:5);
1 T. ~; b: v8 Z, {9 w- U: yf =  sum((y(:,2)-yexp(:,2)).^2) + sum((y(:,3)-yexp(:,3)).^2)   ...
  M" r2 q/ ~# p# V5 D( _9 Z3 S    + sum((y(:,4)-yexp(:,4)).^2) + sum((y(:,5)-yexp(:,5)).^2)   ...
! Q) x  w$ ]  O4 W; _, V" o    + sum((y(:,6)-yexp(:,6)).^2) ;
7 B$ i) [8 i9 G9 a6 g, o& k: }

* O$ I) p& t- f( o7 f* h( X$ @1 q

function dxdt = KineticEqs(t,x,k)
/ q$ C' X5 @" u: }* @& I# TdGldt = 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);
9 y1 J- m% A6 U* a: @! W$ v) D# ^2 DdFadt = k(3)*x(1)+k(5)*x(2)+(k(6)+k(7))*x(5);
: y1 S( ~( O( rdLadt = k(7)*x(5);
dHmdt = k(4)*x(2)-(k(6)+k(7)+k(10))*x(5);
dxdt = [dGldt; dFrdt; dFadt; dLadt; dHmdt];


8 _1 A+ `$ ^; M$ N& t) o# z