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

箫剑→残念

function parafit
) O; v7 t2 i- \%  k1->k-1，k2->k1，k3->k2,k4->k3,k5->k4
% k6->k6 k7->k7
% dGlcdt = k-1*C(Fru)-(k1+k2)*C(Glc);
) q4 c7 J) c+ K) n, [" V2 c  f% dFrudt = k1*C(Glc)-(k-1+k3+k4)C(Fru);
6 V' l$ n1 a1 x% dFadt = k(2)*C(Glc)+k4*C(Fru)+（k6+k7)*C(Hmf);
$ x3 |& u$ X# R& {% dLadt = k(7)*C(Hmf);
%dHmfdt = k(3)*C(Fru)-（k6+k7)*C(Hmf);
clear all
clc
format long
%        t/min   Glc    Fru        Fa   La   HMF/ mol/L
* J6 g4 O: Y& I- i5 u3 X/ Q# j  Kinetics=[0    0.25    0           0    0       0
          15    0.2319    0.01257    0.0048    0    2.50E-04
. b1 v, F0 t0 M/ w          30    0.19345    0.027    0.00868    0    7.00E-04
' J' F/ J6 @$ ?; _  W6 C" c7 ?- h          45    0.15105    0.06975    0.02473    0    0.0033
( Q( L. |) n& b- X+ q          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
          180    0.04488    0.0682    0.07135    0.0091    0.03623
- H( O- w4 C9 l$ M          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];
, f' q6 k7 p$ k$ p5 b: qk0 = [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];                  % 参数下限
- m1 O- f6 w1 K6 h. Y) }' Y* tub = [1  1  1  1  1  1  1  1  1  1];    % 参数上限
x0 = [0.25  0  0  0  0];
yexp = Kinetics;                 % yexp: 实验数据[x1        x4        x5        x6]
- u+ E% `/ V  L+ Y% warning off
* ~0 s  Y* T$ i$ u( P% 使用函数 ()进行参数估计
[k,fval,flag] = fmincon(@ObjFunc7Fmincon,k0,[],[],[],[],lb,ub,[],[],x0,yexp);
fprintf('\n使用函数fmincon()估计得到的参数值为:\n')
# t' Z6 y, k3 v& x( m' S3 f$ i& c# zfprintf('\tk1 = %.11f\n',k(1))
, @5 E7 L# T% d6 R! Ffprintf('\tk2 = %.11f\n',k(2))
: E* V3 ~9 W  H# L, sfprintf('\tk3 = %.11f\n',k(3))
, l2 }: E. \. |2 b% A' V6 Ifprintf('\tk4 = %.11f\n',k(4))
fprintf('\tk5 = %.11f\n',k(5))
fprintf('\tk6 = %.11f\n',k(6))
1 n& n! c- x% r" h. L- Lfprintf('\tk7 = %.11f\n',k(7))
7 m5 L1 ~7 c' t1 D0 J& Q7 @. qfprintf('\tk8 = %.11f\n',k(8))
& i5 M/ H+ O7 ^1 ifprintf('\tk9 = %.11f\n',k(9))
+ J! f- v! ~0 j3 x8 Gfprintf('\tk10 = %.11f\n',k(10))
) }" M8 P' W5 z$ L7 x/ e. ]fprintf('  The sum of the squares is: %.1e\n\n',fval)
k_fm= k;
) U5 I6 |  R. J0 \3 X% warning off
7 ]0 s- O5 T( k% 使用函数lsqnonlin()进行参数估计
$ Q+ v  D" a3 X7 u, b4 N[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc7LNL,k0,lb,ub,[],x0,yexp);      
4 q, N$ ]' Z% z! b5 Z9 P) Z5 U. Oci = nlparci(k,residual,jacobian);
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
! S4 U: i6 \$ U% Bfprintf('\tk1 = %.11f\n',k(1))
fprintf('\tk2 = %.11f\n',k(2))
- c* E% G& Y2 D, ^- Kfprintf('\tk3 = %.11f\n',k(3))
fprintf('\tk4 = %.11f\n',k(4))
" _, G7 [& E! q/ U  u1 w9 Mfprintf('\tk5 = %.11f\n',k(5))
fprintf('\tk6 = %.11f\n',k(6))
fprintf('\tk7 = %.11f\n',k(7))
* c! }5 H- U- @4 \" \7 zfprintf('\tk8 = %.11f\n',k(8))
fprintf('\tk9 = %.11f\n',k(9))
, ^% Y) y& [3 b2 w. X3 P% [' Gfprintf('\tk10 = %.11f\n',k(10))
fprintf('  The sum of the squares is: %.1e\n\n',resnorm)
4 n. r' f0 u# l" c4 uk_ls = k;
0 j, L& c  e0 A+ `- Goutput
2 ?- t. O, R4 d" wwarning off
7 [) h. H6 f0 W, J; ^( ]% 以函数fmincon()估计得到的结果为初值，使用函数lsqnonlin()进行参数估计
; Y2 Z( c% C! Y' I# D9 kk0 = k_fm;
7 A1 A1 O) t4 e2 I' s8 Y/ v; {[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
9 W! U/ o( Z1 U( ?2 y    lsqnonlin(@ObjFunc7LNL,k0,lb,ub,[],x0,yexp);      
& `3 E3 w* m6 s& @$ Kci = nlparci(k,residual,jacobian);
fprintf('\n\n以fmincon()的结果为初值，使用函数lsqnonlin()估计得到的参数值为:\n')
. l6 @2 N& u& c8 U8 zfprintf('\tk1 = %.11f\n',k(1))
' B/ k2 k" {6 z3 W+ P1 Y! Lfprintf('\tk2 = %.11f\n',k(2))
/ ?" O) ^4 k- ~8 [) cfprintf('\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))
6 x% L  r+ x' E: y$ J/ ufprintf('\tk8 = %.11f\n',k(8))
fprintf('\tk9 = %.11f\n',k(9))
3 [3 h8 Y7 X! ?% d# b. qfprintf('\tk10 = %.11f\n',k(10))
1 X% A4 C; u: P& jfprintf('  The sum of the squares is: %.1e\n\n',resnorm)
k_fmls = k;
output
tspan = [0 15 30 45 60 90 120 180 240 300 360];
' }5 Y' T6 U/ [3 H! X4 N3 L[t x] = ode45(@KineticEqs,tspan,x0,[],k_fmls);
8 Y- @- N+ I- e3 Dfigure;
plot(t,x(:,1),t,yexp(:,2),'*');legend('Glc-pr','Glc-real')
; u* l/ ?9 |4 a  P, o6 e  v4 ufigure;plot(t,x(:,2:5));
p=x(:,1:5)
+ s, W- _4 T2 xhold on
plot(t,yexp(:,3:6),'o');legend('Fru-pr','Fa-pr','La-pr','HMF-pr','Fru-real','Fa-real','La-real','HMF-real')
: i1 z. m1 c- V' J$ p+ X/ i8 f

( u0 L) V# S! A  g
function f = ObjFunc7LNL(k,x0,yexp)
2 Y0 Z3 z. K0 Xtspan = [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);
f1 = y(:,2) - yexp(:,2);
8 z) j6 C8 ~# _; `f2 = y(:,3) - yexp(:,3);
- i6 S' L, |; H* z: Mf3 = y(:,4) - yexp(:,4);
9 |4 l. I( W" U( t0 of4 = y(:,5) - yexp(:,5);
; F/ o/ `+ V* B* |% l) i- Tf5 = y(:,6) - yexp(:,6);
6 S7 u' m- s: {f = [f1; f2; f3; f4; f5];

$ O4 d/ j; _, I( z

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);   
8 p8 z/ E: Z; xy(:,2) = x(:,1);
y(:,3:6) = x(:,2:5);
f =  sum((y(:,2)-yexp(:,2)).^2) + sum((y(:,3)-yexp(:,3)).^2)   ...
1 e2 e% n. m+ Z+ b* ?- F3 ]    + sum((y(:,4)-yexp(:,4)).^2) + sum((y(:,5)-yexp(:,5)).^2)   ...
0 {5 o$ Z$ S, ], x! T+ v. g; L, Q    + sum((y(:,6)-yexp(:,6)).^2) ;
. S% R- R) u4 l
. }. E* L- {% F/ O


/ \, p" u1 w2 ~function dxdt = KineticEqs(t,x,k)
dGldt = k(1)*x(2)-(k(2)+k(3)+k(8))*x(1);
+ |/ A. w/ {" t4 d* T; BdFrdt = 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);
dLadt = k(7)*x(5);
7 X4 b. e( Q% Z! o5 w$ ^4 {dHmdt = k(4)*x(2)-(k(6)+k(7)+k(10))*x(5);
dxdt = [dGldt; dFrdt; dFadt; dLadt; dHmdt];
3 d+ c; T* i8 F- J
1 y# r8 a% K! e
4 F3 [; v" m* F