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

箫剑→残念

function parafit
%  k1->k-1，k2->k1，k3->k2,k4->k3,k5->k4
% k6->k6 k7->k7
( ~) d" L0 z$ n% 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);
5 ]/ T2 n. J4 a+ ]% y% dLadt = k(7)*C(Hmf);
%dHmfdt = k(3)*C(Fru)-（k6+k7)*C(Hmf);
clear all
clc
  l7 t3 f% w" A9 P- Pformat long
% h- W7 D6 D# c7 ]! N3 n1 Q%        t/min   Glc    Fru        Fa   La   HMF/ mol/L
7 q9 z1 W, R, K8 ?* ~8 e3 K' K  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
          45    0.15105    0.06975    0.02473    0    0.0033
          60    0.13763    0.07397    0.02615    0    0.00428
. [- h0 K# g( `: F( d7 ]          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
  |/ _  K1 K9 W          240    0.03653    0.06488    0.08945    0.01828    0.05452
8 ~7 k/ [1 m# |4 @' o& Y          300    0.02738    0.05448    0.09098    0.0227    0.0597
          360    0.01855    0.04125    0.09363    0.0239    0.06495];
# d- d. q) i8 S7 L) a0 s+ uk0 = [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];                  % 参数下限
+ c  ?5 a0 d; j; ?( iub = [1  1  1  1  1  1  1  1  1  1];    % 参数上限
: o; t% ]0 v  f$ G. {( ]; Bx0 = [0.25  0  0  0  0];
yexp = Kinetics;                 % yexp: 实验数据[x1        x4        x5        x6]
/ e+ C3 ?# _2 w. {* _; w; n% warning off
% 使用函数 ()进行参数估计
6 `/ x1 S) V! `- X2 A% ]6 \+ Z! ^+ u[k,fval,flag] = fmincon(@ObjFunc7Fmincon,k0,[],[],[],[],lb,ub,[],[],x0,yexp);
fprintf('\n使用函数fmincon()估计得到的参数值为:\n')
fprintf('\tk1 = %.11f\n',k(1))
9 v! {2 T3 @; ^" a; Jfprintf('\tk2 = %.11f\n',k(2))
fprintf('\tk3 = %.11f\n',k(3))
- J  w2 g. D3 _1 efprintf('\tk4 = %.11f\n',k(4))
4 w4 b$ w9 a) zfprintf('\tk5 = %.11f\n',k(5))
fprintf('\tk6 = %.11f\n',k(6))
, Q# P+ G4 E1 |& `8 Vfprintf('\tk7 = %.11f\n',k(7))
fprintf('\tk8 = %.11f\n',k(8))
$ n. G2 `* q) b. yfprintf('\tk9 = %.11f\n',k(9))
fprintf('\tk10 = %.11f\n',k(10))
  j. Y+ f: k% y5 h# L: Rfprintf('  The sum of the squares is: %.1e\n\n',fval)
k_fm= k;
% warning off
% 使用函数lsqnonlin()进行参数估计
- j1 C( q) E& P. p  ?[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc7LNL,k0,lb,ub,[],x0,yexp);      
ci = nlparci(k,residual,jacobian);
/ M- P" O. f6 O& Ofprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\tk1 = %.11f\n',k(1))
fprintf('\tk2 = %.11f\n',k(2))
- R) z9 _7 q" t: H- Ifprintf('\tk3 = %.11f\n',k(3))
fprintf('\tk4 = %.11f\n',k(4))
fprintf('\tk5 = %.11f\n',k(5))
9 |9 c9 N; X. A1 k$ ]fprintf('\tk6 = %.11f\n',k(6))
fprintf('\tk7 = %.11f\n',k(7))
8 [" b- N& d# }0 ]' N1 |  A1 @fprintf('\tk8 = %.11f\n',k(8))
8 {) `1 \( B* \8 ^5 hfprintf('\tk9 = %.11f\n',k(9))
& A9 }8 D2 i/ W/ Ufprintf('\tk10 = %.11f\n',k(10))
fprintf('  The sum of the squares is: %.1e\n\n',resnorm)
7 ^& O6 s7 s% X& H* S6 Tk_ls = k;
output
/ y. x7 w0 R" H8 cwarning off
% 以函数fmincon()估计得到的结果为初值，使用函数lsqnonlin()进行参数估计
k0 = k_fm;
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc7LNL,k0,lb,ub,[],x0,yexp);      
! \3 K3 w/ K8 R8 K" Bci = nlparci(k,residual,jacobian);
2 m% b8 |7 S5 j  {2 n- `% F3 n+ kfprintf('\n\n以fmincon()的结果为初值，使用函数lsqnonlin()估计得到的参数值为:\n')
3 {- Q7 G" {% ^fprintf('\tk1 = %.11f\n',k(1))
fprintf('\tk2 = %.11f\n',k(2))
/ ~  \! |4 R  H: j6 s0 `# ]fprintf('\tk3 = %.11f\n',k(3))
$ X+ {- Q9 u9 _( h' }8 t" 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))
# N5 O3 S% r' ]0 a+ H6 r5 b. Rfprintf('  The sum of the squares is: %.1e\n\n',resnorm)
% g4 M  D& F" Ok_fmls = k;
' p% p: b! c6 a) H  `( B3 b6 P# [0 N5 |5 doutput
tspan = [0 15 30 45 60 90 120 180 240 300 360];
[t x] = ode45(@KineticEqs,tspan,x0,[],k_fmls);
figure;
5 b  u' [" {2 z( W' cplot(t,x(:,1),t,yexp(:,2),'*');legend('Glc-pr','Glc-real')
figure;plot(t,x(:,2:5));
" n( a6 x/ C" I& J# pp=x(:,1:5)
$ X! k7 N# Z" [5 W& N8 ~$ s) qhold on
plot(t,yexp(:,3:6),'o');legend('Fru-pr','Fa-pr','La-pr','HMF-pr','Fru-real','Fa-real','La-real','HMF-real')
5 e- e6 m: }) r4 e


function f = ObjFunc7LNL(k,x0,yexp)
1 P; S; l9 D. m9 F2 t' S1 s& r  G% Y: Ktspan = [0 15 30 45 60 90 120 180 240 300 360];
# N8 h& s# H; i$ q[t, x] = ode45(@KineticEqs,tspan,x0,[],k);   
" N) }" l3 s; ?' W0 i! u# Ry(:,2) = x(:,1);
, ~9 y/ l& I3 t0 }y(:,3:6) = x(:,2:5);
, ]+ z0 d$ M6 M5 k  Qf1 = y(:,2) - yexp(:,2);
) p8 }0 U& V$ N& J- Bf2 = y(:,3) - yexp(:,3);
/ t! R6 ]3 \3 J2 S# g1 Uf3 = y(:,4) - yexp(:,4);
7 b$ E4 j+ ^* q( [8 W) z" K8 Z7 pf4 = y(:,5) - yexp(:,5);
f5 = y(:,6) - yexp(:,6);
) i, p6 Y( ~% s/ Wf = [f1; f2; f3; f4; f5];
& L9 _( Z% z5 ]7 l
. W; u) a" Y! a4 _( d7 d" k

