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

箫剑→残念

function parafit
9 P& b3 t4 @( m- ^9 v( \% Z%  k1->k-1，k2->k1，k3->k2,k4->k3,k5->k4
% k6->k6 k7->k7
% dGlcdt = k-1*C(Fru)-(k1+k2)*C(Glc);
. a0 h- B: N* m( ~+ _  \7 {% dFrudt = k1*C(Glc)-(k-1+k3+k4)C(Fru);
+ J: _% C1 f1 [% dFadt = k(2)*C(Glc)+k4*C(Fru)+（k6+k7)*C(Hmf);
$ }! D% e/ ^8 e9 w4 F* ~4 n4 v% dLadt = k(7)*C(Hmf);
4 C8 @! E2 p9 K( ?: S' ^3 R: P%dHmfdt = k(3)*C(Fru)-（k6+k7)*C(Hmf);
clear all
8 y2 \0 A' A+ C5 L1 \clc
+ l# h5 }* F% L$ hformat long
%        t/min   Glc    Fru        Fa   La   HMF/ mol/L
; P! I6 [1 y, c  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
, |5 ?' _. ^) v9 Z          45    0.15105    0.06975    0.02473    0    0.0033
8 X+ R, _& O; o2 R" I  S          60    0.13763    0.07397    0.02615    0    0.00428
% C" ?- C# F! o5 e# G/ n          90    0.08115    0.07877    0.07485    0    0.01405
4 H8 s& b3 v* \& B# q0 p          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
          360    0.01855    0.04125    0.09363    0.0239    0.06495];
( R/ I6 m0 c( a, o# l# pk0 = [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];    % 参数上限
& {9 s* _4 [; d2 m7 px0 = [0.25  0  0  0  0];
: D1 J" N4 D7 i% d  Ayexp = Kinetics;                 % yexp: 实验数据[x1        x4        x5        x6]
/ s4 t  o4 d, Z5 l% u) [$ ]% warning off
" [, k+ K( m( T% ]& M6 W% 使用函数 ()进行参数估计
; r. Q2 V4 F* \[k,fval,flag] = fmincon(@ObjFunc7Fmincon,k0,[],[],[],[],lb,ub,[],[],x0,yexp);
. r' o1 W& Y+ \fprintf('\n使用函数fmincon()估计得到的参数值为:\n')
7 e% r0 o1 K9 Q3 e" H8 C) _7 n1 P% ]fprintf('\tk1 = %.11f\n',k(1))
. K' [2 G5 r; D: k% x2 m' C. Y, Efprintf('\tk2 = %.11f\n',k(2))
/ F/ o8 O9 |9 Xfprintf('\tk3 = %.11f\n',k(3))
+ x; y0 i, n# I9 g% v$ Q+ gfprintf('\tk4 = %.11f\n',k(4))
! T/ f7 {- a. e% A2 f# v1 \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))
  |0 u  z5 k* r- @fprintf('\tk10 = %.11f\n',k(10))
4 B' _/ p' e8 B9 M8 wfprintf('  The sum of the squares is: %.1e\n\n',fval)
* M' I$ f" @0 L- `  K* Ak_fm= k;
( X/ U4 F  O; k" [1 ~# i% warning off
" ^& [$ c' W/ y1 i& h5 x% 使用函数lsqnonlin()进行参数估计
) P* F8 O" M# z! {$ [9 L. d[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
: c- k5 P! T/ m& L: C# b/ W    lsqnonlin(@ObjFunc7LNL,k0,lb,ub,[],x0,yexp);      
) k- \; q+ b$ B4 V% q8 ~: Lci = nlparci(k,residual,jacobian);
$ o" R0 v# H9 N: n, D! Ufprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\tk1 = %.11f\n',k(1))
! ~$ o7 t+ Q: t% y7 Q- @) sfprintf('\tk2 = %.11f\n',k(2))
5 _' E8 w# m- f! V3 o9 t9 Sfprintf('\tk3 = %.11f\n',k(3))
" h  m2 s, d1 A& f& }4 Yfprintf('\tk4 = %.11f\n',k(4))
6 ]; ]' [( z7 qfprintf('\tk5 = %.11f\n',k(5))
# F8 {7 I' I6 a: hfprintf('\tk6 = %.11f\n',k(6))
$ v" }7 }/ k4 I/ _! q3 zfprintf('\tk7 = %.11f\n',k(7))
fprintf('\tk8 = %.11f\n',k(8))
( l: b4 P5 p7 S; d& Mfprintf('\tk9 = %.11f\n',k(9))
fprintf('\tk10 = %.11f\n',k(10))
fprintf('  The sum of the squares is: %.1e\n\n',resnorm)
k_ls = k;
9 P% S! b# c: Y- Poutput
warning off
% 以函数fmincon()估计得到的结果为初值，使用函数lsqnonlin()进行参数估计
4 t  C2 |2 X+ Zk0 = k_fm;
. J0 n- G9 |0 Y1 w[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc7LNL,k0,lb,ub,[],x0,yexp);      
, a) a- ~, P( x8 R* z% lci = nlparci(k,residual,jacobian);
* w6 i" E1 p9 R. X) ?fprintf('\n\n以fmincon()的结果为初值，使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\tk1 = %.11f\n',k(1))
fprintf('\tk2 = %.11f\n',k(2))
& y, ]# ]7 p9 N6 [7 a* n8 l, p) Xfprintf('\tk3 = %.11f\n',k(3))
) c8 g5 W: N  c, P2 Y; Ofprintf('\tk4 = %.11f\n',k(4))
fprintf('\tk5 = %.11f\n',k(5))
# k1 V0 E2 t# P( b/ Ufprintf('\tk6 = %.11f\n',k(6))
9 a7 M8 W0 v0 }6 k2 Mfprintf('\tk7 = %.11f\n',k(7))
fprintf('\tk8 = %.11f\n',k(8))
fprintf('\tk9 = %.11f\n',k(9))
fprintf('\tk10 = %.11f\n',k(10))
9 P, h. [8 ?4 Q8 ]/ d/ o' zfprintf('  The sum of the squares is: %.1e\n\n',resnorm)
k_fmls = k;
output
4 n) v5 o9 A) G) L' Z& H# W: ~tspan = [0 15 30 45 60 90 120 180 240 300 360];
[t x] = ode45(@KineticEqs,tspan,x0,[],k_fmls);
1 W2 w. M0 d2 Yfigure;
plot(t,x(:,1),t,yexp(:,2),'*');legend('Glc-pr','Glc-real')
figure;plot(t,x(:,2:5));
( s: G# [! C, M9 Xp=x(:,1:5)
& H$ j' ^2 P; \: Z8 ohold 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 g6 i4 G9 |2 F+ X6 p6 v  V

9 z2 H$ B- @, G7 j# {# ]function f = ObjFunc7LNL(k,x0,yexp)
% m0 }  _% _# Ntspan = [0 15 30 45 60 90 120 180 240 300 360];
[t, x] = ode45(@KineticEqs,tspan,x0,[],k);   
y(:,2) = x(:,1);
" X2 }5 u% j) M* f+ ry(:,3:6) = x(:,2:5);
  m$ `  A  b8 A7 yf1 = y(:,2) - yexp(:,2);
f2 = y(:,3) - yexp(:,3);
! d3 ~6 h$ B9 U0 o; j& |f3 = y(:,4) - yexp(:,4);
f4 = y(:,5) - yexp(:,5);
( E( t& b9 U7 ?7 ?4 f) i3 Gf5 = y(:,6) - yexp(:,6);
f = [f1; f2; f3; f4; f5];
0 s. w9 b9 }' X; L
$ H  f, n/ W0 I: I: T9 T1 x
$ y) E: h% K7 O& {- S
3 k2 F3 Y+ w( f9 N+ ]% @function f = ObjFunc7Fmincon(k,x0,yexp)
tspan = [0 15 30 45 60 90 120 180 240 300 360];
) }5 C9 f" L1 b9 z7 t[t x] = ode45(@KineticEqs,tspan,x0,[],k);   
y(:,2) = x(:,1);
y(:,3:6) = x(:,2:5);
' r7 v3 Y7 J4 j3 S. k* xf =  sum((y(:,2)-yexp(:,2)).^2) + sum((y(:,3)-yexp(:,3)).^2)   ...
+ Q- D1 ?$ p4 S( j, M    + sum((y(:,4)-yexp(:,4)).^2) + sum((y(:,5)-yexp(:,5)).^2)   ...
4 T4 W( Q* {3 K$ g3 X! r4 o( l    + sum((y(:,6)-yexp(:,6)).^2) ;
' @6 V" Q( B( w2 D$ s

5 c# g* H* h( x' v
! d9 p4 z6 U: `- Z, {: V5 ]) D
function dxdt = KineticEqs(t,x,k)
+ n4 }8 n7 k! mdGldt = k(1)*x(2)-(k(2)+k(3)+k(8))*x(1);
9 i. d/ a, k1 O, J( JdFrdt = k(2)*x(1)-(k(1)+k(4)+k(5)+k(9))*x(2);
: i3 X  b+ g- z5 ^! h+ P$ T, KdFadt = k(3)*x(1)+k(5)*x(2)+(k(6)+k(7))*x(5);
dLadt = k(7)*x(5);
dHmdt = k(4)*x(2)-(k(6)+k(7)+k(10))*x(5);
+ a5 s% t* \8 f6 d, b: y! A6 ^$ Wdxdt = [dGldt; dFrdt; dFadt; dLadt; dHmdt];
6 g* n5 W% z4 N6 u( L3 Y' D5 V/ V# }

! y0 h) P8 l: U4 A; _, n+ @& o

董事长之友

蒙的一比 大哥