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

箫剑→残念

function parafit
$ O7 o( L7 f$ I, @& {%  k1->k-1，k2->k1，k3->k2,k4->k3,k5->k4
% k6->k6 k7->k7
- g. w. C- A0 \; S% j) v% Z, C% 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);
% dLadt = k(7)*C(Hmf);
%dHmfdt = k(3)*C(Fru)-（k6+k7)*C(Hmf);
clear all
/ T% f- ?9 Q- c- B  Xclc
format long
5 {1 b  S; c2 B' M* d" M# N%        t/min   Glc    Fru        Fa   La   HMF/ mol/L
$ M9 N5 I0 f$ ~! ~' T  Kinetics=[0    0.25    0           0    0       0
; P7 H! z7 m' B: L, ?/ |" o. x          15    0.2319    0.01257    0.0048    0    2.50E-04
4 Q/ D0 g" w- O5 w) g- y3 B          30    0.19345    0.027    0.00868    0    7.00E-04
! D* L& Z  E* z9 c% i          45    0.15105    0.06975    0.02473    0    0.0033
( p2 |3 [! i. j          60    0.13763    0.07397    0.02615    0    0.00428
          90    0.08115    0.07877    0.07485    0    0.01405
2 n4 [7 E8 A0 e* u6 P7 {" V' ^7 E9 Q          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
$ x; f  w9 J6 [% ]9 q& f& a8 U          360    0.01855    0.04125    0.09363    0.0239    0.06495];
" d# H" [# L5 R- e) P( I& kk0 = [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];                  % 参数下限
" y2 y& g4 W- k. f( w; V# u5 `ub = [1  1  1  1  1  1  1  1  1  1];    % 参数上限
+ i- E% o$ f  l' _: Q  j% {! z- ~x0 = [0.25  0  0  0  0];
yexp = Kinetics;                 % yexp: 实验数据[x1        x4        x5        x6]
, _, i9 {" n* v/ a. ^% warning off
! ~! ^' U8 I$ E/ p, q8 k% 使用函数 ()进行参数估计
# S1 V& i7 ]( W! S0 R[k,fval,flag] = fmincon(@ObjFunc7Fmincon,k0,[],[],[],[],lb,ub,[],[],x0,yexp);
# @9 Y! A! G8 z# N+ z1 |2 xfprintf('\n使用函数fmincon()估计得到的参数值为:\n')
' P- o! f8 w  C+ v' `, sfprintf('\tk1 = %.11f\n',k(1))
( U8 E) \7 Z( `: s0 tfprintf('\tk2 = %.11f\n',k(2))
fprintf('\tk3 = %.11f\n',k(3))
fprintf('\tk4 = %.11f\n',k(4))
fprintf('\tk5 = %.11f\n',k(5))
fprintf('\tk6 = %.11f\n',k(6))
fprintf('\tk7 = %.11f\n',k(7))
# L  R2 x4 H0 efprintf('\tk8 = %.11f\n',k(8))
* d. z8 v2 d. s2 |fprintf('\tk9 = %.11f\n',k(9))
fprintf('\tk10 = %.11f\n',k(10))
# n6 @" ^: V4 h1 R& F1 `fprintf('  The sum of the squares is: %.1e\n\n',fval)
k_fm= k;
% warning off
8 p/ O& N, N$ {' s  S; P" G8 |% 使用函数lsqnonlin()进行参数估计
6 }6 B* M4 r4 S6 k[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc7LNL,k0,lb,ub,[],x0,yexp);      
ci = nlparci(k,residual,jacobian);
( N6 G1 W1 \5 D7 g5 Afprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\tk1 = %.11f\n',k(1))
3 r, o  ~* u" Z  x7 Y( Lfprintf('\tk2 = %.11f\n',k(2))
fprintf('\tk3 = %.11f\n',k(3))
* x0 x3 g2 k5 gfprintf('\tk4 = %.11f\n',k(4))
fprintf('\tk5 = %.11f\n',k(5))
$ j# W2 A5 N# [fprintf('\tk6 = %.11f\n',k(6))
fprintf('\tk7 = %.11f\n',k(7))
fprintf('\tk8 = %.11f\n',k(8))
# x( d$ Z$ Z4 q: E/ k( z. Bfprintf('\tk9 = %.11f\n',k(9))
# W# [1 a+ d5 }' Hfprintf('\tk10 = %.11f\n',k(10))
: V  b- w7 V, F& @" ufprintf('  The sum of the squares is: %.1e\n\n',resnorm)
9 ?2 k: x# O4 S; t& P* lk_ls = k;
: F% [1 [' W/ Y- q) Z8 goutput
warning off
3 H& d: Y2 y* `6 i% 以函数fmincon()估计得到的结果为初值，使用函数lsqnonlin()进行参数估计
k0 = k_fm;
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
) O6 q' H9 F: n5 L    lsqnonlin(@ObjFunc7LNL,k0,lb,ub,[],x0,yexp);      
ci = nlparci(k,residual,jacobian);
& u# K" P  \& v7 Jfprintf('\n\n以fmincon()的结果为初值，使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\tk1 = %.11f\n',k(1))
; W# X% r) c( l- o1 P# ufprintf('\tk2 = %.11f\n',k(2))
2 H2 e( x3 }4 W8 y) ufprintf('\tk3 = %.11f\n',k(3))
fprintf('\tk4 = %.11f\n',k(4))
fprintf('\tk5 = %.11f\n',k(5))
fprintf('\tk6 = %.11f\n',k(6))
( |  Q1 y# H* H& `! A& @fprintf('\tk7 = %.11f\n',k(7))
fprintf('\tk8 = %.11f\n',k(8))
fprintf('\tk9 = %.11f\n',k(9))
+ [$ \1 ^9 Y$ k0 Q0 Y6 o: yfprintf('\tk10 = %.11f\n',k(10))
) x+ Q6 |6 e) R- S8 Vfprintf('  The sum of the squares is: %.1e\n\n',resnorm)
k_fmls = k;
# e6 [3 x* Q0 l6 e8 Zoutput
1 G. W2 J# Z0 g& ?tspan = [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));
. }& ^, J/ y2 I5 C% k# dp=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')


# T1 q' ?2 A% t7 |" ?2 e$ x' m
1 F' v2 {$ A; O4 O' B7 b  z# yfunction f = ObjFunc7LNL(k,x0,yexp)
" v/ z5 E* x( E" P7 ?0 Rtspan = [0 15 30 45 60 90 120 180 240 300 360];
[t, x] = ode45(@KineticEqs,tspan,x0,[],k);   
y(:,2) = x(:,1);
# ]3 C% |$ K6 D+ _/ Uy(:,3:6) = x(:,2:5);
/ t2 b2 r8 h2 T) i2 Kf1 = y(:,2) - yexp(:,2);
f2 = y(:,3) - yexp(:,3);
: K7 p3 N) z6 {8 y& P  \( bf3 = y(:,4) - yexp(:,4);
8 F% H! v$ {! S( r! v6 W& c0 b; mf4 = y(:,5) - yexp(:,5);
f5 = y(:,6) - yexp(:,6);
: @) a1 ~7 g' L- {f = [f1; f2; f3; f4; f5];



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);   
y(:,2) = x(:,1);
y(:,3:6) = x(:,2:5);
7 k5 b9 Z2 E0 F9 O" Cf =  sum((y(:,2)-yexp(:,2)).^2) + sum((y(:,3)-yexp(:,3)).^2)   ...
    + sum((y(:,4)-yexp(:,4)).^2) + sum((y(:,5)-yexp(:,5)).^2)   ...
* T& W0 H: j  m1 h5 C* `    + sum((y(:,6)-yexp(:,6)).^2) ;

& U7 l. s7 `6 F) ]# F6 T) k
- S; V1 o% b) v! v/ q8 W

  W( F4 y9 ^4 ufunction dxdt = KineticEqs(t,x,k)
dGldt = k(1)*x(2)-(k(2)+k(3)+k(8))*x(1);
- M9 `: @$ [# B- o1 k2 ?9 P! B3 F! TdFrdt = 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);
" q# f+ {7 Z. C2 P% ?dLadt = k(7)*x(5);
: r8 I% O6 B/ _) j9 }& Y4 F+ UdHmdt = k(4)*x(2)-(k(6)+k(7)+k(10))*x(5);
dxdt = [dGldt; dFrdt; dFadt; dLadt; dHmdt];
7 p, b" w9 u5 ^% u
+ c" g+ L: H# J$ {1 T! Z