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

箫剑→残念

function parafit
%  k1->k-1，k2->k1，k3->k2,k4->k3,k5->k4
% 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);
) E, J  p% O( t0 j$ ~, B% dLadt = k(7)*C(Hmf);
' g& o+ }0 v8 V$ X# j%dHmfdt = k(3)*C(Fru)-（k6+k7)*C(Hmf);
; I) x7 M* |2 Q/ e! P) z! R: dclear all
' ^- D" g$ c" N. @clc
8 t! H  |0 Q9 h3 i, eformat long
+ s) k7 w4 `6 I, l1 i) q. B/ e%        t/min   Glc    Fru        Fa   La   HMF/ mol/L
5 q0 L. `1 W! p( p  Kinetics=[0    0.25    0           0    0       0
          15    0.2319    0.01257    0.0048    0    2.50E-04
          30    0.19345    0.027    0.00868    0    7.00E-04
% X7 B4 n9 V$ e' T( Q* Z6 L% X2 w, }          45    0.15105    0.06975    0.02473    0    0.0033
          60    0.13763    0.07397    0.02615    0    0.00428
/ }' C  C% K4 c- B) C3 |          90    0.08115    0.07877    0.07485    0    0.01405
          120    0.0656    0.07397    0.07885    0.00573    0.02143
. [  t, f3 r+ r  D" F1 A          180    0.04488    0.0682    0.07135    0.0091    0.03623
          240    0.03653    0.06488    0.08945    0.01828    0.05452
# P6 [4 E7 U" z3 ^! I- G4 E          300    0.02738    0.05448    0.09098    0.0227    0.0597
- v& m' G& H# Z8 K3 b- v          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];        % 参数初值
+ L; Q6 P. i6 M, x4 Mlb = [0  0  0  0  0  0  0  0  0  0];                  % 参数下限
ub = [1  1  1  1  1  1  1  1  1  1];    % 参数上限
0 a4 T& P5 e* r9 B8 t+ @$ U' i* `" t, }x0 = [0.25  0  0  0  0];
  l% u1 l% E# ^; P  f: W8 Q7 y# qyexp = Kinetics;                 % yexp: 实验数据[x1        x4        x5        x6]
' V1 \' s  o. Z' h; M, R% warning off
( x; u! ^+ p9 w0 N% 使用函数 ()进行参数估计
, Y  ]" \4 g3 i% s5 I* i' {[k,fval,flag] = fmincon(@ObjFunc7Fmincon,k0,[],[],[],[],lb,ub,[],[],x0,yexp);
fprintf('\n使用函数fmincon()估计得到的参数值为:\n')
5 u7 }6 D  V3 H( w! V" ]fprintf('\tk1 = %.11f\n',k(1))
/ i1 c& o! Y  N5 Z# Dfprintf('\tk2 = %.11f\n',k(2))
fprintf('\tk3 = %.11f\n',k(3))
fprintf('\tk4 = %.11f\n',k(4))
" @$ @4 {+ `$ @! r: tfprintf('\tk5 = %.11f\n',k(5))
) Q8 O  I1 G1 ], O. L0 {fprintf('\tk6 = %.11f\n',k(6))
5 R( A8 |& z: ~4 Q1 W+ M" |) D  ~5 L$ yfprintf('\tk7 = %.11f\n',k(7))
fprintf('\tk8 = %.11f\n',k(8))
$ x- b9 m/ t3 N* n5 [- N, y7 j- ofprintf('\tk9 = %.11f\n',k(9))
4 F* k( ]# v0 b, xfprintf('\tk10 = %.11f\n',k(10))
fprintf('  The sum of the squares is: %.1e\n\n',fval)
k_fm= k;
/ B+ i' f/ W- w8 N% warning off
% 使用函数lsqnonlin()进行参数估计
& }$ _2 S3 ]+ a  w; ~[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
0 a6 Q$ w5 R5 M& e# D    lsqnonlin(@ObjFunc7LNL,k0,lb,ub,[],x0,yexp);      
" ^! b- ]* T1 pci = nlparci(k,residual,jacobian);
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\tk1 = %.11f\n',k(1))
fprintf('\tk2 = %.11f\n',k(2))
* @2 A/ B, ], J& h  v$ U9 S' cfprintf('\tk3 = %.11f\n',k(3))
3 f# {: y  q; Ofprintf('\tk4 = %.11f\n',k(4))
' Q- U8 ?4 U4 C9 Zfprintf('\tk5 = %.11f\n',k(5))
fprintf('\tk6 = %.11f\n',k(6))
& U" O, K, L8 j: F- e: g( p- g: ~fprintf('\tk7 = %.11f\n',k(7))
  q2 v; I6 _8 g4 |! @fprintf('\tk8 = %.11f\n',k(8))
+ |- R% W9 V$ jfprintf('\tk9 = %.11f\n',k(9))
fprintf('\tk10 = %.11f\n',k(10))
. a( K( R$ P, G: p0 e  |: tfprintf('  The sum of the squares is: %.1e\n\n',resnorm)
8 E1 c, l5 l, d; j0 Ek_ls = k;
" S4 ?. r7 }! v# O( r) F" L! Soutput
# X. m; Q  r/ T6 c/ q* j/ K& V, a+ Zwarning off
, b1 q& `/ G  k$ I* v% 以函数fmincon()估计得到的结果为初值，使用函数lsqnonlin()进行参数估计
3 q6 I. o" \# I9 D8 gk0 = k_fm;
9 ]" ]( \4 b* _* h% f3 l3 N[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
  N) Z3 }8 {9 L2 a# G9 k    lsqnonlin(@ObjFunc7LNL,k0,lb,ub,[],x0,yexp);      
ci = nlparci(k,residual,jacobian);
fprintf('\n\n以fmincon()的结果为初值，使用函数lsqnonlin()估计得到的参数值为:\n')
/ [, r! e4 E5 d* Dfprintf('\tk1 = %.11f\n',k(1))
fprintf('\tk2 = %.11f\n',k(2))
7 J7 w" N( O4 M. V7 q) Cfprintf('\tk3 = %.11f\n',k(3))
" M4 Y; E6 i4 yfprintf('\tk4 = %.11f\n',k(4))
fprintf('\tk5 = %.11f\n',k(5))
fprintf('\tk6 = %.11f\n',k(6))
6 D! L- l4 r4 U1 Y: J' \( bfprintf('\tk7 = %.11f\n',k(7))
fprintf('\tk8 = %.11f\n',k(8))
fprintf('\tk9 = %.11f\n',k(9))
9 n; C$ i8 J+ ~1 `fprintf('\tk10 = %.11f\n',k(10))
fprintf('  The sum of the squares is: %.1e\n\n',resnorm)
  f8 g0 D) U$ D' G' l$ `6 C* Xk_fmls = k;
+ [; f, u" w+ @, eoutput
tspan = [0 15 30 45 60 90 120 180 240 300 360];
+ q0 Y5 G( y: x2 X" T8 e0 Y. ~[t x] = ode45(@KineticEqs,tspan,x0,[],k_fmls);
9 q# ~8 J: Q8 m2 A: C- Bfigure;
4 i: j  B; T$ ~' Y8 Uplot(t,x(:,1),t,yexp(:,2),'*');legend('Glc-pr','Glc-real')
1 G1 M# f* e; ~3 q. g. Tfigure;plot(t,x(:,2:5));
0 Y1 \9 e$ K* c- U: z6 l: M+ _p=x(:,1:5)
' g9 Z8 K5 Z' B$ khold on
plot(t,yexp(:,3:6),'o');legend('Fru-pr','Fa-pr','La-pr','HMF-pr','Fru-real','Fa-real','La-real','HMF-real')
/ W" o' E( n6 X" @


function f = ObjFunc7LNL(k,x0,yexp)
tspan = [0 15 30 45 60 90 120 180 240 300 360];
[t, x] = ode45(@KineticEqs,tspan,x0,[],k);   
y(:,2) = x(:,1);
4 \6 R6 ^) G: ~( z( e) Ry(:,3:6) = x(:,2:5);
f1 = y(:,2) - yexp(:,2);
- [0 v+ W0 X, v  r) m: h0 |& c. s0 ~f2 = y(:,3) - yexp(:,3);
f3 = y(:,4) - yexp(:,4);
/ w- N- Z* U# |) }. l& l* w2 \" @f4 = y(:,5) - yexp(:,5);
/ m' E- f4 V" F4 I8 U) E9 r. L7 z" pf5 = y(:,6) - yexp(:,6);
" }* v: ]7 C; }0 d* Z1 n8 \8 Df = [f1; f2; f3; f4; f5];

# b+ E) i8 J, u1 A$ `8 [

function f = ObjFunc7Fmincon(k,x0,yexp)
4 [! w" X- y/ M8 \tspan = [0 15 30 45 60 90 120 180 240 300 360];
7 N# W; m( p. E0 D( w. O* w5 @[t x] = ode45(@KineticEqs,tspan,x0,[],k);   
y(:,2) = x(:,1);
' i2 o; ~2 t( m) m7 ^! ey(:,3:6) = x(:,2:5);
( ]& E# N0 Z+ i* b! `. R; wf =  sum((y(:,2)-yexp(:,2)).^2) + sum((y(:,3)-yexp(:,3)).^2)   ...
! d* G# C, i2 g4 q    + sum((y(:,4)-yexp(:,4)).^2) + sum((y(:,5)-yexp(:,5)).^2)   ...
    + sum((y(:,6)-yexp(:,6)).^2) ;

" I( A9 ]! |; F4 I" }5 o
& S0 h& H6 G: G( b

function dxdt = KineticEqs(t,x,k)
dGldt = k(1)*x(2)-(k(2)+k(3)+k(8))*x(1);
8 h2 l# X. i5 W, ]- L$ M- zdFrdt = k(2)*x(1)-(k(1)+k(4)+k(5)+k(9))*x(2);
) S4 p6 N( h7 V# [: O& ?4 k) JdFadt = 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);
( E- d- T. D( x0 \dxdt = [dGldt; dFrdt; dFadt; dLadt; dHmdt];
5 y+ o4 j1 J1 y& R3 u  s

# i. a% _0 d: {, ~$ ~! p