多元模型回归与分析教案.ppt

多元模型回归与分析 Still waters run deep.流静水深流静水深,人静心深人静心深 Where there is life,there is hope。有生命必有希望。有生命必有希望一、实验数据分析一、实验数据分析q由实验数据回归模型，得到模型参数前，对数据自变量间的由实验数据回归模型，得到模型参数前，对数据自变量间的线性相关性进行检验，是发现回归模型应用的可靠性和准确线性相关性进行检验，是发现回归模型应用的可靠性和准确性受限制的有效方法。性受限制的有效方法。因自变量间的线性相关性，使得无法区分它们对因变量的作用因自变量间的线性相关性，使得无法区分它们对因变量的作用;回归模型参数时会遇到几乎是奇异的数据矩阵，这样的模型参数有很回归模型参数时会遇到几乎是奇异的数据矩阵，这样的模型参数有很大的不确定性大的不确定性(95的参数置信度范围宽的参数置信度范围宽)。q例：例：回归二氧化硫的催化氧化速率方程：回归二氧化硫的催化氧化速率方程：装有载铂氧化铝催化剂颗粒的微分固定床反应器中，测定二氧装有载铂氧化铝催化剂颗粒的微分固定床反应器中，测定二氧化硫的催化氧化速率。总压为化硫的催化氧化速率。总压为790 mmHg时，记录流体相的组时，记录流体相的组成分压，有下表所示的速率结果，通过这些数据求取二氧化成分压，有下表所示的速率结果，通过这些数据求取二氧化硫的催化氧化速率方程。硫的催化氧化速率方程。2二氧化硫的催化氧化速率二氧化硫的催化氧化速率 r分 压(atm)mol/g.hSO3SO2O20.020.04280.02550.1860.040.03310.03530.1900.060.02720.04090.1930.080.02360.04430.1950.100.02140.04640.1960.120.02010.04760.197表表8 82 2 二氧化硫的催化氧化速率二氧化硫的催化氧化速率 3两种模型的非线性回归两种模型的非线性回归q1、一般的指数速率方程形式、一般的指数速率方程形式(8.2.1)k=0.517113.3；a=-1.987.02；b=-0.2164.556；c=6.078124.7 拟合结果：拟合结果：参数的参数的9595置信度太宽，模型参数不可靠置信度太宽，模型参数不可靠。q2、根据原子氧的吸附机理，得到的速率方程式（、根据原子氧的吸附机理，得到的速率方程式（Smith，Chemical Engineering Kinetics,3rd Ed.,1981,McGraw-Hill,P.374）(8.2.2)K=73K=73，为反应平衡常数，为反应平衡常数A=0.1017A=0.10170.09580.0958；B=16.02B=16.024.334.33 拟合结果：拟合结果：与方程与方程(8.2.1)(8.2.1)相比，方程参数的置信度有了显著改善。相比，方程参数的置信度有了显著改善。4对速率方程的进一步分析对速率方程的进一步分析 如果把方程如果把方程(8.2.2)改写为：改写为：(8.2.3)将模型参数代入计算并以方程左边为横坐标、右边为纵坐标将模型参数代入计算并以方程左边为横坐标、右边为纵坐标作图。作图。结果并不是斜结果并不是斜率为率为-1-1的直线。的直线。说明表所给的说明表所给的速率数据没有速率数据没有足够的信息来足够的信息来表明速率方程表明速率方程中的逆反应贡中的逆反应贡献。献。如将如将SO3SO3分压分压对对O2O2分压作分压作图，这两分图，这两分压间有近似压间有近似线性关系。线性关系。所以方程所以方程（8.2.18.2.1）的）的置信区间范置信区间范围大。围大。5二、回归模型的选择（二、回归模型的选择（1）q例：例：水饱和蒸汽压的模型回归水饱和蒸汽压的模型回归 水的蒸汽压数据选用的温度范围为水的蒸汽压数据选用的温度范围为0120 q三参数的三参数的Antoine方程：方程：q四参数的四参数的Riedel回归方程：回归方程：q五参数回归方程（参考五参数回归方程（参考Thek-Stiel的蒸汽压预测方程提出）：的蒸汽压预测方程提出）：(8.2.4)(8.2.5)(8.2.6)6水饱和蒸汽压的模型回归结果水饱和蒸汽压的模型回归结果参数Antoine方程改进Thek-Stiel方程A18.5587.5132B-3973.2-10.449C-39.9832.8683D-.064796E-6.8475R20.99999981.0表表8-38-3 水饱和蒸汽压的方程拟合结果水饱和蒸汽压的方程拟合结果 拟合度十分接近拟合度十分接近1 1，表明拟合是成功的，但实际上用，表明拟合是成功的，但实际上用AntoineAntoine方程来拟合回归得到的结果不理想，说明仅从拟方程来拟合回归得到的结果不理想，说明仅从拟合度上来判断结果的好坏是不够的。为什么呢？合度上来判断结果的好坏是不够的。为什么呢？7因变量与残差关系图因变量与残差关系图 q残差定义：残差定义：(8.2.7)q考察模型参数估计方法的两个基本假设：考察模型参数估计方法的两个基本假设：参数估计的误差相互不相关联，是随机的。参数估计的误差相互不相关联，是随机的。估计误差符合正态分布。估计误差符合正态分布。检查模型适合体系数据程度的最有效方法之一检查模型适合体系数据程度的最有效方法之一是对因变量与残差作图，观察其分布情况。是对因变量与残差作图，观察其分布情况。8Antoine方程拟合的残差方程拟合的残差 残差虽然很小，但其分残差虽然很小，但其分布不是随机的。布不是随机的。残差的分布同正态分布相残差的分布同正态分布相比，有较大的差距。比，有较大的差距。两方面的结果充分说明了拟合回归的两方面的结果充分说明了拟合回归的AntoineAntoine方程还不能充方程还不能充分反映蒸汽压与温度间的关系，造成残差间存在关联。分反映蒸汽压与温度间的关系，造成残差间存在关联。采用采用RiedelRiedel方程拟合得到的也是类似的结果。方程拟合得到的也是类似的结果。9改进改进Thek-Stiel方程方程的拟合结果方程方程的拟合结果拟合误差比拟合误差比AntoineAntoine方程方程小了近一个数量级，而且小了近一个数量级，而且残差分布是随机分布的。残差分布是随机分布的。误差分布基本符合正态误差分布基本符合正态分布。分布。改进改进Thek-StielThek-Stiel方程方程描述水饱和蒸汽压的合适模型。方程方程描述水饱和蒸汽压的合适模型。10二、回归模型的选择（二、回归模型的选择（2）q前面说明了前面说明了模型参数较少模型参数较少时会出现拟合残差时会出现拟合残差的分布不是随机的，而是呈现某种分布，相的分布不是随机的，而是呈现某种分布，相互关联。互关联。q在模型回归拟合数据的过程中，如在模型回归拟合数据的过程中，如模型参数模型参数过多过多会出现什么情况？如何判断回归拟合模会出现什么情况？如何判断回归拟合模型中有过多的参数呢？型中有过多的参数呢？11丙烷在氢型丝光沸石上的吸附平衡丙烷在氢型丝光沸石上的吸附平衡q例：例：选用不同吸附方程拟合丙烷在氢型丝光沸石体系选用不同吸附方程拟合丙烷在氢型丝光沸石体系303K的吸附平衡数据。的吸附平衡数据。q目标：目标：说明如何对模型拟合结果进行统计分析，确定模型拟合的好坏、模说明如何对模型拟合结果进行统计分析，确定模型拟合的好坏、模型参数的可靠性和准确性，从而进行拟合模型的选择。型参数的可靠性和准确性，从而进行拟合模型的选择。P,kPaq,mmol/gP,kPaq,mmol/gP,kPaq,mmol/gP,kPaq,mmol/g0.100.09 1.080.4812.670.81115.891.140.140.12 1.470.5116.700.85140.071.170.220.18 1.510.5324.810.90158.901.190.330.24 2.270.5934.280.95176.761.200.410.30 3.220.6443.850.98193.371.220.490.31 4.720.6954.621.02206.811.240.570.36 5.060.7065.791.040.770.41 7.390.7573.191.060.990.4410.260.7994.661.09表表8-4 303K8-4 303K时丙烷在氢型丝光沸石上的吸附平衡数据时丙烷在氢型丝光沸石上的吸附平衡数据 12具有代表性的、也是适用性较广的模型具有代表性的、也是适用性较广的模型 q1、Lanmuir(L)双参数方程：双参数方程：q2、Freundlich(F)双参数方程：双参数方程：(8.2.8)q3、BET双参数方程：双参数方程：q4、Langmuir-Freundlich(LF)三参数方程：三参数方程：q5、三参数方程：、三参数方程：q6、Toth三参数方程：三参数方程：q7、扩展的、扩展的LF方程（五参数）：方程（五参数）：q8、(14)式的特殊形式（四参数）式的特殊形式（四参数）：(8.2.9)(8.2.10)(8.2.11)(8.2.12)(8.2.13)(8.2.14)(8.2.15)13各模型的计算结果各模型的计算结果 Eq.(8)Eq.(9)Eq.(10)Eq.(11)Eq.(12)Eq.(13)Eq.(14)Eq.(15)nm1.0840.051/0.9760.0250.4380.0534.62317.221.5350.2570.7580.0880.7690.068a0.5530.1310.4460.034/1.3820.107/1.3290.4961.4270.144b/0.200.018/0.4940.0620.6781.6580.5490.0760.9420.1070.9750.058c/812.1116.3/20.9652.650.3410.0591.9071.5930.0170.003d/0.0240.0480.9120.046e/1.6380.586/s29.09010-27.72110-25.13410-23.79910-22.9933.15410-21.04810-28.68610-3R20.987900.991270.996140.997950.998580.998590.999860.99990表表8 85 5 吸附等温线关联的参数值、方差和回归系数吸附等温线关联的参数值、方差和回归系数 从表中可看出，方程从表中可看出，方程(8(814)14)拟合方差逐渐减少，回归系数更拟合方差逐渐减少，回归系数更接近接近1 1（方程（方程(12)(12)是通过压力数据来拟合的，故拟合方差和其是通过压力数据来拟合的，故拟合方差和其它方程的结果不是在同一数量级上它方程的结果不是在同一数量级上)。由方程。由方程(14)(14)的五参数形的五参数形式改进的方程式改进的方程(15)(15)式获得的结果最好，实验数据点几乎完全落式获得的结果最好，实验数据点几乎完全落在方程（在方程（1515）式的曲线上）式的曲线上(见下图见下图)。14方程方程(13)和方程和方程(15)的拟合结果的拟合结果 方程方程(15)(15)式获得的结果最好，实验数据点几乎完式获得的结果最好，实验数据点几乎完全落在方程（全落在方程（1515）式的曲线上。）式的曲线上。15判断模型参数是否过少的依据判断模型参数是否过少的依据q通过对方程通过对方程(13)和五参数方程和五参数方程(15)的残差进行分析，的残差进行分析，方程方程(13)因参数过少，吸附量的计算误差与实验吸附量之间存在着某种分布。因参数过少，吸附量的计算误差与实验吸附量之间存在着某种分布。方程方程(15)计算误差在零的两边是随机分布的，看不出规律性。因此，拟合计算计算误差在零的两边是随机分布的，看不出规律性。因此，拟合计算误差有无规律性的分布是判断模型参数是否过少的依据。误差有无规律性的分布是判断模型参数是否过少的依据。因此，拟合计算误差有无规律性的分布是判断模型参数是因此，拟合计算误差有无规律性的分布是判断模型参数是否过少的依据。否过少的依据。方程方程(13)(13)的拟合误差的拟合误差 方程方程(15)(15)的拟合误差的拟合误差 16判断模型参数是否过多的依据判断模型参数是否过多的依据Eq.(14)Eq.(15)nm0.7580.0880.7690.068a1.3290.4961.4270.144b0.9420.1070.9750.058c1.9071.5930.0170.003d0.0240.0480.9120.046e1.6380.586/s21.04810-28.68610-3R20.999860.99990q在方程在方程(14)的计算结果中，有的计算结果中，有些参数些参数95%的置信度较大，说的置信度较大，说明这些参数之间有联系，不是明这些参数之间有联系，不是独立的。独立的。q而对于方程而对于方程(14)的五参数形式，的五参数形式，即方程即方程(15)，其所有参数的，其所有参数的95%置信度都较小。事实上，置信度都较小。事实上，方程方程(15)就是据此分析对吸附就是据此分析对吸附平衡理论作进一步研究而获得平衡理论作进一步研究而获得的。的。因此，拟合参数因此，拟合参数95%95%的置信度是否较大是判断模型参数是否的置信度是否较大是判断模型参数是否过多的依据。过多的依据。17回归模型的选择总结回归模型的选择总结q模型参数较少模型参数较少时会出现拟合残差的分布不是时会出现拟合残差的分布不是随机的，而是呈现某种分布，相互关联。残随机的，而是呈现某种分布，相互关联。残差的分布偏离正态分布较远。差的分布偏离正态分布较远。q模型参数过多模型参数过多会出现某些参数会出现某些参数95%的置信度的置信度较大，说明这些参数之间有联系，不是独立较大，说明这些参数之间有联系，不是独立的。的。18习题习题q研究二硫化碳饱和蒸汽压的模型回归问题。（研究二硫化碳饱和蒸汽压的模型回归问题。（P266，Ex8.3）二硫化碳的基本性质：二硫化碳的基本性质：临界温度为临界温度为273.05273.05 临界压力为临界压力为72.87 atm72.87 atm。温度，蒸汽压，mmHg温度，蒸汽压，mmHg-701.610198.0-603.520297.5-507.130432.7-4014.040616.7-3026.250995.6-2046.5601170.4-1078.8701558.00127.319Statistica的非线性估计的非线性估计q非线性估计方法非线性估计方法User-Specified Regression,least square 可以计算95%置信区间20“least square”与与“Custom Loss”比较比较q“least square”计算结果（采用计算结果（采用Levenberg-Marquardt方法）方法）q“Custom Loss”计算结果（采用计算结果（采用Quasi-Newton法）法）Matrix ill conditioned;cannot compute standard errors.EConf21“Custom Loss”的方差分析与迭代步骤的方差分析与迭代步骤 “Custom Custom Loss Loss”无迭代无迭代历史纪录，协历史纪录，协方差分析结果方差分析结果已出现病态。已出现病态。Covariance matrix cannot be computed.22Statistica非线性估计的残差分析非线性估计的残差分析残差分布情况残差分布情况残差对预测值作图残差对预测值作图23对方程对方程(15)的残差分析的残差分析Histogram of Histogram of residualsresidualsResidual vs.Residual vs.PredictedPredicted24误差正态分布图误差正态分布图“least square”计算的误差正态分布图计算的误差正态分布图“Custom Loss”计算的误差正态分布图计算的误差正态分布图qAntoine 方程拟合结果方程拟合结果25非线性函数的管理非线性函数的管理26三、三、MATLAB的拟合函数的拟合函数q多项式拟合函数多项式拟合函数polyfitq非线性最小二乘法非线性最小二乘法lsqnonlin()非线性最小二乘（优化问题）非线性最小二乘（优化问题）lsqcurvefit()非线性最小二乘曲线拟合非线性最小二乘曲线拟合nlinfit()前两种的简化版本前两种的简化版本nlparci()计算参数的置信区间计算参数的置信区间nlpredci()计算预测值的置信区间计算预测值的置信区间nlintool()nlinfit()、nlparci()、nlpredci()的集成图形用的集成图形用户界面拟合函数户界面拟合函数q注意：不同的拟合函数命令，其优化目标函数定义以注意：不同的拟合函数命令，其优化目标函数定义以及调用形式不同，注意区分！及调用形式不同，注意区分！27（一）（一）Polyfit的使用的使用qp=polyfit(x0,y0,n)其中其中x0和和y0分别为观察节点和观察值向量；分别为观察节点和观察值向量；n表示插值多项式的次数；表示插值多项式的次数；输出值输出值p表示插值多项式的系数。表示插值多项式的系数。q例：某实验中测得一组数据，其值如下：例：某实验中测得一组数据，其值如下：x l 2 3 4 5 y 1.3 1.8 2.2 2.9 3.5 已知已知x和和y成线性关系，即成线性关系，即y=kx+b，求系数，求系数k和和b x=1 2 3 4 5;y=1.3 1.8 2.2 2.9 3.5;p=polyfit(x,y,1)y1=polyval(p,x);plot(x,y1)hold on;plot(x,y,b*);p=0.55 0.69也就是，也就是，k为为0.55，b为为0.6928（二）（二）lsqcurvefit的使用的使用q方程（目标函数）方程（目标函数）Find coefficients beta that best fit the equationq调用形式调用形式beta=lsqcurvefit(fun,beta0,xdata,ydata)beta=lsqcurvefit(fun,beta0,xdata,ydata,lb,ub)beta=lsqcurvefit(fun,beta0,xdata,ydata,lb,ub,options)xdata和和ydata分别为观察节点和观察值向量；分别为观察节点和观察值向量；fun自定义的非线性拟合模型；自定义的非线性拟合模型；beta0拟合参数的初始值；拟合参数的初始值；beta拟合模型中的参数；拟合模型中的参数；lb,ub拟合参数初值的边界值拟合参数初值的边界值,lb=beta 1%two output arguments J=.%Jacobian of the function evaluated at betaend options=optimset(Jacobian,on)qFor examplefunction F=myfun(beta,xdata)F=beta(1)*exp(beta(2)*xdata);32lsqcurvefit应用示例应用示例%Assume you determined xdata and ydata experimentallyxdata=0.9 1.5 13.8 19.8 24.1 28.2 35.2 60.3 74.6 81.3;ydata=455.2 428.6 124.1 67.3 43.2 28.1 13.1-0.4-1.3-1.5;beta0=100;-1%Starting guessbeta,resnorm=lsqcurvefit(myfun,beta0,xdata,ydata)function F=myfun(beta,xdata)F=beta(1)*exp(beta(2)*xdata);beta=498.8309-0.1013 resnorm=9.504933（三）（三）lsqnonlin的使用的使用q方程（目标函数）方程（目标函数）Find coefficients beta that best fit the equationq调用形式调用形式beta=lsqnonlin(fun,beta 0)beta=lsqnonlin(fun,beta 0,lb,ub)beta=lsqnonlin(fun,beta 0,lb,ub,options)fun自定义的优化函数；自定义的优化函数；beta 0优化参数的初始值；优化参数的初始值；beta拟合模型中的参数；拟合模型中的参数；lb,ub拟合参数初值的边界值拟合参数初值的边界值,lb=beta s=dsolve(Dy1=-k1*y1,Dy2=k1*y1-k2*y2,y1(0)=1,y2(0)=0)s=y1:1x1 sym y2:1x1 sym s.y1 s.y2 simplify(s.y2)ans=exp(-k1*t)ans=-(-k1+k2)/(k1-k2)2*k1*exp(-k2*t)-1/(k1-k2)*k1*exp(-k1*t)ans=k1*(exp(-k2*t)-exp(-k1*t)/(k1-k2)41对积分式进行参数估计的源程序对积分式进行参数估计的源程序 function seqcurvefit11clear all;load seqdata;beta0=0.005 0.001;lb=0 0;ub=inf inf;beta,resnorm,residual,exitflag,output,lambda,jacobian=.lsqnonlin(seqfun,beta0,lb,ub,t,c);ci=nlparci(beta,residual,jacobian);function y=seqfun(beta,t,c)ca=exp(-beta(1)*t);cb=beta(1)/(beta(2)-beta(1)*(exp(-beta(1)*t)-exp(-beta(2)*t);y=ca-c(:,1)cb-c(:,2);42输出计算结果的源程序输出计算结果的源程序%print resultfprintf(n Estimated Parameters by Lsqnonlin():n)fprintf(t k1=%.4f%.4fn,beta(1),ci(1,2)-beta(1)fprintf(t k2=%.4f%.4fn,beta(2),ci(2,2)-beta(2)fprintf(The sum of the squares is:%.1enn,sum(residual.2)%plot of fit resultstc=linspace(0,max(t),200);y_row,y_col=size(c);zeroc=zeros(200,y_col);yc=seqfun(beta,tc,zeroc);plot(t,c(:,1),ro,tc,yc(:,1),r-);hold onplot(t,c(:,2),b+,tc,yc(:,2),b-);xlabel(Time);ylabel(Concentration);hold off用函数求拟合值。用函数求拟合值。注意转置，与函数定义注意转置，与函数定义的矩阵维数一致！的矩阵维数一致！43参数拟合结果参数拟合结果 seqcurvefit11Optimization terminated:relative function value changing by less than OPTIONS.TolFun.Estimated Parameters by Lsqnonlin():k1=0.0412 0.0018 k2=0.0121 0.0008 The sum of the residual squares is:2.7e-00344方法方法2：对微分式进行参数估计：对微分式进行参数估计function seqcurvefit21clear all;load seqdata;y_row,y_col=size(c);beta0=0.005 0.001;c0=1 0;lb=0 0;ub=inf inf;beta,resnorm,residual,exitflag,output,lambda,jacobian=.lsqnonlin(seqfun,beta0,lb,ub,t,c,y_col,c0);ci=nlparci(beta,residual,jacobian);45拟合模型和目标函数的定义拟合模型和目标函数的定义function y=seqfun(beta,t,c,y_col,c0)%Objective functiontspan=0 max(t);tt yy=ode45(modeleqs,tspan,c0,beta);for col=1:y_col yc(:,col)=spline(tt,yy(:,col),t);endy=c(:,1)-yc(:,1);c(:,2)-yc(:,2);function dydt=modeleqs(t,y,beta)%Model equationdydt=-beta(1)*y(1);beta(1)*y(1)-beta(2)*y(2);beta,resnorm,residual,exitflag,output,lambda,jacobian=.lsqnonlin(seqfun,beta0,lb,ub,t,c,y_col,c0);46输出计算结果的源程序输出计算结果的源程序%print resultfprintf(n Estimated Parameters by Lsqnonlin():n)fprintf(t k1=%.4f%.4fn,beta(1),ci(1,2)-beta(1)fprintf(t k2=%.4f%.4fn,beta(2),ci(2,2)-beta(2)fprintf(The sum of the squares is:%.1enn,sum(residual.2)%plot of fit resultstspan=0 max(t);tt yc=ode45(modeleqs,tspan,c0,beta);tc=linspace(0,max(t),200);yca=spline(tt,yc(:,1),tc);plot(t,c(:,1),ro,tc,yca,r-);hold onycb=spline(tt,yc(:,2),tc);plot(t,c(:,2),b+,tc,ycb,b-);xlabel(Time);ylabel(Concentration);hold off 首先解微分方程求拟合首先解微分方程求拟合值。然后用样条插值逼值。然后用样条插值逼近。近。47参数拟合结果参数拟合结果 seqcurvefit21Optimization terminated:relative function value changing by less than OPTIONS.TolFun.Estimated Parameters by Lsqnonlin():k1=0.0412 0.0018 k2=0.0121 0.0008 The sum of the residual squares is:5.0e-003 48使用使用MATLAB进行参数估计示例进行参数估计示例q例例2：青霉素发酵过程动力学参数：青霉素发酵过程动力学参数估计：估计：在间歇发酵罐中研究青霉素发酵在间歇发酵罐中研究青霉素发酵过程动力学，微生物过程动力学，微生物Penicillium chrysogenum在一定的控制条件在一定的控制条件下生长，细胞生长速率可以用逻下生长，细胞生长速率可以用逻辑模型描述：辑模型描述：青霉素的生产速率模型为：青霉素的生产速率模型为：q初始条件为：初始条件为：t=0，y1=0.29，y2=0。TimehoursCellconcentrationdry weightPenicillinconcentrationunits/ml00.180100.120220.480.0089341.460.0642461.560.2266581.730.4373701.990.6943822.621.2459942.881.43151063.432.04021183.371.92781303.922.18481423.962.42041543.582.46151663.582.2831783.342.70781903.472.654249使用使用MATLAB的的dsolve求积分求积分 s=dsolve(Dy1=k1*y1*(1-y1/k1),Dy2=k3*y1-k4*y2,y1(0)=0.29,y2(0)=0)s=y1:1x1 sym y2:1x1 sym s.y1 s.y2 simplify(s.y2)ans=k1/(1+1/29*exp(-k1*t)*(100*k1-29)ans=Int(k3*k1/(1+1/29*exp(-k1*_z1)*(100*k1-29)*exp(k4*_z1),_z1=0.t)*exp(-k4*t)ans=29*k3*k1*Int(exp(k4*_z1)/(29+100*exp(-k1*_z1)*k1-29*exp(-k1*_z1),_z1=0.t)*exp(-k4*t)没有确定的积分表达式没有确定的积分表达式50参数估计的源程序参数估计的源程序function PenicilliumEstclear all;load PenicilliumData%Load experimental data y_row,y_col=size(y);%Nonlinear least square estimate using lsqnonlin()beta0=0.1 4.0 0.02 0.02;y0=0.29 0.0;lb=0 0.01 0 0;ub=inf inf inf inf;beta,resnorm,residual,exitflag,output,lambda,jacobian=.lsqnonlin(Func,beta0,lb,ub,x,y,y_col,y0);ci=nlparci(beta,residual,jacobian);51拟合模型和目标函数的定义拟合模型和目标函数的定义%=function f=Func(beta,x,y,y_col,y0)%Define objective functiontspan=0 max(x);tt yy=ode45(ModelEqs,tspan,y0,beta);for col=1:y_col yc(:,col)=spline(tt,yy(:,col),x);endf1=y(:,1)-yc(:,1);f2=y(:,2)-yc(:,2);f=f1;f2;%=function dydt=ModelEqs(t,y,beta)%Model equationsdydt=beta(1)*y(1)*(1-y(1)/beta(2);beta(3)*y(1)-beta(4)*y(2);52输出计算结果的源程序输出计算结果的源程序%resultfprintf(n Estimated Parameters by Lsqnonlin():n)fprintf(t k1=%.4f%.4fn,beta(1),ci(1,2)-beta(1)fprintf(t k2=%.4f%.4fn,beta(2),ci(2,2)-beta(2)fprintf(t k3=%.4f%.4fn,beta(3),ci(3,2)-beta(3)fprintf(t k4=%.4f%.4fn,beta(4),ci(4,2)-beta(4)fprintf(The sum of the residual squares is:%.1enn,sum(residual.2)53输出计算结果的源程序输出计算结果的源程序%plotnamex=Time,hours;namey(1,1:length(Cell concentration,%dry weight)=.Cell concentration,%dry weight;namey(2,1:length(Pencillin concentration,units/ml)=.Pencillin concentration,units/ml;tspan=0 max(x);tt yy=ode45(ModelEqs,tspan,y0,beta);xc=linspace(0,max(x),200);for col=1:y_col yc(:,col)=spline(tt,yy(:,col),xc);figure(col)plot(x,y(:,col),o,xc,yc(:,col)xlabel(namex)ylabel(namey(col,:)end 54参数拟合结果参数拟合结果 penicilliumestOptimization terminated:relative function value changing by less than OPTIONS.TolFun.Estimated Parameters by Lsqnonlin():k1=0.0423 0.0040 k2=3.6850 0.1949 k3=0.0186 0.0081 k4=0.0235 0.0142 The sum of the residual squares is:1.4e+000 55参数拟合结果参数拟合结果56