多元模型回归和分析

多元模型回归和分析2一、实验数据分析由实验数据回归模型，得到模型参数前，对数据自变量间的线性相关性进行检验，是发现回归模型应用的可靠性和准确性受限制的有效方法。因自变量间的线性相关性，使得无法区分它们对因变量的作用;回归模型参数时会遇到几乎是奇异的数据矩阵，这样的模型参数有很大的不确定性(95％的参数置信度范围宽)。例：回归二氧化硫的催化氧化速率方程：装有载铂氧化铝催化剂颗粒的微分固定床反应器中，测定二氧化硫的催化氧化速率。总压为790mmHg时，记录流体相的组成分压，有下表所示的速率结果，通过这些数据求取二氧化硫的催化氧化速率方程。3二氧化硫的催化氧化速率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－2二氧化硫的催化氧化速率4两种模型的非线性回归1、一般的指数速率方程形式(8.2.1)k=0.517±113.3；a=-1.98±7.02；b=-0.216±4.556；c=6.078±124.7

拟合结果：参数的95％置信度太宽，模型参数不可靠。

2、根据原子氧的吸附机理，得到的速率方程式（Smith，ChemicalEngineeringKinetics,3rdEd.,1981,McGraw-Hill,P.374）(8.2.2)K=73，为反应平衡常数A=0.1017±0.0958；B=16.02±4.33

拟合结果：与方程(8.2.1)相比，方程参数的置信度有了显著改善。

5对速率方程的进一步分析如果把方程(8.2.2)改写为：(8.2.3)将模型参数代入计算并以方程左边为横坐标、右边为纵坐标作图。结果并不是斜率为-1的直线。说明表所给的速率数据没有足够的信息来表明速率方程中的逆反应贡献。如将SO3分压对O2分压作图，这两分压间有近似线性关系。所以方程（）的置信区间范围大。6二、回归模型的选择（1）例：水饱和蒸汽压的模型回归水的蒸汽压数据选用的温度范围为0～120℃三参数的Antoine方程：四参数的Riedel回归方程：五参数回归方程（参考Thek-Stiel的蒸汽压预测方程提出）：(8.2.4)(8.2.5)(8.2.6)7水饱和蒸汽压的模型回归结果参数Antoine方程改进Thek-Stiel方程A18.5587.5132B-3973.2-10.449C-39.9832.8683D-.064796E-6.8475R20.99999981.0表8-3 水饱和蒸汽压的方程拟合结果

拟合度十分接近1，表明拟合是成功的，但实际上用Antoine方程来拟合回归得到的结果不理想，说明仅从拟合度上来判断结果的好坏是不够的。为什么呢？

8因变量与残差关系图残差定义：(8.2.7)考察模型参数估计方法的两个基本假设：参数估计的误差相互不相关联，是随机的。估计误差符合正态分布。检查模型适合体系数据程度的最有效方法之一是对因变量与残差作图，观察其分布情况。

9Antoine方程拟合的残差残差虽然很小，但其分布不是随机的。残差的分布同正态分布相比，有较大的差距。

两方面的结果充分说明了拟合回归的Antoine方程还不能充分反映蒸汽压与温度间的关系，造成残差间存在关联。采用Riedel方程拟合得到的也是类似的结果。

10改进Thek-Stiel方程方程的拟合结果拟合误差比Antoine方程小了近一个数量级，而且残差分布是随机分布的。误差分布基本符合正态分布。改进Thek-Stiel方程方程描述水饱和蒸汽压的合适模型。

11二、回归模型的选择（2）前面说明了模型参数较少时会出现拟合残差的分布不是随机的，而是呈现某种分布，相互关联。在模型回归拟合数据的过程中，如模型参数过多会出现什么情况？如何判断回归拟合模型中有过多的参数呢？12丙烷在氢型丝光沸石上的吸附平衡例：选用不同吸附方程拟合丙烷在氢型丝光沸石体系303K的吸附平衡数据。目标：说明如何对模型拟合结果进行统计分析，确定模型拟合的好坏、模型参数的可靠性和准确性，从而进行拟合模型的选择。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-4303K时丙烷在氢型丝光沸石上的吸附平衡数据

