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

箫剑→残念

function parafit
%  k1->k-1，k2->k1，k3->k2,k4->k3,k5->k4
5 X1 N' p  \+ v' q5 {2 g8 c" C! u) T% k6->k6 k7->k7
. U+ G2 f8 R  l* ]. K1 j% dGlcdt = k-1*C(Fru)-(k1+k2)*C(Glc);
' r/ v; q% {: O# ?4 @) ^% dFrudt = k1*C(Glc)-(k-1+k3+k4)C(Fru);
1 E' g8 c1 D; g; X5 `% A2 J% 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
4 }; b% v- l; q9 J. D, h% A/ Mclc
  v* x9 f3 ?8 _- aformat 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
          30    0.19345    0.027    0.00868    0    7.00E-04
" v0 ]( M( _8 u5 n0 V5 x- t8 [          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
          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
          360    0.01855    0.04125    0.09363    0.0239    0.06495];
: }8 d+ r$ n1 h. ]$ jk0 = [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];                  % 参数下限
+ B: G5 |# {8 k0 B) [ub = [1  1  1  1  1  1  1  1  1  1];    % 参数上限
5 f, u1 _" N: a; _x0 = [0.25  0  0  0  0];
) c6 |4 [8 B- T+ @) x6 a1 Iyexp = Kinetics;                 % yexp: 实验数据[x1        x4        x5        x6]
% warning off
% 使用函数 ()进行参数估计
[k,fval,flag] = fmincon(@ObjFunc7Fmincon,k0,[],[],[],[],lb,ub,[],[],x0,yexp);
4 H& F8 b  ?( t9 `! G# d( Afprintf('\n使用函数fmincon()估计得到的参数值为:\n')
fprintf('\tk1 = %.11f\n',k(1))
0 {, ~$ o1 w7 w% n+ \fprintf('\tk2 = %.11f\n',k(2))
fprintf('\tk3 = %.11f\n',k(3))
8 P% u# ?- R& j# _( K4 p2 @fprintf('\tk4 = %.11f\n',k(4))
fprintf('\tk5 = %.11f\n',k(5))
& G5 l/ P7 C# f  Yfprintf('\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))
1 r- S! B1 V+ Pfprintf('  The sum of the squares is: %.1e\n\n',fval)
" P2 w$ }: N: c) F* O: _k_fm= k;
% warning off
  f, i" m3 j: W7 z% 使用函数lsqnonlin()进行参数估计
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc7LNL,k0,lb,ub,[],x0,yexp);      
ci = nlparci(k,residual,jacobian);
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\tk1 = %.11f\n',k(1))
+ r" p! n4 D; B( {& ?3 rfprintf('\tk2 = %.11f\n',k(2))
' }  m6 r3 h, f& h1 t/ d. Cfprintf('\tk3 = %.11f\n',k(3))
fprintf('\tk4 = %.11f\n',k(4))
fprintf('\tk5 = %.11f\n',k(5))
% L3 X3 D% L# wfprintf('\tk6 = %.11f\n',k(6))
* H$ J3 Y) h) a+ J" s: R0 P2 Mfprintf('\tk7 = %.11f\n',k(7))
fprintf('\tk8 = %.11f\n',k(8))
0 w( ]) e/ y6 x# |0 v, E1 I# c0 ^/ Nfprintf('\tk9 = %.11f\n',k(9))
5 ~+ y( ^! @: ?- R- G, e+ Zfprintf('\tk10 = %.11f\n',k(10))
fprintf('  The sum of the squares is: %.1e\n\n',resnorm)
k_ls = k;
output
! m4 ~* o8 i5 M5 bwarning off
% 以函数fmincon()估计得到的结果为初值，使用函数lsqnonlin()进行参数估计
k0 = k_fm;
6 p6 g. q( N6 ^- `[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc7LNL,k0,lb,ub,[],x0,yexp);      
ci = nlparci(k,residual,jacobian);
. |9 ?! x' y! Z/ u# Yfprintf('\n\n以fmincon()的结果为初值，使用函数lsqnonlin()估计得到的参数值为:\n')
1 Y8 B; B, Z2 |. \2 nfprintf('\tk1 = %.11f\n',k(1))
3 w% L6 l$ H+ Pfprintf('\tk2 = %.11f\n',k(2))
fprintf('\tk3 = %.11f\n',k(3))
fprintf('\tk4 = %.11f\n',k(4))
. B: q* e7 i' \" j( |9 kfprintf('\tk5 = %.11f\n',k(5))
fprintf('\tk6 = %.11f\n',k(6))
8 Z3 m7 ?& e4 f7 _- \+ Kfprintf('\tk7 = %.11f\n',k(7))
, M. h  b: R' o' f+ h$ [fprintf('\tk8 = %.11f\n',k(8))
' h5 ~; ?. N7 G/ e5 jfprintf('\tk9 = %.11f\n',k(9))
1 v9 X9 Q7 O: k2 Z3 `- Jfprintf('\tk10 = %.11f\n',k(10))
( z& Y. c- }1 }7 {1 l1 K6 X! m: ~fprintf('  The sum of the squares is: %.1e\n\n',resnorm)
k_fmls = k;
, f3 w8 d8 w" f- z8 B0 T" \6 boutput
7 A4 ]+ {3 }& B; O$ h5 S0 Ktspan = [0 15 30 45 60 90 120 180 240 300 360];
[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)
. l$ k5 T+ E2 `0 \$ a3 Q% R/ Ghold on
+ N7 a: A8 ?; ^4 f) W7 Gplot(t,yexp(:,3:6),'o');legend('Fru-pr','Fa-pr','La-pr','HMF-pr','Fru-real','Fa-real','La-real','HMF-real')


5 p( _3 F# k7 l9 o( l
function f = ObjFunc7LNL(k,x0,yexp)
1 t% d4 K3 O- j, Wtspan = [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);
, `4 r9 S' f6 Y# cf1 = y(:,2) - yexp(:,2);
- W# W7 M: F8 Gf2 = y(:,3) - yexp(:,3);
f3 = y(:,4) - yexp(:,4);
f4 = y(:,5) - yexp(:,5);
1 h: U' d2 }  Y9 E0 n6 yf5 = y(:,6) - yexp(:,6);
f = [f1; f2; f3; f4; f5];
% P/ N! f. J! b# n/ l' y


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);   
- i1 m% r1 |4 E: y+ D6 l# Sy(:,2) = x(:,1);
y(:,3:6) = x(:,2:5);
) c8 q9 G9 R0 x/ }/ jf =  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) ;

0 J- @. L3 t# t% p+ b7 e& Y. ^
2 V1 y" n3 r+ V" r! Z; E

1 [& w) ^+ y7 Xfunction 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);
, c: o4 X* @" n4 A8 N( c; W) fdFadt = k(3)*x(1)+k(5)*x(2)+(k(6)+k(7))*x(5);
dLadt = k(7)*x(5);
: V% x3 @6 ^. X. \) [4 h4 [1 tdHmdt = k(4)*x(2)-(k(6)+k(7)+k(10))*x(5);
7 F, N; Z. d1 U6 a3 N! L; c1 @5 Wdxdt = [dGldt; dFrdt; dFadt; dLadt; dHmdt];

" k+ Z3 {8 d& W. z- b/ @" w$ S