最小二乘法在数学模型建立与检验中的运用.docVIP

最小二乘法在数学模型建立与检验中的应用信息与计算科学专业2008级 周建勤摘要：本文主要研究了最小二乘法在建立数学模型中的参数学模型中的参数估计数估计，模型检验中的应用。通过给出最小二乘法在Matlab中的代码计算模型参数，误差精确度，并给出检验模型是否具有多重共线，异方差性，序列相关性方法。关键词：最小二乘法；参数估计；误差精确度；多重共线性；异方差；自相关。Application of Least-Square Method on establish and test mathematical modelZhou Jianqin ,Grade 2008,Information and Computing ScienceAbstract: In this text we main consider application of Least-Square Method in use of parameter estimation and model checking in mathematical models. By giving the least squares methods code in Matlab to find the model Heteroscedasticity, autocorrelation method parameters,Error accuracy and Test whether the model with multiple collinear heteroscedasticity, autocorrelation method.Keywords: least squares method; parameter estimation; error accuracy; multicollinearity; heteroscedasticity; autocorrelation.背景介绍 最小二乘方法最早是有高斯提出的，他用这种方法解决了天文学方面的问题，特别是确定了某些行星和彗星的天体轨迹。最小二乘法普遍适用于各个科学领域，它在解决实际问题中发挥了重要的作用。它在生产实践、科学实验及经济活动中均有广泛应用。 近半个多世纪以来，随着计算机技术的迅速发展，数学的应用不仅在工程技术、自然科学等领域发挥着越来越重要的作用，而且以空前的广度和深度向经济、管理、金融、生物、医学、环境、地质、人口、交通等新的领域渗透，所谓数学技术已经成为当代高新技术的重要组成部分 。数学模型是一种模拟，是用数学符号、数学式子、程序、图形等对实际课题本质属性的抽象而又简洁的刻划，它或能解释某些客观现象，或能预测未来的发展规律，或能为控制某一现象的发展提供某种意义下的最优策略或较好策略。数学模型一般并非现实问题的直接翻版，它的建立常常既需要人们对现实问题深入细微的观察和分析，又需要人们灵活巧妙地利用各种数学知识。这种应用知识从实际课题中抽象、提炼出数学模型的过程就称为数学建模。不论是用数学方法在科技和生产领域解决哪类实际问题，还是与其它学科相结合形成交叉学科，首要的和关键的一步是建立研究对象的数学模型，并加以计算求解。 在已知系统模型结构时，用系统的输入和输出数据计算系统模型参数的过程。18世纪末德国数学家C.F.高斯首先提出参数估计的方法，他用最小二乘法计算天体运行的轨道。20世纪60年代,随着电子计算机的普及,参数估计有了飞速的发展。参数估计有多种方法，有最小二乘法、极大似然法、极大验后法、最小风险法和极小化极大熵法等。在一定条件下，后面三个方法都与极大似然法相同。最基本的方法是最小二乘法和极大似然法。 由于测量仪器的精度不完善和人为因素及外界条件的影响，测量误差总是不可避 免的。为了提高成果的质量，处理好这些测量中存在的误差问题，观测值的个数往往要多于确定未知量所必须观测的个数，也就是要进行多余观测。有了多余观测，势必在观测结果之间产生矛盾，测量平差的目的就在于消除这些矛盾而求得观测量的最可靠结果并评定测量成果的精度。测量平差采用的原理就是“最小二乘法”。 由于误差的客观存在，真值一般是无法测得的。测量次数无限多时，根据正负误差出现的概率相等的误差分布定律，在不存在系统误差的情况下，它们的平均值极为接近真值。故在实验科学中真值的定义为无限多次观测值的平均值。但实际测定的次数总是有限的，由有限次数求出的平均值，只能近似地接近于真值，可称此平均值为最佳值。 所谓多重共线性（Multicollinearity）是指线性回归模型中的解释变量之间由于存在精确相关关系或高度相关关系而使模型估计失真或难以估计准确。一般来说，由于经济数据的限制使得模型设计不当，导致设计矩阵中解释变量间存在普遍的相关关系。 经济学是非实验型科学，经济数据是被动生成和由从事经济研究的人员被动获得，而且经济数据的获得是不可控的，大多数情况下，人们并不能按照自己的设计与要求获得相应的经济数据。所以，为建模研究而取得的样本数据常常不能提供足够的信息，以至于导致多重共线性的产生。异方差性（heteroscedasticity ）是为了保证回归参数估计量具有良好的统计性质，经典线性回归模型的一个重要假定是：总体回归函数中的随机误差项满足同方差性，即它们都有相同的方差。如果这一假定不满足，则称线性回归模型存在异方差性。 线性回归模型中随机误差项存在序列相关的原因很多，但主要是经济变量自身特点、数据特点、变量选择及模型函数形式选择引起的。1.经济变量惯性的作用引起随机误差项自相关2. 经济行为的滞后性引起随机误差项自相关3. 一些随机因素的干扰或影响引起随机误差项自相关4. 模型设定误差引起随机误差项自相关5.观测数据处理引起随机误差项序列相关 文章对经典的最小二乘方法的应用背景、原理与算法进行了介绍，给出了它们在线性模型参数估计中的MATLAB实现，以经典单方程数学模型为对象，介绍建立数学模型的过程，主要采用回归分析的方法的数学模型，采用最小二乘法计算模型中参数，误差估计，并检验模型的多重共线性，异方差，自相关性。知识预备最小二乘法最小二乘法（又称最小平方法）是一种数学优化技术。它通过最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以简便地求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方和为最小。最小二乘法还可用于曲线拟合。其他一些优化问题也可通过最小化能量或最大化熵用最小二乘法来表达 。最小二乘法原理在我们研究两个变量（x,y）之间的相互关系时，通常可以得到一系列成对的数据（x1,y1.x2,y2. xm,ym）；将这些数据描绘在x -y直角坐标系中，若发现这些点在一条直线附近，可以令这条直线方程如（式1-1）。Y计= a0 + a1 X （式1-1)其中：a0、a1 是任意实数为建立这直线方程就要确定a0和a1，应用最小二乘法原理，将实测值Yi与利用（式1-1）计算值（Y计=a0+a1X）的离差（Yi-Y计）的平方和（Yi - Y计）2最小为“优化判据”。令： = （Yi - Y计）（式1-2)把（式1-1）代入（式1-2）中得： = （Yi - a0 - a1 Xi)2 （式1-3)当（Yi-Y计）平方最小时，可用函数 对a0、a1求偏导数，令这两个偏导数等于零。（式1-4)（式1-5)亦即：m a0 + （Xi ) a1 = Yi （式1-6)（Xi ) a0 + （Xi2 ) a1 = （Xi,Yi) （式1-7)得到的两个关于a0、 a1为未知数的两个方程组，解这两个方程组得出：a0 = （Yi) / m - a1（Xi) / m （式1-8)a1 = mXi Yi - （Xi Yi) / mXi2 - （Xi)2 ) （式1-9)这时把a0、a1代入（式1-1）中， 此时的(式1-1）就是我们回归的元线性方程即：数学模型。在回归过程中，回归的关联式是不可能全部通过每个回归数据点（x1,y1. x2,y2.xm,ym），为了判断关联式的好坏，可借助相关系数“R”，统计量“F”，剩余标准偏差“S”进行判断；“R”越趋近于 1 越好；“F”的绝对值越大越好；“S”越趋近于 0 越好。R = XiYi - m （Xi / m）（Yi / m)/ SQRXi2 - m （Xi / m)2Yi2 - m （Yi / m)2 （式1-10) *在（式1-1）中，m为样本容量，即实验次数；Xi、Yi分别任意一组实验X、Y的数值。编辑本段最小二乘法公式最小二乘法公式注：以下“平”是指某参数的算数平均值。如：X平x的算术平均值。1、（X-X平）（Y-Y平）=（XY-X平Y-XY平+X平Y平）=XY-X平Y-Y平X+nX平Y平=XY-nX平Y平-nX平Y平+nX平Y平=XY-nX平Y平；2、（X -X平）2=（X2-2XX平+X平2)=X2-2nX平2+nX平2=X2-nX平2；3、Y=kX+bk=（XY）平-X平*Y平）/（X2）平-(X平）2），b=Y平-kX平；X平=1/nXi，(XY）平=1/nXiYi；编辑本段最小二乘法拟合对给定数据点(Xi，Yi)(i=0,1，m），在取定的函数类 中，求p(x），使误差的平方和E2最小，E2=p(Xi)-Yi2。从几何意义上讲，就是寻求与给定点 (Xi，Yi)(i=0,1，m）的距离平方和为最小的曲线y=p(x）。函数p(x）称为拟合函数或最小二乘解，求拟合函数p(x）的方法称为曲线拟合的最小二乘法。最小二乘法的矩阵形式Ax=b，其中A为nxk的矩阵，x为kx1的列向量，b为nx1的列向量，nk。这个方程系统称为Over Determined System，如果n x=-2.5 -1.7 -1.1 -0.8 0 0.1 1.5 2.7 3.6; y=-192.9 -85.50 -36.15 -26.52 -9.10 -8.43 -13.12 6.50 68.04;plot(x,y,r*),legend(实验数据(xi,yi)xlabel(x), ylabel(y),title(例7.2.1的数据点(xi,yi)的散点图)运行后屏幕显示数据的散点图（略）. （3）编写下列MATLAB程序计算在处的函数值，即输入程序 syms a1 a2 a3 a4x=-2.5 -1.7 -1.1 -0.8 0 0.1 1.5 2.7 3.6; fi=a1.*x.3+ a2.*x.2+ a3.*x+ a4运行后屏幕显示关于a1,a2, a3和a4的线性方程组fi = -125/8*a1+25/4*a2-5/2*a3+a4, -4913/1000*a1+289/100*a2-17/10*a3+a4, -1331/1000*a1+121/100*a2-11/10*a3+a4, -64/125*a1+16/25*a2-4/5*a3+a4, a4, 1/1000*a1+1/100*a2+1/10*a3+a4, 27/8*a1+9/4*a2+3/2*a3+a4, 19683/1000*a1+729/100*a2+27/10*a3+a4, 5832/125*a1+324/25*a2+18/5*a3+a4编写构造误差平方和的MATLAB程序 y=-192.9 -85.50 -36.15 -26.52 -9.10 -8.43 -13.12 6.50 68.04;fi=-125/8*a1+25/4*a2-5/2*a3+a4, -4913/1000*a1+289/100*a2-17/10*a3+a4, -1331/1000*a1+121/100*a2-11/10*a3+a4, -64/125*a1+16/25*a2-4/5*a3+a4, a4, 1/1000*a1+1/100*a2+1/10*a3+a4, 27/8*a1+9/4*a2+3/2*a3+a4, 19683/1000*a1+729/100*a2+27/10*a3+a4, 5832/125*a1+324/25*a2+18/5*a3+a4;fy=fi-y; fy2=fy.2; J=sum(fy.2)运行后屏幕显示误差平方和如下J= (-125/8*a1+25/4*a2-5/2*a3+a4+1929/10)2+(-4913/1000*a1+289/100*a2-17/10*a3+a4+171/2)2+(-1331/1000*a1+121/100*a2-11/10*a3+a4+723/20)2+(-64/125*a1+16/25*a2-4/5*a3+a4+663/25)2+(a4+91/10)2+(1/1000*a1+1/100*a2+1/10*a3+a4+843/100)2+(27/8*a1+9/4*a2+3/2*a3+a4+328/25)2+(19683/1000*a1+729/100*a2+27/10*a3+a4-13/2)2+(5832/125*a1+324/25*a2+18/5*a3+a4-1701/25)2为求使达到最小，只需利用极值的必要条件 ，得到关于的线性方程组，这可以由下面的MATLAB程序完成，即输入程序 syms a1 a2 a3 a4J=(-125/8*a1+25/4*a2-5/2*a3+a4+1929/10)2+(-4913/1000*a1+289/100*a2-17/10*a3+a4.+171/2)2+(-1331/1000*a1+121/100*a2-11/10*a3+a4+723/20)2+(-64/125*a1+16/25*a2-4/5*a3+a4+663/25)2+(a4+91/10)2+(1/1000*a1+1/100*a2+1/10*a3+a4+843/100)2+(27/8*a1+9/4*a2+3/2*a3+a4+328/25)2+(19683/1000*a1+729/100*a2+27/10*a3+a4-13/2)2+(5832/125*a1+324/25*a2+18/5*a3+a4-1701/25)2; Ja1=diff(J,a1); Ja2=diff(J,a2); Ja3=diff(J,a3); Ja4=diff(J,a4);Ja11=simple(Ja1), Ja21=simple(Ja2), Ja31=simple(Ja3), Ja41=simple(Ja4),运行后屏幕显示J分别对a1, a2 ,a3 ,a4的偏导数如下Ja11=56918107/10000*a1+32097579/25000*a2+1377283/2500*a3+23667/250*a4-8442429/625Ja21 = 32097579/25000*a1+1377283/2500*a2+23667/250*a3+67*a4+767319/625Ja31 =1377283/2500*a1+23667/250*a2+67*a3+18/5*a4-232638/125Ja41 =23667/250*a1+67*a2+18/5*a3+18*a4+14859/25解线性方程组Ja11 =0，Ja21 =0，Ja31 =0，Ja41 =0，输入下列程序A=56918107/10000, 32097579/25000, 1377283/2500, 23667/250; 32097579/25000, 1377283/2500, 23667/250, 67; 1377283/2500, 23667/250, 67, 18/5; 23667/250, 67, 18/5, 18;B=8442429/625, -767319/625, 232638/125, -14859/25; C=B/A, f=poly2sym(C)运行后屏幕显示拟合函数f及其系数C如下C = 5.0911 -14.1905 6.4102 -8.2574f=716503695845759/140737488355328*x3-7988544102557579/562949953421312*x2+1804307491277693/281474976710656*x-4648521160813215/562949953421312 故所求的拟合曲线为.（4）编写下面的MATLAB程序估计其误差，并作出拟合曲线和数据的图形.输入程序 xi=-2.5 -1.7 -1.1 -0.8 0 0.1 1.5 2.7 3.6; y=-192.9 -85.50 -36.15 -26.52 -9.10 -8.43 -13.12 6.50 68.04;n=length(xi); f=5.0911.*xi.3-14.1905.*xi.2+6.4102.*xi -8.2574;x=-2.5:0.01: 3.6; F=5.0911.*x.3-14.1905.*x.2+6.4102.*x -8.2574;fy=abs(f-y); fy2=fy.2; Ew=max(fy), E1=sum(fy)/n, E2=sqrt(sum(fy2)/n)plot(xi,y,r*), hold on, plot(x,F,b-), hold offlegend(数据点(xi,yi),拟合曲线y=f(x), xlabel(x), ylabel(y),title(例7.2.1的数据点(xi,yi)和拟合曲线y=f(x)的图形)运行后屏幕显示数据与拟合函数f的最大误差Ew，平均误差E1和均方根误差E2及其数据点和拟合曲线y=f(x)的图形（略）.Ew = E1 = E2 =3.105 4 0.903 4 1.240 9多重共线性PLS偏最小二乘法经典的最小二乘估计，必需满足一些假设条件，多重共线性就是其中的一种。实际上，解释变量间完全不相关的情形是非常少见的，大多数变量都在某种程度上存在着一定的共线性，而存在着共线性会给模型带来许多不确定性的结果。设回归模型如果矩阵的列向量存在一组不全为零的数使得, =1,2,则称其存在完全共线性,如果, =1,2,则称其存在近似的多重共线性。(偏最小二乘法) H.Wold在1975年提出的 偏最小二乘法近年来引起广泛的关注，在解决多重共线性方面能很好的达到目的，偏最小二乘法集中了最小二乘法、主成分分析法和典型相关分析的的优点克服了两种方法的缺点。偏最小二乘法吸取了主成分回归提取主成分的思想，但不同的是主成分回归只是从自变量中去寻找主成分与因变量无关，因而主成分与因变量在算法上关系不密切，从而导致最后主成分在实际应用中无法更好的进一步拟合因变量，偏最小二乘法则是从因变量出发，选择与因变量相关性较强而又能方便运算的自变量的线性组合。 偏最小二乘回归分析偏最小二乘回归分析考虑p个因变量L与m个自变量的建模问题。偏最小二乘回归的基本作法是首先在自变量集中提出第一成分(是的线性组合，且尽可能多地提取原自变量集中的变异信息）；同时在因变量集中也提取第一成分并要求与相关程度达到最大。然后建立因变量与的回归，如果回归方程已达到满意的精度，则算法中止。否则继续第二对成分的提取，直到能达到满意的精度为止。若最终对自变量集提取r个成分，偏最小二乘回归将通过建立与的回归式，然后再表示为与原自变量的回归方程式，即偏最小二乘回归方程式。偏最小二乘回归MATLAB程序代码单因变量function y=pls(pz)row,col=size(pz);aver=mean(pz);stdcov=std(pz); %求均值和标准差rr=corrcoef(pz); %求相关系数矩阵%data=zscore(pz); %数据标准化stdarr = ( pz - aver(ones(row,1),:) )./ stdcov( ones(row,1),:); % 标准化数据结果与zscore()一致x0=pz(:,1:col-1);y0=pz(:,end); %提取原始的自变量、因变量数据e0=stdarr(:,1:col-1);f0=stdarr(:,end); %提取标准化后的自变量、因变量数据num=size(e0,1);%求样本点的个数temp=eye(col-1);%对角阵for i=1:col-1%以下计算 w，w*和 t 的得分向量， w(:,i)= ( e0* f0 )/ norm( e0*f0 ); t(:,i)=e0*w(:,i) %计算成分 ti 的得分 alpha(:,i)=e0*t(:,i)/(t(:,i)*t(:,i) %计算 alpha_i ,其中(t(:,i)*t(:,i)等价于norm(t(:,i)2 e=e0-t(:,i)*alpha(:,i) %计算残差矩阵 e0=e; %计算w*矩阵 if i=1 w_star(:,i)=w(:,i); else for j=1:i-1 temp=temp*(eye(col-1)-w(:,j)*alpha(:,j); end w_star(:,i)=temp*w(:,i); end%以下计算 ss(i)的值 beta=t(:,1:i),ones(num,1)f0 %求回归方程的系数 beta(end,:)=; %删除回归分析的常数项 cancha=f0-t(:,1:i)*beta; %求残差矩阵 ss(i)=sum(sum(cancha.2); %求误差平方和%以下计算 press(i) for j=1:num t1=t(:,1:i);f1=f0; she_t=t1(j,:);she_f=f1(j,:); %把舍去的第 j个样本点保存起来 t1(j,:)=;f1(j,:)=; %删除第j个观测值 beta1=t1,ones(num-1,1)f1; %求回归分析的系数