13具有代表性的、也是适用性较广的模型1、Lanmuir(L)双参数方程：2、Freundlich(F)双参数方程：(8.2.8)3、BET双参数方程：4、Langmuir-Freundlich(LF)三参数方程：5、三参数方程：6、Toth三参数方程：7、扩展的LF方程（五参数）：8、(14)式的特殊形式（四参数）：(8.2.9)(8.2.10)(8.2.11)(8.2.12)(8.2.13)(8.2.14)(8.2.15)14各模型的计算结果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－5吸附等温线关联的参数值、方差和回归系数

从表中可看出，方程(8→14)拟合方差逐渐减少，回归系数更接近1（方程(12)是通过压力数据来拟合的，故拟合方差和其它方程的结果不是在同一数量级上)。由方程(14)的五参数形式改进的方程(15)式获得的结果最好，实验数据点几乎完全落在方程（15）式的曲线上(见下图)。15方程(13)和方程(15)的拟合结果方程(15)式获得的结果最好，实验数据点几乎完全落在方程（15）式的曲线上。16判断模型参数是否过少的依据通过对方程(13)和五参数方程(15)的残差进行分析，方程(13)因参数过少，吸附量的计算误差与实验吸附量之间存在着某种分布。方程(15)计算误差在零的两边是随机分布的，看不出规律性。因此，拟合计算误差有无规律性的分布是判断模型参数是否过少的依据。因此，拟合计算误差有无规律性的分布是判断模型参数是否过少的依据。

方程(13)的拟合误差

方程(15)的拟合误差

17判断模型参数是否过多的依据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.99990在方程(14)的计算结果中，有些参数95%的置信度较大，说明这些参数之间有联系，不是独立的。而对于方程(14)的五参数形式，即方程(15)，其所有参数的95%置信度都较小。事实上，方程(15)就是据此分析对吸附平衡理论作进一步研究而获得的。因此，拟合参数95%的置信度是否较大是判断模型参数是否过多的依据。18回归模型的选择总结模型参数较少时会出现拟合残差的分布不是随机的，而是呈现某种分布，相互关联。残差的分布偏离正态分布较远。模型参数过多会出现某些参数95%的置信度较大，说明这些参数之间有联系，不是独立的。19习题研究二硫化碳饱和蒸汽压的模型回归问题。（P266，Ex8.3）二硫化碳的基本性质：临界温度为273.05℃临界压力为72.87atm。温度，℃蒸汽压，mmHg温度，℃蒸汽压，mmHg-701.610198.0-603.520297.5-507.130432.7-4014.040616.7-3026.250995.6-2046.5601170.4-1078.8701558.00127.320Statistica的非线性估计非线性估计方法User-SpecifiedRegression,leastsquare

可以计算95%置信区间21“leastsquare”与“CustomLoss

”比较“leastsquare”计算结果（采用Levenberg-Marquardt方法）“CustomLoss

”计算结果（采用Quasi-Newton法） Matrixillconditioned;cannotcomputestandarderrors.E±Conf22“CustomLoss

”的方差分析与迭代步骤

“CustomLoss”无迭代历史纪录，协方差分析结果已出现病态。Covariancematrixcannotbecomputed.23Statistica非线性估计的残差分析残差分布情况残差对预测值作图24对方程(15)的残差分析HistogramofresidualsResidualvs.Predicted25误差正态分布图“leastsquare”计算的误差正态分布图“CustomLoss

