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

箫剑→残念

function parafit
%  k1->k-1，k2->k1，k3->k2,k4->k3,k5->k4
  i6 a" h: P7 q0 L1 ?1 j9 V9 d# `) Y% k6->k6 k7->k7
% dGlcdt = k-1*C(Fru)-(k1+k2)*C(Glc);
# S4 I! ^, a  `6 C% dFrudt = k1*C(Glc)-(k-1+k3+k4)C(Fru);
0 v& y  q! X/ u( g! {3 G- k, r! p% dFadt = k(2)*C(Glc)+k4*C(Fru)+（k6+k7)*C(Hmf);
( i% j  x( d2 C! X% dLadt = k(7)*C(Hmf);
%dHmfdt = k(3)*C(Fru)-（k6+k7)*C(Hmf);
clear all
' G+ b. W! ^7 f( {$ Xclc
format long
%        t/min   Glc    Fru        Fa   La   HMF/ mol/L
5 d+ j7 \! c2 Y/ {; p3 U  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
2 i( y; t) [* n# s4 s. @          45    0.15105    0.06975    0.02473    0    0.0033
9 v2 B0 o) k$ w( P: R9 x4 W! ?          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
6 f( d& D7 {1 ^7 S6 }" l" N8 d& E          180    0.04488    0.0682    0.07135    0.0091    0.03623
          240    0.03653    0.06488    0.08945    0.01828    0.05452
# @1 N; L6 M4 A  @1 [4 j' R4 o          300    0.02738    0.05448    0.09098    0.0227    0.0597
7 m) p" b3 h" p5 |/ X          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];        % 参数初值
5 V7 O3 y% }5 N: Z$ W9 {lb = [0  0  0  0  0  0  0  0  0  0];                  % 参数下限
. Z! q5 I' K5 f6 D5 cub = [1  1  1  1  1  1  1  1  1  1];    % 参数上限
: b  x; d+ O/ m9 y! [x0 = [0.25  0  0  0  0];
yexp = Kinetics;                 % yexp: 实验数据[x1        x4        x5        x6]
7 j1 b3 B8 I" j, ?* a  w9 W% warning off
% 使用函数 ()进行参数估计
[k,fval,flag] = fmincon(@ObjFunc7Fmincon,k0,[],[],[],[],lb,ub,[],[],x0,yexp);
fprintf('\n使用函数fmincon()估计得到的参数值为:\n')
3 @& G& D6 u( T. I6 h3 |1 P7 Mfprintf('\tk1 = %.11f\n',k(1))
2 a9 M# c8 V+ R1 y1 {fprintf('\tk2 = %.11f\n',k(2))
fprintf('\tk3 = %.11f\n',k(3))
fprintf('\tk4 = %.11f\n',k(4))
fprintf('\tk5 = %.11f\n',k(5))
1 F! _* u* c1 h( O9 X) [fprintf('\tk6 = %.11f\n',k(6))
fprintf('\tk7 = %.11f\n',k(7))
fprintf('\tk8 = %.11f\n',k(8))
' R; ?: \8 T4 P2 U9 w6 ^fprintf('\tk9 = %.11f\n',k(9))
fprintf('\tk10 = %.11f\n',k(10))
+ l/ R8 n4 S4 g, @* F8 Efprintf('  The sum of the squares is: %.1e\n\n',fval)
k_fm= k;
, y+ y6 L. ^7 x% warning off
% 使用函数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))
4 R$ `! N# D' C/ T  e* pfprintf('\tk2 = %.11f\n',k(2))
fprintf('\tk3 = %.11f\n',k(3))
9 }# {4 a! v& w' K% ?fprintf('\tk4 = %.11f\n',k(4))
+ X0 B- }2 O3 Tfprintf('\tk5 = %.11f\n',k(5))
% m& B. _, E. |5 B' w2 \( q# efprintf('\tk6 = %.11f\n',k(6))
: m' \# [5 ]9 {8 p* ~/ O- @fprintf('\tk7 = %.11f\n',k(7))
5 z& I4 l& C% F# ?fprintf('\tk8 = %.11f\n',k(8))
. q2 ^/ ]! k2 e3 b4 t: bfprintf('\tk9 = %.11f\n',k(9))
fprintf('\tk10 = %.11f\n',k(10))
' P8 E) e+ V7 X0 F2 tfprintf('  The sum of the squares is: %.1e\n\n',resnorm)
k_ls = k;
' ?$ ^8 K" S  j, u7 routput
( U+ ]+ t% D2 d4 d( jwarning off
% 以函数fmincon()估计得到的结果为初值，使用函数lsqnonlin()进行参数估计
+ w% w: R5 S$ ^1 p- |4 g/ r$ u5 {" }k0 = k_fm;
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
5 Q6 z& G8 ?6 o    lsqnonlin(@ObjFunc7LNL,k0,lb,ub,[],x0,yexp);      
/ J( o$ Z5 Z. Y. \% k" sci = nlparci(k,residual,jacobian);
fprintf('\n\n以fmincon()的结果为初值，使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\tk1 = %.11f\n',k(1))
fprintf('\tk2 = %.11f\n',k(2))
% ?! [4 D% `; [0 X/ E: ifprintf('\tk3 = %.11f\n',k(3))
fprintf('\tk4 = %.11f\n',k(4))
fprintf('\tk5 = %.11f\n',k(5))
fprintf('\tk6 = %.11f\n',k(6))
- H$ y" _6 |) d6 p: [3 Afprintf('\tk7 = %.11f\n',k(7))
, f7 _, R9 `4 S/ V# s  g. H6 E; r$ Hfprintf('\tk8 = %.11f\n',k(8))
fprintf('\tk9 = %.11f\n',k(9))
fprintf('\tk10 = %.11f\n',k(10))
5 y) F( ~5 f, l$ [. ~4 Kfprintf('  The sum of the squares is: %.1e\n\n',resnorm)
/ S& ]% o! r& V6 k0 Rk_fmls = k;
* T/ J# T  K# ~6 ?output
1 r  h* }5 }$ Ftspan = [0 15 30 45 60 90 120 180 240 300 360];
[t x] = ode45(@KineticEqs,tspan,x0,[],k_fmls);
  D( h! D" H8 W3 ?, G2 h( A7 yfigure;
& t( W$ ]* V# M( I; J9 v/ Xplot(t,x(:,1),t,yexp(:,2),'*');legend('Glc-pr','Glc-real')
9 @& _5 c: m. @3 D+ W4 Q$ @$ [figure;plot(t,x(:,2:5));
/ ~' f- n" ?; w8 g+ ]; F4 x: p6 g! np=x(:,1:5)
hold on
plot(t,yexp(:,3:6),'o');legend('Fru-pr','Fa-pr','La-pr','HMF-pr','Fru-real','Fa-real','La-real','HMF-real')
8 O! r4 e) H5 p5 S. `
$ _9 A1 c* ~& W6 L5 B  U
5 Z+ x" o+ \' b+ N* `: M2 ~
function f = ObjFunc7LNL(k,x0,yexp)
tspan = [0 15 30 45 60 90 120 180 240 300 360];
$ P' l6 Y1 r& W" \4 O* \# V[t, x] = ode45(@KineticEqs,tspan,x0,[],k);   
y(:,2) = x(:,1);
y(:,3:6) = x(:,2:5);
" g0 I( g5 U/ _# s- j8 U3 Uf1 = y(:,2) - yexp(:,2);
f2 = y(:,3) - yexp(:,3);
f3 = y(:,4) - yexp(:,4);
( n0 |7 g4 Q2 a* [  Ff4 = y(:,5) - yexp(:,5);
9 r( Y9 s1 U: h- o- ff5 = y(:,6) - yexp(:,6);
f = [f1; f2; f3; f4; f5];


& g" ]9 f# C6 k# S* m' D* M
function f = ObjFunc7Fmincon(k,x0,yexp)
/ F* g% p  d) R- V% c$ \+ xtspan = [0 15 30 45 60 90 120 180 240 300 360];
+ P" f; P/ M0 j[t x] = ode45(@KineticEqs,tspan,x0,[],k);   
0 ^9 C! T+ X: ry(:,2) = x(:,1);
y(:,3:6) = x(:,2:5);
f =  sum((y(:,2)-yexp(:,2)).^2) + sum((y(:,3)-yexp(:,3)).^2)   ...
    + sum((y(:,4)-yexp(:,4)).^2) + sum((y(:,5)-yexp(:,5)).^2)   ...
" O9 l8 i: T/ m" {$ N+ j: i    + sum((y(:,6)-yexp(:,6)).^2) ;
3 e7 u0 i+ _# K" C+ t& c
1 V& k) ]6 a0 b* O/ g3 ?# u


2 L0 G5 B4 w: I4 bfunction dxdt = KineticEqs(t,x,k)
dGldt = k(1)*x(2)-(k(2)+k(3)+k(8))*x(1);
6 L. t4 D+ _, s. H; v+ @& t" vdFrdt = k(2)*x(1)-(k(1)+k(4)+k(5)+k(9))*x(2);
7 {4 \' t- G4 J* b5 ?dFadt = k(3)*x(1)+k(5)*x(2)+(k(6)+k(7))*x(5);
" v, }1 t/ ^, v5 S1 cdLadt = k(7)*x(5);
dHmdt = k(4)*x(2)-(k(6)+k(7)+k(10))*x(5);
dxdt = [dGldt; dFrdt; dFadt; dLadt; dHmdt];
/ `8 |" g, W& K* j' G
" k4 Q  k% J$ j# N1 r' u7 ?: ?2 t
# K! l9 Z9 D4 r3 R: V6 B* L