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

箫剑→残念

function parafit
%  k1->k-1，k2->k1，k3->k2,k4->k3,k5->k4
$ b7 V2 F1 ]8 R8 @  }8 ~9 P% g% k6->k6 k7->k7
' S; X! E! a, J9 y% v% 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);
%dHmfdt = k(3)*C(Fru)-（k6+k7)*C(Hmf);
clear all
clc
6 G% B: [3 `, O# t5 N. pformat long
$ {* w# D$ B, w$ Z9 S%        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
8 R: x' s8 u5 _" w% L          30    0.19345    0.027    0.00868    0    7.00E-04
          45    0.15105    0.06975    0.02473    0    0.0033
          60    0.13763    0.07397    0.02615    0    0.00428
( [: Y2 h$ h' A- }          90    0.08115    0.07877    0.07485    0    0.01405
( ^7 l+ ?+ U* m' T4 i3 t8 j! D          120    0.0656    0.07397    0.07885    0.00573    0.02143
          180    0.04488    0.0682    0.07135    0.0091    0.03623
5 e4 \, B& e/ e/ f7 |          240    0.03653    0.06488    0.08945    0.01828    0.05452
) `  y4 D7 r. r, ?3 y' ^; X          300    0.02738    0.05448    0.09098    0.0227    0.0597
          360    0.01855    0.04125    0.09363    0.0239    0.06495];