”计算的误差正态分布图Antoine方程拟合结果26非线性函数的管理27三、MATLAB的拟合函数多项式拟合函数polyfit非线性最小二乘法lsqnonlin()——非线性最小二乘（优化问题）lsqcurvefit()——非线性最小二乘曲线拟合nlinfit()——前两种的简化版本nlparci()——计算参数的置信区间nlpredci()——计算预测值的置信区间nlintool()——nlinfit()、nlparci()、nlpredci()的集成图形用户界面拟合函数注意：不同的拟合函数命令，其优化目标函数定义以及调用形式不同，注意区分！！！28（一）Polyfit的使用p=polyfit(x0,y0,n)其中x0和y0分别为观察节点和观察值向量；n表示插值多项式的次数；输出值p表示插值多项式的系数。例：某实验中测得一组数据，其值如下：

xl2345y1.31.82.22.93.5已知x和y成线性关系，即y=kx+b，求系数k和b

x=[12345];y=[1.31.82.22.93.5];p=polyfit(x,y,1)y1=polyval(p,x);plot(x,y1)holdon;plot(x,y,'b*');p=0.550.69也就是，k为0.55，b为0.6929（二）lsqcurvefit的使用方程（目标函数）Findcoefficientsbetathatbestfittheequation调用形式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<=ub。30lsqcurvefit的使用（续）调用形式（续）[beta,resnorm]=lsqcurvefit(...)[beta,resnorm,residual]=lsqcurvefit(...)[beta,resnorm,residual,exitflag]=lsqcurvefit(...)[beta,resnorm,residual,exitflag,output]=lsqcurvefit(...)[beta,resnorm,residual,exitflag,output,lambda]=lsqcurvefit(...)[beta,resnorm,residual,exitflag,output,lambda,jacobian]=lsqcurvefit(...)resnorm——返回beta处的残差平方和，sum((fun(x,xdata)-ydata).^2)residual——返回解beta处的残差，

fun(x,xdata)-ydata

exitflag——退出方式output——返回优化信息的输出结果，iterations、funcCount、algorithm、stepsize等lambda——解beta处的Lagrange乘子jacobian——返回函数在解beta处的Jacobian矩阵31lsqcurvefit的拟合函数的定义拟合函数的定义：functionF=myfun(beta,xdata) F=...%Computefunctionvaluesatxbeta=lsqcurvefit(@myfun,beta0,xdata,ydata)Note:(1)funshouldreturnfun(x,xdata),andnotthesum-of-squaressum((fun(x,xdata)-ydata).^2).(2)Thealgorithmimplicitlysquaresandsumsfun(x,xdata)-ydata.32lsqcurvefit的拟合函数的定义（续）IftheJacobiancanalsobecomputedbyuser-definedfunction[F,J]=myfun(beta,xdata)F=...%objectivefunctionvaluesatbetaifnargout>1%twooutputargumentsJ=...%Jacobianofthefunctionevaluatedatbetaendoptions=optimset('Jacobian','on')ForexamplefunctionF=myfun(beta,xdata) F=beta(1)*exp(beta(2)*xdata);33lsqcurvefit应用示例%Assumeyoudeterminedxdataandydataexperimentallyxdata=[0.91.513.819.824.128.235.260.374.681.3];ydata=[455.2428.6124.167.343.228.113.1-0.4-1.3-1.5];beta0=[100;-1]%Startingguess[beta,resnorm]=lsqcurvefit(@myfun,beta0,xdata,ydata)functionF=myfun(beta,xdata)F=beta

(1)*exp(beta(2)*xdata);beta=498.8309-0.1013resnorm=9.504934（三）lsqnonlin的使用方程（目标函数）Findcoefficientsbetathatbestfittheequation调用形式beta=lsqnonlin(fun,beta0)beta=lsqnonlin(fun,beta0,lb,ub)beta=lsqnonlin(fun,beta0,lb,ub,options)fun——自定义的优化函数；beta

