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

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

信息与计算科学专业2008级 周建勤

摘要：本文主要研究了最小二乘法在建立数学模型中的参数学模型中的参数估计数估计，模型检验中的应用。通过给出最小二乘法在Matlab 中的代码计算模型参数，误差精确度，并给出检验模型是否具有多重共线，异方差性，序列相关性方法。

关键词：最小二乘法；参数估计；误差精确度；多重共线性；异方差；自相关。 Application of Least-Square Method on establish and test mathematical model

Zhou Jianqin ,Grade 2008,Information and Computing Science

Abstract: In this text we main consider application of Least-Square Method in use of p arameter estimation and model checking in mathematical models. By giving the least squares method's code in Matlab to find the m odel 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 = [m∑Xi Yi - （∑Xi ∑Yi)] / [m∑Xi2 - （∑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)]/ SQR{[∑Xi2 - m （∑Xi / m)2][∑Yi2 - 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=

∑（X^2--2XX平+X平^2)=

∑X^2--2nX平^2+nX平^2=∑X^2--nX平^2；

3、Y=kX+b k=（（XY ）平--X 平*Y平）/（（X^2）平--(X平）^2），

b=Y平--kX 平； X 平=1/n∑Xi ， (XY）平=1/n∑XiYi ；编辑本段最小二乘法拟合

对给定数据点{(Xi，Yi)}(i=0,1，…，m ），在取定的函数类Φ 中，求p(x）∈Φ，使误差的平方和E^2最小，E^2=∑[p(Xi)-Yi]^2。从几何意义上讲，就是寻求与给定点 {(Xi，Yi)}(i=0,1，…，m ）的距离平方和为最小的曲线y=p(x）。函数p(x）称为拟合函数或最小二乘解，求拟合函数p(x）的方法称为曲线拟合的最小二乘法。最小二乘法的矩阵形式

Ax=b，其中A 为nxk 的矩阵，x 为kx1的列向量，b 为nx1的列向量，n>k。这个方程系统称为Over Determined System，如果n

数学模型

模型准备

了解问题的实际背景，明确其实际意义，掌握对象的各种信息。用数学语言来描述问题。 模型假设

根据实际对象的特征和建模的目的，对问题进行必要的简化，并用精确的语言提出一些恰当的假设。

模型建立

在假设的基础上，利用适当的数学工具来刻划各变量之间的数学关系，建立相应的数学结构（尽量用简单的数学工具）。

模型求解

利用获取的数据资料，对模型的所有参数做出计算（或近似计算）。

模型分析

对所得的结果进行数学上的分析。

模型检验

将模型分析结果与实际情形进行比较，以此来验证模型的准确性、合理性和适用性。如果模型与实际较吻合，则要对计算结果给出其实际含义，并进行解释。如果模型与实际吻合较差，则应该修改假设，再次重复建模过程。

模型应用

应用方式因问题的性质和建模的目的而异。

多重共线 （1）经济变量相关的共同趋势

（2）滞后变量的引入

（3）样本资料的限制

多重共线性的影响

（1）完全共线性下参数估计量不存在

（2）近似共线性下OLS 估计量非有效 多重共线性使参数估计值的方差增大，1/(1-r2)为(Variance Inflation Factor, VIF)

（3）参数估计量经济含义不合理

（4）变量的显著性检验失去意义，可能将重要的解释变量排除在模型之外

（5）模型的预测功能失效。变大的方差容易使区间预测的“区间”变大，使预测失去意义。 需要注意：即使出现较高程度的多重共线性，OLS 估计量仍具有线性性等良好的统计性质。但是OLS 法在统计推断上无法给出真正有用的信息。

异方差性

回归模型的随机扰动项ui 在不同的观测值中的方差不等于一个常数，Var(ui)= 常数（i=1，2，…，n ），或者Var （u ） Var （u ）（i j ），这时我们就称随机扰动项ui 具有异方差性（Heteroskedasticity ）。 在实际经济问题中，随机扰动项ui 往往是异方差的，但主要在截面数据分析中出现。 例如 （1）调查不同规模公司的利润，发现大公司的利润波动幅度比小公司的利润波动幅度大； （2）分析家庭支出时发现高收入家庭支出变化比低收入家庭支出变化大。 在分析家庭支出模型时，我们会发现高收入家庭通常比低收入

家庭对某些商品支出有更大的方差；图5-1显示了一元线性回归中随机变量的方差ui 随着解释变量 的增加而变化的情况。 异方差性破坏了古典模型的基本假定，如果我们直接应用最小二乘法估计回归模型，将得不到准确、有效的结果。 1．模型中缺少某些解释变量，从而随机扰动项产生系统模式 由于随机扰动项ui 包含了所有无法用解释变量表示的各种因素对被解释变量的影响，即模型中略去的经济变量对被解释变量的影响。如果其中被略去的某一因素或某些因素随着解释变量观测值的不同而对被解释变量产生不同的影响，就会使ui 产生异方差性。 例如，以某一时间截面上不同收入家庭的数据为样本，研究家庭对某一消费品（如服装、食品等）的需求，设其模型为： （5-1） 其中Qi 表示对某一消费品的需求量，Ii 为家庭收入，ui 为随机扰动项。ui 包括除家庭收入外其他因素对Qi 的影响。如：消费习惯、偏好、季节、气候等因素，ui 的方差就表示这些因素的影响可能使得Qi 偏离均值的程度。在气候异常时，高收入家庭就会拿出较多的钱来购买衣服，而低收入的家庭购买衣服的支出就很有限，这时对于不同的收入水平Ii ，Qi 偏离均值的程度是不同的，Var(ui) 常数，于是就存在异方差性了。 再比如，以某一时间截面上不同地区的数据为样本，研究某行业的产出随投入要素的变化而变化的关系，建立如下模型： （5-2） 其中Yi 表示某行业的产出水平。

Li 表示劳动力对产出的影响。Ki 表示资本对产出的影响，ui 表示除劳动力和资本外其他因素对产出水平的影响，诸如地理位置、国家政策等。显然，对于不同的行业 ，这些因素对产出 的影响程度是不 同的，引起 偏离零均值的程度也是不同的，这就出现了异方差。 异方差性容易出现在截面数据中，这是因为在截面数据中通常涉及某一确定时点上的总体单位。比如个别的消费者及其家庭、不同行业或者农村、城镇等区域的划分，这些单位各自有不同的规模或水平，一般情况下用截面数据作样本时出现异方差性的可能性较大。

自相关性

随机误差项的自相关性可以有多种形式，其中最常见的类型是随机误差项之间存在一阶自相关性或一阶自回归形式，即随机误差项只与它的前一期值相关：cov(ut,u t-1) =E（ut,u t-1) =/= 0, 或者u t=f(u t-1),则称这种关系为一阶自相关。 一阶自相关性可以表示为 ut= p1 u i-1 + p2 u i-2 + p3 u i-3 + …… p p u t-p + v t 称之为p 阶自回归形式，或模型 存在 p 阶自相关 由于无法观察到误差项 u t,只能通过残差项 e t来判断 u t 的行为。如果 u t或 e t呈出下图(a) -(d) 形式，则表示u t 存在自相关，如果 ut 或et 呈现图中 (e) 形式，

则 表示 u t不存在自相关

自相关性

线性回归模型中的随机误差项的序列相关问题较为普遍，特别是在应用时间序列资料时，随机误差项的序列相关经常发生。 自相关性产生的原因： 线性回归模型中随机误差项存在序列相关的原因很多，但主要是经济变量自身特点、数据特点、变量选择及模型函数形式选择引起的。 1. 经济变量惯性的作用引起随机误差项自相关 2. 经济行为的滞后性引起随机误差项自相关 3. 一些随机因素的干扰或影响引起随机误差项自相关

4. 模型设定误差引起随机误差项自相关 5. 观测数据处理引起随机误差项序列相关 自相关的后果： 线性相关模型的随机误差项存在自相关的情况下，用OLS （普通最小二乘法）进行参数估计，会造成以下几个方面的影响。 从高斯-马尔可夫定理的证明过

程中可以看出，只有在同方差和非自相关性的条件下，OLS 估计才具有最小方差性。当模型存在自相关性时，OLS 估计仍然是无偏估计，但不再具有有效性。这与存在异方差性时的情况一样，说明存在其他的参数估计方法，其估计误差小于OLS 估计的误差；也就是说，对于存在自相关性的模型，应该改用其他方法估计模型中的参数。 1. 自相关不影响OLS 估计量的线性和无偏性，但使之失去有效性 2. 自相关的系数估计量将有相当大的方差

3. 自相关系数的T 检验不显著 4. 模型的预测功能失效

MATLAB

MATLAB 是矩阵实验室（Matrix Laboratory）的简称，是美国MathWorks 公司出品的商业数学软件，

用于算法开发、数据可视化、

数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MATLAB 和Simulink 两大部分。

MATLAB 是由美国mathworks 公司发布的主要面对科学计算、可视化以及交互式程序设

统非交互式程序设计语言（如C 、Fortran ）的编辑模式，代表了当今国际科学计算软件的

先进水平。

matlab 开发工作界面 接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。 MATLAB 的基本数据单位是矩阵，它的指令表达式与数学、工程中常用的形式十分相似，故用MATLAB 来解算问题要比用C ，FORTRAN 等语言完成相同的事情简捷得多，并且MATLAB 也吸收了像Maple 等软件的优点，使MATLAB 成为一个强大的。在新的版本中也加入了对，，，的支持。可以直接调用, 用户也可以将自己编写的实用程序导入到MATLAB 函数库中方便自己以后调用，此外许多的MATLAB 爱好者都编写了一些经典的程序，用户可以直接进行下载就可以用。编辑本段

基本应用

MATLAB 产品族可以用来进行以下各种工作： ●数值分析 ●数值和符号计算

技术 ●通讯系统设计与仿真

MATLAB 在通讯系统设计与仿真的应用

●财务与金融工程 MATLAB 的应用范围非常广，包括信号和图像处理、通讯、控制系统设计、测试和测量、财务建模和分析以及计算生物学等众多应用领域。附加的工具箱（单独提供的专用MATLAB 函数集）扩展了MATLAB 环境，以解决这些应用领域内特定类型的问题。

模型的参数估计

模型的参数估计是在建立了理论模型并收集整理了符合模型的样本数据之后，就可以选择适当的方法估计模型，得到模型的参数估计量，本文主要研究线性回归模型，即用回归分析方法建立线性模型。回归分析的主要目的是要通过样本的回归函数尽可能准确的估计总体回归函数。其中使用最广泛的是普通最小二乘法。现以简单的一元线性回归模型为例运用最小二乘法进行模型的参数估计。

1、什么是一元线性回归模型?

一元回归模型, 是最简单的计量经济学模型. 在模型中, 只有一个解释变量, 被解释变量与解释变量之间存在线性关系.

2、一元回归模型的一般形式:

y i =β0+β1x i +μi , i=1,2,…n (1)

其中,y 是被解释变量,x 是解释变量, β0、β1是待定参数, μ是随机误差项,y i ,x i 是随机抽取的n 组样本观测值.

该方程满足如下条件:

E(μi )=0

Var(μi )=б2

Cov(μi , μj )=0

Cov(xi ,u i )=0

i=1,2,...,n j=1,2,…,n i j

3、模型参数估计的任务

(1)一是求得反映变量之间数量关系的参数(即一元线性回归模型y i =β0+β1x i +μi , i=1,2,…n 中的β0, β1) 的估计量;

(2)二是求得随机误差项的分布参数.

4、模型参数估计的普通最小二乘法

普通最小二乘法, 是应用最多的参数估计方法.

(1)什么是最小二乘原理

在已经获得样本观测值y i ,x i (i=1,2,…,n) 的情况下, 假如参数估计量已经求得记为β0, β1, 我们可以得到直线方程:

y i =β0+β1x i i=1,2,…,n (2) 其中, y i 是被解释变量的估计值, 它由参数估计量和解释变量的观测值计算得来.

被解释变量的观测值与估计值, 在总体上越接近越好. 判断的标准是: 二者之差平方和

Q =∑(y i -y i ) 最小.

1n ^2^^^^^^

这就是最小二乘原理.

[思考]为什么用平方和, 而不直接将二者的差简单相加?

(2) 从最小二乘原理, 根据样本观测值, 具体求参数估计值.

由于Q =∑(y i -y i ) , (又 y i =β0+β1x i )

1n ^2^^^

=∑[y i -(β0+β1x i )]

1n ^^2

我们可以知道,Q 是β0, β1二次函数并且是非负数. 所以Q 的极小值总是存在的.(为什么?)

根据极值存在的必要条件知,

^^

⎧⎪∂Q

^=0

⎪⎨∂β0

⎪∂Q

⎪^=0

⎩∂β1

(为什么不是充分条件?)

由此, 不难推得:

⎧

⎨∑(β^^

⎪0+β1x i -y i ) =0

⎪⎩∑(β^+β^

01x i -y i ) x i =0

进而得到:

⎧

⎨y ^^

⎪∑i =n β0+β1∑x i

⎪⎩∑y =β^^ 2

i x i 0∑x i +β1∑x i

于是解得

⎧⎪β^x 2

i y i -x i y i x i

0=

⎪⎨n ∑x 2

i -(∑x i ) 2

⎪^n

⎪β1=y i x i -y i x i

⎩n ∑x 2

i -(∑x i ) 2

另外, 可以将公式(6)简化变形得

⎧∙∙

⎪β^=∑x i y i

⎪1

⎨x ∙2

i ⎪∑ (7)

⎪^__

⎩β^

0=__y -β1x

其中,

x ∙__

i =x i -x ; y ∙=y __

i -y (4) (5) (6)

__x ∑x =n i y ∑; y =__i

n

(3)求随机误差项方差的估计值.

记e i =y i -y i

为第i 个样本观测值的残差. 即被解释变量的观测值与估计值之差. 则随机误差项方差的估计值为: ^

e i ∑=σμn -2（8） 22

证明从略.

至此, 普通最小二乘法一元线性回归模型的参数估计问题得到解决。

M atlab是一种功能强大的系统分析和仿真工

我们选用它作为实现曲线拟合的软件工具

Matlab 语言编程实现最小二乘法的思路:

(1)输入各参量x 、y 的测量值(以数组形式

, 这样便于在计算过程中引用);

(2)用M atlab语言中的plot 函数x 、y 的曲线

图, 以此图对比典型曲线图, 选择合适的经验

;

(3)按照上例中的方法, 选一个系数a, 求

b) 对它的偏导数, 求出其计算表达式;

(4)编写M atlab的M 函数, 用来完成经验公

待定系数a 的计算, 该函数输入量为x 、y 、

量为a 、Q, 按照由最小二乘法推导出的公式

代入数值由x 、y 、b 计算a 、Q;

(5)改变b 的取值, 多次调用该M 函数, 比较

结果中的Q 值, 最小的Q 值所对应的a 、b 值即为

所求。

? 改变b 的取值#这部分工作也可编一个循

环函数, 输入b 可能取的区间, 计算不同b 对应的

Q, 再进行比较, 保留使Q 最小的b 及对应的a 。

但通常b 的改变对Q 的影响不是线性的, 为方便

观察结果并选择适当的b,? 改变b 的取值#这部

分工作最好还是编程者自己完成。

最后, 只要将得到的函数图像和x 、y 的曲线

关系图进行对比, 就可以直观的看到拟合的效

果了。

另外,Matlab 语言提供了一个函数, 可以完成

线性曲线拟合, 这就是函数polyfi 。t 函数polyfit

的输入量为x 、y 、n, 其中x 、y 即为需要建立相互

关系的两个量的测量值, 以数组的形式输入,n 为

多项式的次数; 输出的是多项式系数的行向量, 而

得到的多项式是降幂的。

例1

给出一组数据点(x i , y i ) 列入表1-1中，试用线性最小二乘法求拟合曲线，作出拟合曲线.

表1-1 例1的一组数据(

x , y )

解 （1）在MATLAB 工作窗口输入程序

>> 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）编写下列MA TLAB 程序计算f (x ) 在(x i , y i ) 处的函数值，即输入程序

>> syms a1 a2 a3 a4

x=[-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

运行后屏幕显示关于a 1,a 2, a3和a 4的线性方程组

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

为求a 1, a 2, a 3, a 4使J 达到最小，只需利用极值的必要条件∂J =0 (k =1, 2, 3, 4) ，∂a k

得到关于a 1, a 2, a 3, a 4的线性方程组，这可以由下面的MA TLAB 程序完成，即输入程序

>> syms a1 a2 a3 a4

J=(-125/8*a1+25/4*a2-5/2*a3+a4+1929/10)^2+(-4913/1000*a1+28

9/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 分别对a 1, a2 ,a3 ,a4的偏导数如下

Ja11=

56918107/10000*a1+32097579/25000*a2+1377283/2500*a3+23

667/250*a4-8442429/625

Ja21 =

32097579/25000*a1+1377283/2500*a2+23667/250*a3+67*a4+7

67319/625

Ja31 =

1377283/2500*a1+23667/250*a2+67*a3+18/5*a4-232638/125

Ja41 =

23667/250*a1+67*a2+18/5*a3+18*a4+14859/25

解线性方程组Ja 11 =0，Ja 21 =0，Ja 31 =0，Ja 41 =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.2574

f=[**************]/[**************]*x^3

-[**************]9/[**************]*x^2

+[**************]3/[**************]*x

-[**************]5/[**************]

故所求的拟合曲线为

f (x ) =5.0911x 3-14.1905x 2+6.4102x -8.2574.

（4）编写下面的MA TLAB 程序估计其误差，并作出拟合曲线和数据的图形. 输入程序

>> 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 off

legend(' 数据点(xi,yi)', ' 拟合曲线y=f(x)'),

xlabel('x' ), ylabel('y' ),

title(' 例7.2.1的数据点(xi,yi)和拟合曲线y=f(x)的图形' )

运行后屏幕显示数据(x i , y i ) 与拟合函数f 的最大误差E w ，平均误差E 1和均方根误差E 2及其数据点(x i , y i ) 和拟合曲线y =f (x ) 的图形（略）.

Ew = E1 = E2 =

3.105 4 0.903 4 1.240 9

多重共线性——PLS 偏最小二乘法

经典的最小二乘估计，必需满足一些假设条件，多重共线性就是其中的一种。实际上，解释变量间完全不相关的情形是非常少见的，大多数变量都在某种程度上存在着一定的共线性，而存在着共线性会给模型带来许多不确定性的结果。

设回归模型y =β0+β1x 1+β2x 2+⋯+βp x p

为零的数k 0, k 1, k 2⋯k p 使得k 0

存在完全共线性, 如果k 0

近似的多重共线性。 +ε如果矩阵X 的列向量存在一组不全+k 1x i 1+k 2x i 2+⋯+k p x i p =0, i =1,2,…n , 则称其+k 1x i 1+k 2x i 2+⋯+k p x i p ≈0, i =1,2,…n , 则称其存在

PLS (偏最小二乘法)

H.Wold在1975年提出的 偏最小二乘法近年来引起广泛的关注，在解决多重共线性方面能很好的达到目的，偏最小二乘法集中了最小二乘法、主成分分析法和典型相关分析的的优点克服了两种方法的缺点。偏最小二乘法吸取了主成分回归提取主成分的思想，但不同的是主成分回归只是从自变量中去寻找主成分与因变量无关，因而主成分与因变量在算法上关系不密切，从而导致最后主成分在实际应用中无法更好的进一步拟合因变量，偏最小二乘法则是从因变量出发，选择与因变量相关性较强而又能方便运算的自变量的线性组合。 偏最小二乘回归分析

偏最小二乘回归分析

考虑p 个因变量

L 与m 个自变量y p x m 的建模问题。偏最小二

t t x m 的线性组合，且尽可能多地

u 1并要求u 1与t 1相关程乘回归的基本作法是首先在自变量集中提出第一成分1(1是提取原自变量集中的变异信息）；同时在因变量集中也提取第一成分

度达到最大。然后建立因变量y p 与1的回归，如果回归方程已达到满意的精度，则算法中t

止。否则继续第二对成分的提取，直到能达到满意的精度为止。若最终对自变量集提取r 个成分r ，偏最小二乘回归将通过建立t y p 与r 的回归式，然后再表示为t y p 与原自变量的回

归方程式，即偏最小二乘回归方程式。

偏最小二乘回归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; %求回归分析的系数

beta1(end,:)=[]; %删除回归分析的常数项

cancha=she_f-she_t*beta1; %求残差向量

press_i(j)=sum(cancha.^2);

end

press(i)=sum(press_i)

if i>1

Q_h2(i)=1-press(i)/ss(i-1)

else

Q_h2(1)=1

end

if Q_h2(i)

fprintf('提出的成分个数 r=%d',i);

r=i;

break

end

end

beta_z=[t,ones(num,1)]\f0; %求标准化Y 关于主成分得分向量t 的回归系数 beta_z(end,:)=[]; %删除常数项

xishu=w_star*beta_z; %求标准化Y 关于X 的回归系数， 且是针对标准数据的回归系数，每一列是一个回归方程

mu_x=aver(1:col-1);mu_y=aver(end);

sig_x=stdcov(1:col-1);sig_y=stdcov(end);

ch0=mu_y-mu_x./sig_x*sig_y*xishu; %计算原始数据的回归方程的常数项

xish=xishu'./sig_x*sig_y; %计算原始数据的回归方程的系数，每一列是一个回归方程 Rc=corrcoef(x0*xish'+ch0,y0)

sol=[ch0;xish'] %显示回归方程的系数，每一列是一个方程，每一列的第一个数是常数项

多因变量

function y=pls(pz,Xnum,Ynum)

[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),:); % 标准化自变量 n=Xnum;m=Ynum; %n 是自变量的个数,m 是因变量的个数

x0=pz(:,1:n);y0=pz(:,n+1:end); %提取原始的自变量、因变量数据

e0=data(:,1:n);f0=data(:,n+1:end); %提取标准化后的自变量、因变量数据

num=size(e0,1);%求样本点的个数

temp=eye(n);%对角阵

for i=1:n

%以下计算 w，w*和 t 的得分向量，

matrix=e0'*f0*f0'*e0;

[vec,val]=eig(matrix) %求特征值和特征向量

val=diag(val); %提出对角线元素

[val,ind]=sort(val,'descend');

w(:,i)=vec(:,ind(1)) %提出最大特征值对应的特征向量

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(n)-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; %求回归分析的系数

beta1(end,:)=[]; %删除回归分析的常数项

cancha=she_f-she_t*beta1; %求残差向量

press_i(j)=sum(cancha.^2);

end

press(i)=sum(press_i)

if i>1

Q_h2(i)=1-press(i)/ss(i-1)

else

Q_h2(1)=1

end

if Q_h2(i)

fprintf('提出的成分个数 r=%d',i);

r=i;

break

end

end

beta_z=[t(:,1:r),ones(num,1)]\f0; %求标准化Y 关于 t 的回归系数

beta_z(end,:)=[]; %删除常数项

xishu=w_star(:,1:r)*beta_z; %求标准化Y 关于X 的回归系数， 且是针对标准数据的回归系数，每一列是一个回归方程

mu_x=aver(1:n);mu_y=aver(n+1:end);

sig_x=stdcov(1:n);sig_y=stdcov(n+1:end);

for i=1:m

ch0(i)=mu_y(i)-mu_x./sig_x*sig_y(i)*xishu(:,i); %计算原始数据的回归方程的常数项

end

for i=1:m

xish(:,i)=xishu(:,i)./sig_x'*sig_y(i); %计算原始数据的回归方程的系数， 每一列是一个回归方程

end

sol=[ch0;xish] %显示回归方程的系数，每一列是一个方程，每一列的第一个数是常数项

异方差性——WLS 加权最小二乘法

如果模型被检验证明存在异方差性，则需要发展新的方法估计模型，最常用的方法是加权最小二乘法。加权最小二乘法是对原模型加权，使之变成一个新的不存在异方差性的模型，然后采用普通最小二乘法估计其参数。加权的基本思想是：在采用OLS 方法时，对较小的

22e e i i 残差平方赋予较大的权数，对较大的赋予较小的权数，以对残差提供的信息的重要程

度作一番校正，提高参数的精度。

加权最小二乘法具体步骤是：

⑴ 选择普通最小二乘法估计原模型，得到随机误差项的近似估计量e i ；

⑵ 建立e i 的数据序列；

⑶ 选择加权最小二乘法，以e i 序列作为权，进行估计得到参数估计量。实际上是以e i 乘原模型的两边，得到一个新模型，采用普通最小二乘法估计新模型。