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

箫剑→残念

function parafit
3 y4 V! x7 V7 T5 w# d%  k1->k-1，k2->k1，k3->k2,k4->k3,k5->k4
% k6->k6 k7->k7
% dGlcdt = k-1*C(Fru)-(k1+k2)*C(Glc);
# d! R& {& K6 D) l* h% dFrudt = k1*C(Glc)-(k-1+k3+k4)C(Fru);
% dFadt = k(2)*C(Glc)+k4*C(Fru)+（k6+k7)*C(Hmf);
* X! J" t) Z; D' Y$ K% dLadt = k(7)*C(Hmf);
%dHmfdt = k(3)*C(Fru)-（k6+k7)*C(Hmf);
4 ^6 B0 ?) A" A- y! A5 S# Aclear all
. f; D& \3 i" W  Q# Mclc
% l( Z1 T* J, l, A; sformat long
8 K' I: l8 l- R: j$ Z) l%        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
$ M6 h5 m) Y3 G8 v7 H1 y          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
* l6 G$ Z6 b4 X( P( g          90    0.08115    0.07877    0.07485    0    0.01405
: ~% M4 v% \' G* l          120    0.0656    0.07397    0.07885    0.00573    0.02143
          180    0.04488    0.0682    0.07135    0.0091    0.03623
" P, f$ J' S* S" O          240    0.03653    0.06488    0.08945    0.01828    0.05452
          300    0.02738    0.05448    0.09098    0.0227    0.0597
: _' n; W3 o' B$ r          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];                  % 参数下限
8 P8 F: v) k/ i' G9 gub = [1  1  1  1  1  1  1  1  1  1];    % 参数上限
x0 = [0.25  0  0  0  0];
yexp = Kinetics;                 % yexp: 实验数据[x1        x4        x5        x6]
* U7 \6 b8 f5 X) w' h) g% warning off
% 使用函数 ()进行参数估计
[k,fval,flag] = fmincon(@ObjFunc7Fmincon,k0,[],[],[],[],lb,ub,[],[],x0,yexp);
fprintf('\n使用函数fmincon()估计得到的参数值为:\n')
fprintf('\tk1 = %.11f\n',k(1))
fprintf('\tk2 = %.11f\n',k(2))
fprintf('\tk3 = %.11f\n',k(3))
. n$ p2 R0 G# s+ ]* ]$ a( ufprintf('\tk4 = %.11f\n',k(4))
fprintf('\tk5 = %.11f\n',k(5))
fprintf('\tk6 = %.11f\n',k(6))
# k; g! g7 g7 T' i- E( Nfprintf('\tk7 = %.11f\n',k(7))
9 h# G" r) [- T9 Hfprintf('\tk8 = %.11f\n',k(8))
. ^9 k! ^% |4 ^9 kfprintf('\tk9 = %.11f\n',k(9))
fprintf('\tk10 = %.11f\n',k(10))
2 m! s9 D; s( G# J$ Y( Hfprintf('  The sum of the squares is: %.1e\n\n',fval)
k_fm= k;
% warning off
% 使用函数lsqnonlin()进行参数估计
2 M2 V& ]7 ~, a. S% Z- S[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))
/ L8 t7 P0 n! N4 a# `0 V; Mfprintf('\tk2 = %.11f\n',k(2))
fprintf('\tk3 = %.11f\n',k(3))
fprintf('\tk4 = %.11f\n',k(4))
. J# w0 l; l! K  k/ ]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))
0 W1 \# b4 i. I5 b0 V5 e0 Wfprintf('  The sum of the squares is: %.1e\n\n',resnorm)
0 n0 F1 {3 `$ \' F, ?k_ls = k;
& C4 D! ~1 U2 H" ~: ooutput
warning off
1 M% r! P# t# ]" ~+ d0 s1 y% 以函数fmincon()估计得到的结果为初值，使用函数lsqnonlin()进行参数估计
k0 = k_fm;
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
) {2 J- X1 \( B6 L) y' o    lsqnonlin(@ObjFunc7LNL,k0,lb,ub,[],x0,yexp);      
ci = nlparci(k,residual,jacobian);
; R! k! O) o9 X1 u  y: dfprintf('\n\n以fmincon()的结果为初值，使用函数lsqnonlin()估计得到的参数值为:\n')
3 R. @5 o+ T" |+ c/ t$ dfprintf('\tk1 = %.11f\n',k(1))
fprintf('\tk2 = %.11f\n',k(2))
fprintf('\tk3 = %.11f\n',k(3))
fprintf('\tk4 = %.11f\n',k(4))
) }5 L7 y/ q2 {1 v& n- rfprintf('\tk5 = %.11f\n',k(5))
  J( O2 H# ], j5 ^8 A% f- ?fprintf('\tk6 = %.11f\n',k(6))
fprintf('\tk7 = %.11f\n',k(7))
fprintf('\tk8 = %.11f\n',k(8))
! w: |0 I) b1 i6 Z( q0 ifprintf('\tk9 = %.11f\n',k(9))
fprintf('\tk10 = %.11f\n',k(10))
8 t+ ]) A7 c8 Lfprintf('  The sum of the squares is: %.1e\n\n',resnorm)
6 G$ g$ H/ b- \* b* d, \, u1 S) |k_fmls = k;
4 J+ l! p& O1 M1 H& y5 D& N9 Z; Soutput
tspan = [0 15 30 45 60 90 120 180 240 300 360];
8 j9 Z# I( i% t; V[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)
& o; h2 s2 D0 l7 w# }; }8 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')
* i6 ]) I+ E. ?* S5 r/ n

; }( ~2 x$ X1 H2 K
function f = ObjFunc7LNL(k,x0,yexp)
tspan = [0 15 30 45 60 90 120 180 240 300 360];
, H- B4 J8 i+ c' Q[t, x] = ode45(@KineticEqs,tspan,x0,[],k);   
y(:,2) = x(:,1);
) X9 [( j7 a3 ?% `1 Q. L* i5 x! Ay(:,3:6) = x(:,2:5);
- B, T. G& U; n6 W# Xf1 = y(:,2) - yexp(:,2);
f2 = y(:,3) - yexp(:,3);
: `+ D5 P" a) R. ^f3 = y(:,4) - yexp(:,4);
$ a8 j  w6 r& U) d: If4 = y(:,5) - yexp(:,5);
- J' N: ~- C7 Q! T* lf5 = y(:,6) - yexp(:,6);
f = [f1; f2; f3; f4; f5];
% C! `# P  b* B: _1 \' y
- e2 M! @8 M& u1 p8 N( I
" G7 F8 ?" i: l- s: z
  G- z% V  L9 Afunction f = ObjFunc7Fmincon(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);
, ?$ u( P; ?" H$ Q. k/ ?& Q" Jy(:,3:6) = x(:,2:5);
2 ?! D3 N% A$ U9 X: nf =  sum((y(:,2)-yexp(:,2)).^2) + sum((y(:,3)-yexp(:,3)).^2)   ...
    + sum((y(:,4)-yexp(:,4)).^2) + sum((y(:,5)-yexp(:,5)).^2)   ...
; \8 L/ n2 ]+ n9 Z    + sum((y(:,6)-yexp(:,6)).^2) ;
4 c3 x: G, i' A# i( |

5 F( m2 R0 n" e3 F  `, O8 b

function dxdt = KineticEqs(t,x,k)
8 L, A7 h% e6 U3 ^* ]$ w1 @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);
: w) [: P; p6 L1 ?% E8 d" GdFadt = k(3)*x(1)+k(5)*x(2)+(k(6)+k(7))*x(5);
" b. K4 P; [( ]- o" n( mdLadt = k(7)*x(5);
% K; z3 N7 ^; \# b0 ydHmdt = k(4)*x(2)-(k(6)+k(7)+k(10))*x(5);
dxdt = [dGldt; dFrdt; dFadt; dLadt; dHmdt];

' j( Y5 @( D0 g1 d1 b* v
) `; n% u* d2 ]: b* Q" P8 Z