0——优化参数的初始值；beta——拟合模型中的参数；lb,ub——拟合参数初值的边界值,lb<=beta<=ub。35lsqnonlin的使用（续）调用形式（续）[x,resnorm]=lsqnonlin(...)[x,resnorm,residual]=lsqnonlin(...)[x,resnorm,residual,exitflag]=lsqnonlin(...)[x,resnorm,residual,exitflag,output]=lsqnonlin(...)[x,resnorm,residual,exitflag,output,lambda]=lsqnonlin(...)[x,resnorm,residual,exitflag,output,lambda,jacobian]=lsqnonlin(...)resnorm——返回beta处的残差平方和，sum((fun(beta,xdata)-ydata).^2)residual——返回解beta处的残差，

fun(beta,xdata)-ydataexitflag——退出方式output——返回优化信息的输出结果，iterations、funcCount、algorithm、stepsizelambda——解beta处的Lagrange乘子jacobian——返回函数在解beta处的Jacobian矩阵36lsqnonlin的用于拟合的目标函数定义拟合函数的定义：functionF=myfun(beta,xdata,ydata)F=fitfun(beta,xdata)-ydata.%Computefunctionvaluesatbetabeta=lsqnonlin(@myfun,beta0,[],xdata,ydata)Note:funshouldreturnthesum-of-squaressum((fun(beta,xdata)-ydata).^2).ForexamplefunctionF=myfun(beta,xdata,ydata) F=beta(1)*exp(beta(2)*xdata)-ydata;37（四）nlinfit的使用nlinfit()是lsqcurvefit()和lsqnonlin()的简化版本调用形式beta=nlinfit(x,y,fun,beta0)[beta,residual,jacobian]=nlinfit(x,y,fun,beta0)[...]=nlinfit(x,y,fun,beta0,options)x和y——分别为观察节点和观察值向量，行数要相同；fun——自定义的非线性拟合模型；beta0——拟合参数的初始值；beta——拟合模型中的参数；residual——残差；jacobian——Jacobian矩阵38（五）nlparci的使用Confidenceintervalsforparametersinnonlinearregression调用ci=nlparci(beta,'covar',sigma)ci=nlparci(beta,resid,'jacobian',J)returnsthe95%confidenceintervalsciforthenonlinearleastsquaresparameterestimatesbeta.ci=nlparci(...,'alpha',alpha)

returns100(1-alpha)%confidenceintervals.Forexampleloadreaction[beta,residual,jacobian]=nlinfit(reactants,rate,'hougen',beta);ci=nlparci(beta,residual,jacobian)ci=-0.74673.2519-0.03770.1632-0.03120.1113-0.06090.2857-0.73813.1208

LowlimitandUpperlimit39（六）nlpredci的使用Confidenceintervalsforpredictionsinnonlinearregression调用[ypred,delta]=nlpredci(modelfun,x,beta,resid,'covar',sigma)[ypred,delta]=nlpredci(modelfun,x,beta,resid,'jacobian',J)returnspredictions,ypred,and95%confidenceintervalhalf-widths,delta,forthenonlinearregressionmodeldefinedbymodelfun,atinputvaluesx.[...]=nlpredci(...,'param1',val1,'param2',val2,...)Forexampleloadreaction;[beta,resid,J]=nlinfit(reactants,rate,@hougen,beta);newX=reactants(1:2,:);[ypred,delta]=nlpredci(@hougen,newX,beta,resid,J);ypred=8.41793.9542delta=0.28050.2474

ypred±delta40使用MATLAB进行参数估计示例例1：等温积分反应器的参数估计 在一等温积分反应器的动力学实验中，发生如下反应：已知：当t=0时，CA0=1，CB0=0，CC0=0。实验数据如下表：ti051015205080100200CAi,exp1.0000.8400.6790.5030.4250.1370.0300.0060.000CBi,exp0.0000.1550.3270.4260.4770.5710.5020.3850.145模型：因为CA+CB+CC=1，所以描述反应器出口反应物浓度变化的模型方程只需要两个：初始条件为：t=0，CA0=1，CB0=0，CC0=0。41使用MATLAB的dsolve求积分>>s=dsolve('Dy1=-k1*y1,Dy2=k1*y1-k2*y2','y1(0)=1,y2(0)=0')