+ |- @5 r: @- c. yk0 = [0.0000000005  0.0000000005  0.0000000005  0.00000000005  0.00005  0.0134  0.00564  0.00001  0.00001  0.00001];        % 参数初值
# a/ q7 T9 y9 U. Jlb = [0  0  0  0  0  0  0  0  0  0];                  % 参数下限
ub = [1  1  1  1  1  1  1  1  1  1];    % 参数上限
7 U% z8 `" x! mx0 = [0.25  0  0  0  0];
' A- x" K" F  i# ?; k. T# hyexp = Kinetics;                 % yexp: 实验数据[x1        x4        x5        x6]
  @5 i. h( R' Q! {/ e+ r% warning off
% 使用函数 ()进行参数估计
) S9 ]7 Z/ s% t: r# g# P[k,fval,flag] = fmincon(@ObjFunc7Fmincon,k0,[],[],[],[],lb,ub,[],[],x0,yexp);
fprintf('\n使用函数fmincon()估计得到的参数值为:\n')
5 s: ?. n3 l# p& ?8 nfprintf('\tk1 = %.11f\n',k(1))
fprintf('\tk2 = %.11f\n',k(2))
fprintf('\tk3 = %.11f\n',k(3))
/ I) l6 ?4 I) D# i$ efprintf('\tk4 = %.11f\n',k(4))
fprintf('\tk5 = %.11f\n',k(5))
fprintf('\tk6 = %.11f\n',k(6))
  g; ]1 o( K- @4 Jfprintf('\tk7 = %.11f\n',k(7))
  E8 s( {; b$ [6 ?: R3 L8 j! @  vfprintf('\tk8 = %.11f\n',k(8))
fprintf('\tk9 = %.11f\n',k(9))
fprintf('\tk10 = %.11f\n',k(10))
7 {/ Q( D" P! q8 gfprintf('  The sum of the squares is: %.1e\n\n',fval)
k_fm= k;
% warning off
% 使用函数lsqnonlin()进行参数估计
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
; Y4 z5 p: Z( O    lsqnonlin(@ObjFunc7LNL,k0,lb,ub,[],x0,yexp);      
' `* ?" l& ]: l+ fci = nlparci(k,residual,jacobian);
2 H- d' U% b  efprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\tk1 = %.11f\n',k(1))
fprintf('\tk2 = %.11f\n',k(2))
fprintf('\tk3 = %.11f\n',k(3))
fprintf('\tk4 = %.11f\n',k(4))
6 P9 m; V% O* Y% [$ Bfprintf('\tk5 = %.11f\n',k(5))
" d- v. h1 r0 i4 u8 Ffprintf('\tk6 = %.11f\n',k(6))
' L0 x* y! l, s- Tfprintf('\tk7 = %.11f\n',k(7))
9 M, p1 u7 O. wfprintf('\tk8 = %.11f\n',k(8))
fprintf('\tk9 = %.11f\n',k(9))
fprintf('\tk10 = %.11f\n',k(10))
2 \' j( {* V& d' g5 I. `% Kfprintf('  The sum of the squares is: %.1e\n\n',resnorm)
k_ls = k;
output
2 u: M: q# ~) @/ Z( owarning off
% 以函数fmincon()估计得到的结果为初值，使用函数lsqnonlin()进行参数估计
1 C) U7 C- X5 j8 s6 hk0 = k_fm;
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
* C8 a0 J7 y$ l0 Q6 A    lsqnonlin(@ObjFunc7LNL,k0,lb,ub,[],x0,yexp);      
ci = nlparci(k,residual,jacobian);
; i' P  b( G" \, x. Gfprintf('\n\n以fmincon()的结果为初值，使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\tk1 = %.11f\n',k(1))
+ n! W, L- G1 H- i. ^fprintf('\tk2 = %.11f\n',k(2))
fprintf('\tk3 = %.11f\n',k(3))
0 M$ X$ r8 q' y" {7 Nfprintf('\tk4 = %.11f\n',k(4))
8 c8 b3 G& D3 T  z1 g  n+ Jfprintf('\tk5 = %.11f\n',k(5))
) F- k! `0 E( D0 ?# p$ f3 Lfprintf('\tk6 = %.11f\n',k(6))
fprintf('\tk7 = %.11f\n',k(7))
fprintf('\tk8 = %.11f\n',k(8))
3 x' p9 v0 w7 t  a2 i- ]: g4 Tfprintf('\tk9 = %.11f\n',k(9))
% S- y6 n% w; bfprintf('\tk10 = %.11f\n',k(10))
fprintf('  The sum of the squares is: %.1e\n\n',resnorm)
. V# ]8 I& D+ h+ P. ]k_fmls = k;
output
7 }# G% y9 p1 x* o( Q/ k6 ]tspan = [0 15 30 45 60 90 120 180 240 300 360];
% v, ^: _6 ]' y& N+ a0 k- C- l[t x] = ode45(@KineticEqs,tspan,x0,[],k_fmls);
% @& v: e3 x" Q6 d. C! w' Z7 Pfigure;
/ ^$ p" S+ v5 c0 `% y; gplot(t,x(:,1),t,yexp(:,2),'*');legend('Glc-pr','Glc-real')
4 s' A$ n9 [# {3 ~$ f, |" |figure;plot(t,x(:,2:5));
3 y  x/ x) t, a0 w+ xp=x(:,1:5)
/ A* W1 m1 b; m7 l9 h' {* E( vhold on
9 j2 _. E! u, s' x% H4 Iplot(t,yexp(:,3:6),'o');legend('Fru-pr','Fa-pr','La-pr','HMF-pr','Fru-real','Fa-real','La-real','HMF-real')
. s9 ?( w# X! f4 a2 H- h- i- Z4 Z) Y# \
( ]$ w" q2 e& w" T

( t6 s3 r$ u8 b# ?function f = ObjFunc7LNL(k,x0,yexp)
6 k' ^5 P& G. Ltspan = [0 15 30 45 60 90 120 180 240 300 360];
[t, x] = ode45(@KineticEqs,tspan,x0,[],k);   
! Z; k; P0 i: ~7 D0 |& ?y(:,2) = x(:,1);
- Y  j3 [# E# S' A5 I; V! g2 ?y(:,3:6) = x(:,2:5);
- E" O8 I6 `! V' j5 `' e5 G  T6 Tf1 = y(:,2) - yexp(:,2);
) h6 f' f8 ~8 C5 J0 \3 Z* Yf2 = y(:,3) - yexp(:,3);
f3 = y(:,4) - yexp(:,4);
/ ^5 `: k$ {% h' Z% `; tf4 = y(:,5) - yexp(:,5);
f5 = y(:,6) - yexp(:,6);
  j- ?  L6 Z+ P4 h# A5 Qf = [f1; f2; f3; f4; f5];

$ L6 ~) [; I. q: M* S6 q

function f = ObjFunc7Fmincon(k,x0,yexp)
tspan = [0 15 30 45 60 90 120 180 240 300 360];
[t x] = ode45(@KineticEqs,tspan,x0,[],k);   
( U. ^3 X5 p8 v1 ?6 Iy(:,2) = x(:,1);
y(:,3:6) = x(:,2:5);
. x1 `9 m0 R+ O  ef =  sum((y(:,2)-yexp(:,2)).^2) + sum((y(:,3)-yexp(:,3)).^2)   ...
    + sum((y(:,4)-yexp(:,4)).^2) + sum((y(:,5)-yexp(:,5)).^2)   ...
    + sum((y(:,6)-yexp(:,6)).^2) ;

6 X9 N2 e! n7 |# Y, A$ F$ b
5 d/ w  z8 K! ]

5 U8 H7 K) O6 C2 sfunction dxdt = KineticEqs(t,x,k)
dGldt = k(1)*x(2)-(k(2)+k(3)+k(8))*x(1);
' R! _* |& I, a. ~2 [3 A8 wdFrdt = 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);
- @! g) m* g3 b" |- o) EdLadt = k(7)*x(5);
dHmdt = k(4)*x(2)-(k(6)+k(7)+k(10))*x(5);
) x3 \! _( d0 x* kdxdt = [dGldt; dFrdt; dFadt; dLadt; dHmdt];