function f = ObjFunc7Fmincon(k,x0,yexp)
# ?- q& t& Z" D  [tspan = [0 15 30 45 60 90 120 180 240 300 360];
+ O2 J4 v* w+ m/ Z* J$ f9 c" g' K[t x] = ode45(@KineticEqs,tspan,x0,[],k);   
y(:,2) = x(:,1);
7 }/ G9 e6 O* T3 v# C  Ay(:,3:6) = x(:,2:5);
f =  sum((y(:,2)-yexp(:,2)).^2) + sum((y(:,3)-yexp(:,3)).^2)   ...
; ?3 t+ p+ i  {    + sum((y(:,4)-yexp(:,4)).^2) + sum((y(:,5)-yexp(:,5)).^2)   ...
    + sum((y(:,6)-yexp(:,6)).^2) ;
3 ^% i: F4 B$ W5 i" S

- @. q( D4 {' f. T, W/ m1 L
9 w0 t, v5 P; q' d5 n8 V
function dxdt = KineticEqs(t,x,k)
dGldt = k(1)*x(2)-(k(2)+k(3)+k(8))*x(1);
" v! c! J, a4 M% T- Z- ~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);
' O, R% i: ?% p; u5 y* x, DdLadt = k(7)*x(5);
0 F( a4 r- ^: D5 V0 y. p; p" z3 kdHmdt = k(4)*x(2)-(k(6)+k(7)+k(10))*x(5);
1 T; C6 H5 }0 D% p+ r( `" _# X1 Ldxdt = [dGldt; dFrdt; dFadt; dLadt; dHmdt];

" s: F( e" c, g# [' q$ L7 x
4 [! T; D0 F" k3 j8 U8 E! T( p

董事长之友

蒙的一比 大哥
) s7 O1 f/ v, z* e4 F; {