多元模型回归与分析.ppt

多元数据模型回归与分析,2,一、实验数据分析,由实验数据回归模型，得到模型参数前，对数据自变量间的线性相关性进行检验，是发现回归模型应用的可靠性和准确性受限制的有效方法。 因自变量间的线性相关性，使得无法区分它们对因变量的作用; 回归模型参数时会遇到几乎是奇异的数据矩阵，这样的模型参数有很大的不确定性(95的参数置信度范围宽)。 例：回归二氧化硫的催化氧化速率方程：,装有载铂氧化铝催化剂颗粒的微分固定床反应器中，测定二氧化硫的催化氧化速率。总压为790 mmHg时，记录流体相的组成分压，有下表所示的速率结果，通过这些数据求取二氧化硫的催化氧化速率方程。,3,二氧化硫的催化氧化速率,表82 二氧化硫的催化氧化速率,4,两种模型的非线性回归,1、一般的指数速率方程形式,(8.2.1),k=0.517113.3；a=-1.987.02； b=-0.2164.556；c=6.078124.7,拟合结果：,参数的95置信度太宽，模型参数不可靠。,2、根据原子氧的吸附机理，得到的速率方程式（Smith，Chemical Engineering Kinetics, 3rd Ed., 1981, McGraw-Hill, P.374）,(8.2.2),K=73，为反应平衡常数 A=0.10170.0958；B=16.024.33,拟合结果：,与方程(8.2.1)相比，方程参数的置信度有了显著改善。,5,对速率方程的进一步分析,如果把方程(8.2.2)改写为：,(8.2.3),将模型参数代入计算并以方程左边为横坐标、右边为纵坐标作图。,结果并不是斜率为-1的直线。说明表所给的速率数据没有足够的信息来表明速率方程中的逆反应贡献。,如将SO3分压对O2分压作图，这两分压间有近似线性关系。 所以方程（8.2.1）的置信区间范围大。,6,二、回归模型的选择（1）,例：水饱和蒸汽压的模型回归 水的蒸汽压数据选用的温度范围为0120,三参数的Antoine方程：,四参数的Riedel回归方程：,五参数回归方程（参考Thek-Stiel的蒸汽压预测方程提出）：,(8.2.4),(8.2.5),(8.2.6),7,水饱和蒸汽压的模型回归结果,表8-3 水饱和蒸汽压的方程拟合结果,拟合度十分接近1，表明拟合是成功的，但实际上用Antoine方程来拟合回归得到的结果不理想，说明仅从拟合度上来判断结果的好坏是不够的。为什么呢？,8,因变量与残差关系图,残差定义：,(8.2.7),考察模型参数估计方法的两个基本假设： 参数估计的误差相互不相关联，是随机的。 估计误差符合正态分布。,检查模型适合体系数据程度的最有效方法之一是对因变量与残差作图，观察其分布情况。,9,Antoine方程拟合的残差,残差虽然很小，但其分布不是随机的。,残差的分布同正态分布相比，有较大的差距。,两方面的结果充分说明了拟合回归的Antoine方程还不能充分反映蒸汽压与温度间的关系，造成残差间存在关联。 采用Riedel方程拟合得到的也是类似的结果。,10,改进Thek-Stiel方程方程的拟合结果,拟合误差比Antoine方程小了近一个数量级，而且残差分布是随机分布的。,误差分布基本符合正态分布。,改进Thek-Stiel方程方程描述水饱和蒸汽压的合适模型。,11,二、回归模型的选择（2）,前面说明了模型参数较少时会出现拟合残差的分布不是随机的，而是呈现某种分布，相互关联。 在模型回归拟合数据的过程中，如模型参数过多会出现什么情况？如何判断回归拟合模型中有过多的参数呢？,12,丙烷在氢型丝光沸石上的吸附平衡,例：选用不同吸附方程拟合丙烷在氢型丝光沸石体系303K的吸附平衡数据。 目标：说明如何对模型拟合结果进行统计分析，确定模型拟合的好坏、模型参数的可靠性和准确性，从而进行拟合模型的选择。,表8-4 303K时丙烷在氢型丝光沸石上的吸附平衡数据,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,各模型的计算结果,表85 吸附等温线关联的参数值、方差和回归系数,从表中可看出，方程(814)拟合方差逐渐减少，回归系数更接近1（方程(12)是通过压力数据来拟合的，故拟合方差和其它方程的结果不是在同一数量级上)。由方程(14)的五参数形式改进的方程(15)式获得的结果最好，实验数据点几乎完全落在方程（15）式的曲线上(见下图)。,15,方程(13)和方程(15)的拟合结果,方程(15)式获得的结果最好，实验数据点几乎完全落在方程（15）式的曲线上。,16,判断模型参数是否过少的依据,通过对方程(13)和五参数方程(15)的残差进行分析， 方程(13)因参数过少，吸附量的计算误差与实验吸附量之间存在着某种分布。 方程(15) 计算误差在零的两边是随机分布的，看不出规律性。因此，拟合计算误差有无规律性的分布是判断模型参数是否过少的依据。,因此，拟合计算误差有无规律性的分布是判断模型参数是否过少的依据。,方程(13)的拟合误差,方程(15)的拟合误差,17,判断模型参数是否过多的依据,在方程(14)的计算结果中，有些参数95%的置信度较大，说明这些参数之间有联系，不是独立的。 而对于方程(14)的五参数形式，即方程(15)，其所有参数的95%置信度都较小。事实上，方程(15)就是据此分析对吸附平衡理论作进一步研究而获得的。,因此，拟合参数95%的置信度是否较大是判断模型参数是否过多的依据。,18,回归模型的选择总结,模型参数较少时会出现拟合残差的分布不是随机的，而是呈现某种分布，相互关联。残差的分布偏离正态分布较远。 模型参数过多会出现某些参数95%的置信度较大，说明这些参数之间有联系，不是独立的。,19,习题,研究二硫化碳饱和蒸汽压的模型回归问题。（P266，Ex8.3）,二硫化碳的基本性质： 临界温度为273.05 临界压力为72.87 atm。,20,Statistica的非线性估计,非线性估计方法 User-Specified Regression, least square,可以计算95%置信区间,21,“least square”与“Custom Loss ”比较,“least square”计算结果（采用Levenberg-Marquardt方法）,“Custom Loss ”计算结果（采用Quasi-Newton法） Matrix ill conditioned; cannot compute standard errors.,EConf,22,“Custom Loss ”的方差分析与迭代步骤,“Custom Loss ”无迭代历史纪录，协方差分析结果已出现病态。,Covariance matrix cannot be computed.,23,Statistica非线性估计的残差分析,残差分布情况,残差对预测值作图,24,对方程(15)的残差分析,Histogram of residuals,Residual vs. Predicted,25,误差正态分布图,“least square”计算的误差正态分布图,“Custom Loss ”计算的误差正态分布图,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表示插值多项式的系数。 例：某实验中测得一组数据，其值如下： 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.69,29,（二）lsqcurvefit的使用,方程（目标函数） Find coefficients beta that best fit the equation,调用形式 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。,30,lsqcurvefit的使用（续）,调用形式（续） 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矩阵,31,lsqcurvefit的拟合函数的定义,拟合函数的定义：,function F = myfun(beta,xdata) F = . % Compute function values at x,beta = lsqcurvefit(myfun,beta0,xdata,ydata),Note: (1)fun should return fun(x,xdata), and not the sum-of-squares sum(fun(x,xdata)-ydata).2). (2)The algorithm implicitly squares and sums fun(x,xdata)-ydata.,32,lsqcurvefit的拟合函数的定义（续）,If the Jacobian can also be computed by user-defined,function F,J = myfun(beta,xdata) F = . % objective function values at beta if nargout 1 % two output arguments J = . % Jacobian of the function evaluated at beta end,options = optimset(Jacobian,on),For example,function F = myfun(beta,xdata) F = beta(1)*exp(beta(2)*xdata);,33,lsqcurvefit应用示例,% Assume you determined xdata and ydata experimentally xdata = 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 guess beta,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.5049,34,（三）lsqnonlin的使用,方程（目标函数） Find coefficients beta that best fit the equation,调用形式 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 = ub。,35,lsqnonlin的使用（续）,调用形式（续） 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)-ydata exitflag退出方式 output返回优化信息的输出结果，iterations、funcCount、algorithm、stepsize lambda解beta处的Lagrange乘子 jacobian返回函数在解beta处的Jacobian矩阵,36,lsqnonlin的用于拟合的目标函数定义,拟合函数的定义：,function F = myfun(beta,xdata,ydata) F = fitfun(beta,xdata)-ydata . % Compute function values at beta,beta = lsqnonlin(myfun, beta0,xdata,ydata),Note: fun should return the sum-of-squares sum(fun(beta,xdata)-ydata).2).,For example,function F = 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残差； jacobianJacobian矩阵,38,（五）nlparci的使用,Confidence intervals for parameters in nonlinear regression 调用 ci = nlparci(beta,covar,sigma) ci = nlparci(beta,resid,jacobian,J) returns the 95% confidence intervals ci for the nonlinear least squares parameter estimates beta. ci = nlparci(.,alpha,alpha) returns 100(1-alpha)% confidence intervals. For example,load reaction beta, residual, jacobian = nlinfit(reactants,rate,hougen,beta); ci = nlparci(beta, residual, jacobian),ci = -0.7467 3.2519 -0.0377 0.1632 -0.0312 0.1113 -0.0609 0.2857 -0.7381 3.1208,Low limit and Upper limit,39,（六）nlpredci的使用,Confidence intervals for predictions in nonlinear regression 调用 ypred, delta = nlpredci(modelfun,x,beta,resid,covar,sigma) ypred, delta = nlpredci(modelfun,x,beta,resid,jacobian,J) returns predictions, ypred, and 95% confidence interval half-widths, delta, for the nonlinear regression model defined by modelfun, at input values x. . = nlpredci(.,param1,val1,param2,val2,.) For example,load reaction; beta,resid,J = nlinfit(reactants,rate,hougen,beta); newX = reactants(1:2,:); ypred, delta = nlpredci(hougen,newX,beta,resid,J);,ypred = 8.4179 3.9542 delta = 0.2805 0.2474,ypreddelta,40,使用MATLAB进行参数估计示例,例1：等温积分反应器的参数估计 在一等温积分反应器的动力学实验中，发生如下反应 ： 已知：当t=0时，CA0=1，CB0=0，CC0=0。实验数据如下表：,模型： 因为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: 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),42,对积分式进行参数估计的源程序,function seqcurvefit11 clear 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);,43,输出计算结果的源程序,% print result fprintf(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 results tc=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 on plot(t,c(:,2),b+,tc,yc(:,2),b-); xlabel(Time); ylabel(Concentration); hold off,用函数求拟合值。 注意转置，与函数定义的矩阵维数一致！,44,参数拟合结果, seqcurvefit11 Optimization 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-003,45,方法2：对微分式进行参数估计,function seqcurvefit21 clear 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);,46,拟合模型和目标函数的定义,function y = seqfun(beta,t,c,y_col,c0) % Objective function tspan = 0 max(t); tt yy = ode45(modeleqs,tspan,c0,beta); for col = 1:y_col yc(:,col) = spline(tt,yy(:,col),t); end y=c(:,1)-yc(:,1); c(:,2)-yc(:,2); function dydt = modeleqs(t,y,beta) % Model equation dydt = -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,输出计算结果的源程序,% print result fprintf(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 results tspan = 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 on ycb = spline(tt,yc(:,2),tc); plot(t,c(:,2),b+,tc,ycb,b-); xlabel(Time); ylabel(Concentration); hold off,首先解微分方程求拟合值。然后用样条插值逼近。,48,参数拟合结果, seqcurvefit21 Optimization 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,49,使用MATLAB进行参数估计示例,例2：青霉素发酵过程动力学参数估计： 在间歇发酵罐中研究青霉素发酵过程动力学，微生物Penicillium chrysogenum在一定的控制条件下生长，细胞生长速率可以用逻辑模型描述： 青霉素的生产速率模型为：,初始条件为：t=0， y1=0.29，y2=0 。,50,使用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),没有确定的积分表达式,51,参数估计的源程序,function PenicilliumEst clear 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);,52,拟合模型和目标函数的定义,% = function f = Func(beta,x,y,y_col,y0) % Define object