第三章 线性最小二乘估计

1、第三章 线性最小二乘估计 最小二乘估计方法是以误差的平方和最小为准则，根据观测数据估计线性模型中未知参数的一种基本参数估计方法。1794年德国数学家C.F.高斯在解决行星轨道猜测问题时首先提出最小二乘法。它的基本思路是选择估计量使模型（包括静态或动态的，线性或非线性的）输出与实测输出之差的平方和达到最小。这种求误差平方和的方式可以避免正负误差相抵，而且便于数学处理（例如用误差的绝对值就不便于处理）。线性最小二乘法是应用最广泛的参数估计方法，它在理论研究和工程应用中都具有重要的作用，同时它又是许多其他更复杂方法的基础。线性最小二乘法是最小二乘法最简单的一种情况，即模型中关于参数的函数是线性的。

2、假设有一测量方程为，其中为已知的测量数据，为待估计量，为测量噪声，我们可以发现与的关系是线性的。最小二乘方法会告诉我们：如果我们不知道传感器的方差，要在已知测量序列的基础上，怎么样得到估计的值？假设已知其方差为，我们如何应用这个方差，会得到更为准确的估计吗？3.1 最小二乘估计方法 若被估计量是维矢量，则每次观测量和观测噪声均为矢量，线性观测方程为，其中，分别为第次观测量和观测噪声，为观测矩阵。若我们可以得到个观测，将这个观测也可以写成向量形式，记 （3.1）则上式可写为。要求构造的估计量使性能指标 （3.2）达到最小，称这种估计为最小二乘估计，记为。利用微分公式，令，可求出所要求的估计量。将

3、（2）式代入该微分公式，可得则我们可以得到即也就是 （3.3）3.2 线性最小二乘加权估计若观察噪声的均值和协方差分别为，其中，如果不相关，则为对角阵，如果相关，则为方阵，此时更复杂，可取。其原因很简单，就是方差越大的，我们将其的权重设置的越小。线性最小二乘加权估计的性能指标是使达到最小。令，则解上述方程得到 （3.4）3.3 线性最小二乘递推估计 线性最小二乘递推估计的任务是：如果我们已经得到第步的估计，当第步的测量到达后，已知情况下，希望构造估计量，使得的平方和最小，即。 下面我们来推导递推最小二乘的公式。假设次的观察方程为 ，若设 ， ，则线性观测方程为，。我们将k次所有的测量写成如下向

4、量形式 根据最小二乘的基本原理，我们知道，即为求达到最小值时的。其求取过程如下： 设加权矩阵为 则我们可得。根据线性最小二乘加权估计方法（3.4），可知估计矢量为 （3.5）设，则 （3.6）还是根据最小二乘加权估计方法（3.4），我们知道，所要求得的为将，代入上式，我们得到 （3.7）其中即 （3.8）由（3.6）式得 （3.9）由（3.7）式得将（3.9）式代入上式得 (3.10)此时我们需要保留，但消去。将（3.8）式代入（3.10）得即得 现将递推最小二乘估计总结如下：已知，为一合适的数 （3.11）其中为增益，为新息。例1. 回到本节开始的问题，假设有一测量方程为，其中为已知的测量数

5、据，共采集1000个数据，如图3.1所示。为待估计量，其真实值为20。为测量噪声，如果我们不知道传感器的方差，要在已知测量序列的基础上，怎么样得到估计的值？图3.1 传感器采集到的1000个数据解：利用（3.1）式，我们可以得到我们所要研究的问题中，,利用公式（3.3）即可得到估计量，为20.0516。例2. 假设已知其方差为，我们如何应用这个方差，会得到更为准确的估计吗？解：根据题意，得到公式（3.4）中的权值W的取值为1/2。利用公式（3.4），我们得到的估计量仍为20.0516。这是为什么呢？这里为什么加权最小二乘没有效果，没有让我们得到更好的估计结果？我们发现，当权值为标量时，最小二乘

6、加权估计和最小二乘的结果是相同的。其运算过程如下：此时我们用小写w表示公式（3.4）中的权值，加权最小二乘估计方法为例3. 如果是实际数据发生突变，如图3.2所示。其真实值在第501个点处由20跳变至30，那么，利用上述方法得到的结果如何呢？图3.2 传感器得到的具有突变特征的数据解：利用上述最小二乘方法，我们可以得到的估计为：24.9623。我们发现结果与真实值相比，有很大的差距，也就是估计方法并不准确了。其原因在于，我们假设的模型中并没有考虑真值的跳变。其实所得数据与真实值20、30是有很大关系，求得是二者的平均值。下面我们使用递推最小二乘估计方法（3.11）。M的初值可以设置为任意大的正

7、数，待估计值的初值设为0，我们得到的估计结果如图3.3所示。我们看到，与最小二乘、递推最小二乘估计相比，具有递推效果的最小二乘估计方法能够发现真实值的变化，并在估计值中反映这种变化。但经过了后500个的递推估计，还是没有得到与真实值30相近的估计。再进一步试想，如果真实值的变化的复杂一些，如很多信号的变化都是满足的一阶马尔科夫过程，那么其估计的结果会更差。 图3.3 具有突变特征变量的估计结果3.4 最小二乘的其他估计形式 最小二乘的估计方法还有几种等价的估计形式，也就是估计结果是一样的，但计算过程有些不同，在这一节中我们就来看一看。在研究具体的估计形式之前，我们先来看一下递推最小二乘的估计方