s=y1:[1x1sym]

y2:[1x1sym]

>>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)42对积分式进行参数估计的源程序functionseqcurvefit11clearall;loadseqdata;beta0=[0.0050.001];

lb=[00];ub=[infinf];[beta,resnorm,residual,exitflag,output,lambda,jacobian]=...lsqnonlin(@seqfun,beta0,lb,ub,[],t,c);

ci=nlparci(beta,residual,jacobian);

functiony=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)];

43输出计算结果的源程序%printresultfprintf('\nEstimatedParametersbyLsqnonlin():\n')fprintf('\tk1=%.4f±

%.4f\n',beta(1),ci(1,2)-beta(1))fprintf('\tk2=%.4f±

%.4f\n',beta(2),ci(2,2)-beta(2))fprintf('Thesumofthesquaresis:%.1e\n\n',sum(residual.^2))%plotoffitresultstc=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-');holdonplot(t,c(:,2),'b+',tc,yc(:,2),'b-');xlabel('Time');ylabel('Concentration');holdoff用函数求拟合值。注意转置，与函数定义的矩阵维数一致！44参数拟合结果>>seqcurvefit11Optimizationterminated:relativefunctionvaluechangingbylessthanOPTIONS.TolFun.EstimatedParametersbyLsqnonlin(): k1=0.0412±0.0018 k2=0.0121±0.0008Thesumoftheresidualsquaresis:2.7e-00345方法2：对微分式进行参数估计functionseqcurvefit21clearall;loadseqdata;[y_row,y_col]=size(c);

beta0=[0.0050.001];c0=[10];lb=[00];ub=[infinf];[beta,resnorm,residual,exitflag,output,lambda,jacobian]=...lsqnonlin(@seqfun,beta0,lb,ub,[],t,c,y_col,c0);ci=nlparci(beta,residual,jacobian);

46拟合模型和目标函数的定义functiony=seqfun(beta,t,c,y_col,c0)%Objectivefunctiontspan=[0max(t)];[ttyy]=ode45(@modeleqs,tspan,c0,[],beta);forcol=1:y_colyc(:,col)=spline(tt,yy(:,col),t);endy=[c(:,1)-yc(:,1);c(:,2)-yc(:,2)];

functiondydt=modeleqs(t,y,beta)%Modelequationdydt=[-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);47输出计算结果的源程序%printresultfprintf('\nEstimatedParametersbyLsqnonlin():\n')fprintf('\tk1=%.4f±

%.4f\n',beta(1),ci(1,2)-beta(1))fprintf('\tk2=%.4f±

%.4f\n',beta(2),ci(2,2)-beta(2))fprintf('Thesumofthesquaresis:%.1e\n\n',sum(residual.^2))%plotoffitresultstspan=[0max(t)];[ttyc]=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-');holdonycb=spline(tt,yc(:,2),tc);plot(t,c(:,2),'b+',tc,ycb,'b-');xlabel('Time');ylabel('Concentration');holdoff

首先解微分方程求拟合值。然后用样条插值逼近。48参数拟合结果>>seqcurvefit21Optimizationterminated:relativefunctionvaluechangingbylessthanOPTIONS.TolFun.EstimatedParametersbyLsqnonlin(): k1=0.0412±0.0018 k2=0.0121±0.0008Thesumoftheresidualsquaresis:5.0e-003

49使用MATLAB进行参数估计示例例2：青霉素发酵过程动力学参数估计：在间歇发酵罐中研究青霉素发酵过程动力学，微生物Penicilliumchrysogenum在一定的控制条件下生长，细胞生长速率可以用逻辑模型描述：青霉素的生产速率模型为：

初始条件为：t=0，y1=0.29，y2=0

。TimehoursCellconcentrationdryweightPenicillinconcentrationunits/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.654250使用MATLAB的dsolve求积分>>s=dsolve('Dy1=k1*y1*(1-y1/k1),Dy2=k3*y1-k4*y2','