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

箫剑→残念

function parafit
%  k1->k-1，k2->k1，k3->k2,k4->k3,k5->k4
2 p& W7 f- t) @* Y% k6->k6 k7->k7
% dGlcdt = k-1*C(Fru)-(k1+k2)*C(Glc);
% dFrudt = k1*C(Glc)-(k-1+k3+k4)C(Fru);
" L% Z( U+ j' o8 t) c; ?% 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);
6 @9 K8 J' p) yclear all
clc
format 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
          45    0.15105    0.06975    0.02473    0    0.0033
" r6 C" t5 q" q          60    0.13763    0.07397    0.02615    0    0.00428
          90    0.08115    0.07877    0.07485    0    0.01405
$ X' }- ~% G# @$ Y# v! a          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
# ~* k' v1 A# X4 h( v% K# v          360    0.01855    0.04125    0.09363    0.0239    0.06495];
3 p) @, _9 _" s7 K& ak0 = [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];                  % 参数下限
ub = [1  1  1  1  1  1  1  1  1  1];    % 参数上限
: y% E- S9 D5 g, H% [! v" x) qx0 = [0.25  0  0  0  0];
yexp = Kinetics;                 % yexp: 实验数据[x1        x4        x5        x6]
4 }" |" B  D& t! D0 i" Q# I5 {8 i% 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))
- _- o2 e) g' F  z. ~4 S2 P8 rfprintf('\tk3 = %.11f\n',k(3))
5 p* b. g: l* \% v) gfprintf('\tk4 = %.11f\n',k(4))
: p) ~" T) d; H* L0 Pfprintf('\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))
fprintf('  The sum of the squares is: %.1e\n\n',fval)
8 B9 d* e  c9 M* a7 G+ Ek_fm= k;
% warning off
% 使用函数lsqnonlin()进行参数估计
; y8 `: \; a, Q+ K7 x[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc7LNL,k0,lb,ub,[],x0,yexp);      
- Y& }2 H( d: P6 p; [" \ci = nlparci(k,residual,jacobian);
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
7 m2 X! `, r+ [8 |  ]fprintf('\tk1 = %.11f\n',k(1))
! _! M5 q" b6 Rfprintf('\tk2 = %.11f\n',k(2))
2 ?) `' g) ~  _- m. `fprintf('\tk3 = %.11f\n',k(3))
0 x6 h% V) U5 E4 M0 Z( [1 q$ P- ]% S' Lfprintf('\tk4 = %.11f\n',k(4))
8 M" f/ P: l( Vfprintf('\tk5 = %.11f\n',k(5))
fprintf('\tk6 = %.11f\n',k(6))
fprintf('\tk7 = %.11f\n',k(7))
, I' ]. s+ J/ K6 G) @$ w0 Tfprintf('\tk8 = %.11f\n',k(8))
( |( ?7 a3 f. v4 g; e0 g* I3 X4 C; sfprintf('\tk9 = %.11f\n',k(9))
fprintf('\tk10 = %.11f\n',k(10))
$ n4 M6 y! S& s' jfprintf('  The sum of the squares is: %.1e\n\n',resnorm)
k_ls = k;
output
- f, A- A$ b( l; }& ^3 _warning off
/ }9 X; M, `* v0 u) U% 以函数fmincon()估计得到的结果为初值，使用函数lsqnonlin()进行参数估计
k0 = k_fm;
5 m+ ^+ {0 Z; Q$ I1 M; g0 F[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc7LNL,k0,lb,ub,[],x0,yexp);      
ci = nlparci(k,residual,jacobian);
) x) E* n* I0 d' vfprintf('\n\n以fmincon()的结果为初值，使用函数lsqnonlin()估计得到的参数值为:\n')
$ p4 e$ x' {9 w  A2 Ffprintf('\tk1 = %.11f\n',k(1))
fprintf('\tk2 = %.11f\n',k(2))
3 d+ u. Y+ d# n5 {$ ]7 A( bfprintf('\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))
fprintf('\tk8 = %.11f\n',k(8))
0 l# }! O7 n- F$ Z" v" lfprintf('\tk9 = %.11f\n',k(9))
fprintf('\tk10 = %.11f\n',k(10))
. T$ n9 H: V7 ~fprintf('  The sum of the squares is: %.1e\n\n',resnorm)
' Q7 v. B. H, t  Q( y& [k_fmls = k;
( \  A0 U8 Y* H/ ?, p! l% _output
tspan = [0 15 30 45 60 90 120 180 240 300 360];
% v; Y" c* S/ S6 r9 I7 K/ J[t x] = ode45(@KineticEqs,tspan,x0,[],k_fmls);
: G, p! z! I6 t/ d6 Xfigure;
plot(t,x(:,1),t,yexp(:,2),'*');legend('Glc-pr','Glc-real')
8 \* ^6 ^. [5 O* qfigure;plot(t,x(:,2:5));
p=x(:,1:5)
hold on
7 H! K" a5 a9 O3 m8 _plot(t,yexp(:,3:6),'o');legend('Fru-pr','Fa-pr','La-pr','HMF-pr','Fru-real','Fa-real','La-real','HMF-real')

* Q3 E7 A" [  E" k

function f = ObjFunc7LNL(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);
0 _1 I( `. u# R5 A7 {' o$ ^y(:,3:6) = x(:,2:5);
f1 = y(:,2) - yexp(:,2);
f2 = y(:,3) - yexp(:,3);
% S- I" \9 V4 v, ^f3 = y(:,4) - yexp(:,4);
f4 = y(:,5) - yexp(:,5);
f5 = y(:,6) - yexp(:,6);
f = [f1; f2; f3; f4; f5];

6 Z2 s3 f0 u- X0 C' j. @. P3 O
9 v' G! \) C/ y- G: o" w5 {
function f = ObjFunc7Fmincon(k,x0,yexp)
: X- o/ `& n  E7 @tspan = [0 15 30 45 60 90 120 180 240 300 360];
[t x] = ode45(@KineticEqs,tspan,x0,[],k);   
1 r$ T( L& D' {y(:,2) = x(:,1);
6 w' a: O; S' ty(:,3:6) = x(:,2:5);
0 v0 A6 c1 g. U% G3 {' B& @f =  sum((y(:,2)-yexp(:,2)).^2) + sum((y(:,3)-yexp(:,3)).^2)   ...
( c/ q2 T$ N/ B1 p3 h+ u$ a- G6 G    + sum((y(:,4)-yexp(:,4)).^2) + sum((y(:,5)-yexp(:,5)).^2)   ...
4 Z& x( m9 L1 K7 C7 k    + sum((y(:,6)-yexp(:,6)).^2) ;



$ D5 H- R8 z& Z4 S, K% J5 G
! \$ G! t! I4 S. Cfunction dxdt = KineticEqs(t,x,k)
! A# K9 l6 s) K1 t# f+ p. DdGldt = k(1)*x(2)-(k(2)+k(3)+k(8))*x(1);
- y7 v6 {" {1 u* O# I; I1 v; e& CdFrdt = 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);
* s+ h# E7 D$ b  L4 q% ^dLadt = k(7)*x(5);
dHmdt = k(4)*x(2)-(k(6)+k(7)+k(10))*x(5);
dxdt = [dGldt; dFrdt; dFadt; dLadt; dHmdt];

董事长之友

蒙的一比 大哥
" d+ |* X# d" _9 k( v* t1 ]8 L