8、差。用来表示滤波增益，并在已知测量方差的情况下可知,最小二乘估计为，其估计方差为 （3.12）上式可以推演得到，前一步的估计方差为 将代入（3.12）式，得到再代入测量方程，考虑到待估计值与不相关，经整理后我们可以得到 （3.13）接下来我们来讨论一个很有趣的问题，如果（3.11）的估计方法能够使最小二乘的性能指标最小，那么会使估计方差（3.13）最小吗？反过来，能使估计方差（3.13）的最小滤波增益又是什么呢？还是由使最小二乘的性能指标得到的增益吗？我们期望我们得到的估计当然是一样的，就是使最小二乘性能指标最小的的取值方法也会使估计方差得到最小值，这当然在某种意义上也是最好的估计了。我们现在

9、这样做，通过求解使（3.13）得到最小值的，并判断其与（3.11）中的是否相同来回答上述问题。过程如下：对求导得到 计算上述等式得到 （3.14）这个答案初看起来有些让人惊讶，这和我们已经得到的滤波增益是这样的不同，但先别急，我们继续下面推导： 先设将得到的（3.14）式代入（3.13）式， （3.15）并展开上式得，合并公式后两项，得到即即（3.16）现在注意到（3.14）中的 表达式，隐藏在上式中，所以我们可以把上式重新写一下，为 （3.17）它是的另一个表达式，形式上相对简单，但是数值计算问题可能导致不是正定的，即使是和都是正定的。其实还有第三种表达式，虽然这第三种表达式由于需要三次求逆

10、矩阵运算，在实际估计时并不经常用到，但我们将会发现一个很有趣的现象，下面就来看这种表达式，它是什么呢？利用重写（3.16），得到对上式两边求逆，根据补充知识2中矩阵的求法（b3.1），可知得到（3.18） 我们发现，如果考虑，（3.18）和（3.11）中的第一个公式相同，也就是、，（3.11）中的实际上是估计方差。接下来我们看看（3.14）中的滤波增益与的关系。 首先利用 修改（3.14）式，我们得到将（3.18）式代入，得到注意到在第一个小括号的右边，我们可以将它乘入括号内，得到现在把第一个括号里面的提出来，得到现在在第一个括号左边乘以，再在括号里面乘以，得到 （3.19）可见，考虑到和是相

11、同的，上式中的和利用最小二乘性能指标得到的滤波增益是一致的。到这里，我们会得到如下结论，利用最小二乘性能指标得到的估计方法 （3.20）可以得到最小的估计方差，其中滤波增益和估计方差有不同的计算方法： （3.21） （3.22） （3.23） （3.24）（3.25）计算有三个公式，（3.23）是一个比较简单的表达式，但是数值计算问题可能导致不正定，即使和都是正定的。（3.25）虽然复杂，但是却可以保证是正定的。（3.24）需要三次矩阵求你计算，因此更加复杂，但是在某些情况下却很有用，这在以后的学习中会用到。3.4 小结本章介绍了线性最小二乘估计方法，包括已知测量噪声方差的加权最小二乘方法和递

12、推的最小二乘方法。通过仿真实例我们发现，已知系统的测量模型很重要，当系统的待估计值是恒值的时候，加权最小二乘方法可以得到很好的效果。但当系统变化时，如果测量模型中没有反映这种变化，则估计的效果就不好。递推最小二乘的优点在于，可以反映真实值的变化，但在估计准确度和速度方面还不够好。尤其是在真实值满足马尔科夫随机过程时，没有对真值变化进行建模的最小二乘方法就不在适用了。习题：1. 如果待估计量是一个随着时间变化的量，其中 是待估计量。系统传感器测量噪声方差为2，我们得到的是从0.01秒到10秒之间的1000个测量数据。使用本节学到的估计方法，设计MATLAB程序，给出估计值 。2. 如果待估计量是

13、一个随着时间变化的量，其中 和都是待估计量。系统传感器测量噪声方差为2，我们得到的是从0.01秒到10秒之间的1000个测量数据。使用本节学到的估计方法，设计MATLAB程序，给出估计值和。3. 美国在1946年至1956年的钢产量为66.6，84.9，88.6，78.0，96.8，105.2，93.2，111.6，88.3，117.0，115.2百万吨，分别用一次线性曲线、二次曲线、三次曲线以及四次曲线来拟合这些数据，预测1957至1960年的钢产量。补充知识1：1.矩阵相关公式：设为常矩阵，我们有，则证明：根据乘法公式：2.矩阵逆的求法 （b3.1）3. 矩阵转置计算 补充知识2：下面的M

14、ATLAB程序包含产生上述例题及习题中数据的部分，以便读者可以方便的获取这些数据，并进行仿真。%产生数据的部分clccleart=0.01:0.01:10;t=t' %需要进行转置，方便产生列向量。 A=20;zA=sqrt(2)*randn(1000,1)+A; %产生例1、2所需的数据 za=sqrt(2)*randn(500,1)+20;sqrt(2)*randn(500,1)+30; %产生例3所需的数据 zx=2.*t+sqrt(2)*randn(1000,1); %产生习题1所需的数据 save LSEdata %该语句将所产生的数据保存为.mat数据文件save LSEdata.txtascii %该语句将所产生的数据保存为.txt 数据文件%最小二乘估计部分H=ones(1000,