第四章参数的最小二乘法估计分解

1、第四章 最小二乘法与组合测量1 概述最小二乘法是用于数据处理和误差估计中的一个很得力的数学工具。 对于从 事精密科学实验的人们来说， 应用最小乘法来解决一些实际问题， 仍是目前必不 可少的手段。 例如，取重复测量数据的算术平均值作为测量的结果， 就是依据了 使残差的平方和为最小的原则， 又如，在本章将要用最小二乘法来解决一类组合 测量的问题。 另外，常遇到用实验方法来拟合经验公式， 这是后面一章回归分析 方法的内容，它也是以最小二乘法原理为基础。最小二乘法的发展已经经历了 200 多年的历史，它最先起源于天文和大地测 量的需要， 其后在许多科学领域里获得了广泛应用， 特别是近代矩阵理论与电子

2、计算机相结合，使最小二乘法不断地发展而久盛不衰。本章只介绍经典的最小二乘法及其在组合测量中的一些简单的应用， 一些深 入的内容可参阅专门的书籍和文献。2 最小二乘法原理最小二乘法的产生是为了解决从一组测量值中寻求最可信赖值的问题。 对某 量 x测量一组数据 x1,x2, ,xn ，假设数据中不存在系统误差和粗大误差，相互独立，服从正态分布，它们的标准偏差依次为： 1, 2, n记最可信赖值为 x ，相应的残差 vi xi x 。测值落入 (xi ,xi dx)的概率。Pi2i 2 exp( 2vii2)dx根据概率乘法定理，测量 x1,x2, ,xn 同时出现的概率为1 1 vPPi1 n e

3、xp 1 ( i )2 (dx)nii ( 2 )n2 i i显然，最可信赖值应使出现的概率 P 为最大，即使上式中页指数中的因子达最小，即2vii2 Mini i2权因子： wio2 即权因子 wi 12 ， 则 ii22 wvvwiv i Min再用微分法，得最可信赖值 xnwi xi即加权算术平均值i1xnwii1这里为了与概率符号区别，以 i 表示权因子。特别是等权测量条件下，有：vvvi2Min以上最可信赖值是在残差平方和或加权残差平方和为最小的意义下求得的， 称之为最小二乘法原理。它是以最小二乘方而得名。为从一组测量数据中求得最佳结果，还可使用其它原理 例如(1) 最小绝对残差和法

4、：vi Min(2) 最小最大残差法： max vi Min(3) 最小广义权差法： maxvi min vi Min以上方法随着电子计算机的应用才逐渐引起注意，但最小二乘法便于解析，至今仍用得最广泛 3. 线性参数最小二乘法 先举一个实际遇到的测量问题，为精密测定三个电容值 : x1,x2,x3 采用的测 量方案是，分别等权、独立测得 x1,x2,x1 x3,x2 x3, 列出待解的数学模型。x1=0.3x2 =-0.4x1 + x3 =0.5x2 + x3 =-0.3 这是一个超定方程组，即方程个数多于待求量个数，不存在唯一的确定解， 事实上，考虑到测量有误差，记它们的测量误差分别为 v1

5、,v2,v3,v4 ，按最小二乘 法原理vi2 Min 分别对 x1,x2 ,x3 求偏导数，令它们等于零， 得如下的确定性方程 组。( x1 -0.3)+( x1+x3 -0.5)=0( x2 +0.4)+( x2 +x3 +0.3)=0( x1+x3 -0.5)+( x2 +x3 +0.3)=0可求出唯一解 x1=0.325, x2 =-0.425, x3 =0.150 这组解称之为原超定方程组 的最小二乘解。以下，一般地讨论线性参数测量方程组的最小二乘解及其精度估计。一、正规方程组设线性测量方程组的一般形式为：yiai1x1 ai2 x2ait xt(i 1,2, ,n)y1 a11 x

6、1a12x2a1t xty 2 a 21 x1a22x2a2t xtyn an1x1 an2x2ant xt式中，有 n个直接测得值 y1,y2, ,yn，t 个待求量 x1,x2, , xt 。nt, 各 yi等权，无系统误差和粗大误差固 yi 含有测量误差， 每个测量方程都不严格成立， 故有相应的测量残差方程tvi yiaij xj (i 1,2, , n)j1yi 实测值x j 待估计量，最佳估计值，最可信赖值taij x j 最可信赖的“ y”值。j1按最小二乘法原理，待求的 x j 应满足n n tvvvi2 yiaij xj 2 Mini1 i1 j 1上式分别对 xj 求偏导数，

7、且令其等于零，经推导得a1a1x1a1a2 x2a1at xta1 ya2a1x1a2a2x2a2 at xta2 yata1x1ata2x2atat xtat y式中，aj， y分别为如下列向量aja1 j a2 jyy1y2an jynalak和aj y分别为如下两列向量的内积：al ak =a1l a1k a2la2kanl ankaj y=a1j y1 a2j y2anj yn正规方程组有如下特点：（1）主对角线系数是测量方程组各列系数的平方和，全为正数。（2）其它系数关于主对角线对称3）方程个数等于待求量个数，有唯一解由此可见，线性测量方程组的最小二乘解归结为对线性正规方程组的求解。

8、为了便于进一步讨论问题，下面借助矩阵工具给出正规方程组的矩阵形式。记列向量Yyy1y2Xx1x2v1L=ll122Vv2ynxtvnln和 n t 阶矩阵a11a12a1tAa21a22a2 tan1an 2ant则测量方程组可记为：AX = Y般意义下的方程组测量残差方程组记为当估计出的 xj 已经是最可信赖的值，则 AX 是 yi 的最佳结果最小二乘原理记为(L- AX)T(L- AX) Min利用矩阵的导数及其性质有(V TV ) 2 V V xx2( LxT(XTxAT)(L- AX)xx2 AT ( L - AX )2ATL- 2ATAX令 （V TV ） 0 ，得正规方程组的矩阵形

9、式。 xAT AX = AL展开系数矩阵 AT A和列向量 ATL ，可得代数形式的正规方程组。上述和矩阵的导数有关，因此，我们来分析“矩阵最小二乘法”、矩阵最小二乘法1. 矩阵的导数设n t 阶矩阵a11a12a1tAia21a22a2tani an2ant(aij ) (A1A2 At)n 阶列向量（ n+1 阶矩阵） V 和 t 阶列向量 Xv1x1Vv2Xx2vnxtV 与 X 的转置（行向量）记为 V T 与 X T .关于向量 X 的标量函数。( X) (x1x2 xt)定义如下几个导数（1）矩阵对标量 x 的导数矩阵内A元素aij是x的函数，对矩阵 AX的导数，定义为各元素对 x

10、的导数，构成新的导数矩阵。若 aij 是变量 x 的函数，则定义dAdaijdx ( dx )(E-1)（2）标量函数对向量的导数标量函数 ，对列向量 X 的导数，等于标量函数xi（i 1 t）的导数组成的列向量（行向量的转置）对向量 X 的组成元素yx ( xy xyxy)Txx1 x2xt(E-2)标量函数 ，对行向量 XT 的导数，等于标量函数对向量 X 的组成元素xi（i 1t） 的导数组成的行向量yxT( y yyx1 x2xt) ( yx)Tx(E-3)（3）行（列）向量对列（行）向量的导数V T （v1v2 vn）行向量V T对列向量 X的导数等于行向量各组成元素，对列向量各组成

11、元素分别求得Vv1T(Txv2TxxvTn )Txvx11v1 xt v2 xtV T TXv2x1Xxvnvnxixt(E-5)v1x2vnx1vi v2xxvnxvnx2(E-4)vn关于矩阵的导数有如下性质：1）矩阵 A和B乘积对标量 x 的导数d(AB) AdB dA Bdx dx dx(E-6)(E-7)(E-8)(VTV)x=2VTVx(E-9)(VXTTV) 2VTVXT(E-10)2）常数阵的导数为零矩阵。dA 0dx3）向量关于自身转置向量的导数为单位方阵X T = dXX = dX T4）向量与向量转置乘积的导数5）关于常数矩阵与向量乘积的导数T(X TA) A(E-11)

12、XT (ATX)= AT(E-12)XTT(VTAV)= 2 V AV XX(E-13)T (V T AV ) = 2V T A T(E-14)XX利用（E-1）、（E-4）、和（ E-5）三个定义式，容易证明式（ E-6）、（E-7）、E-8）、和（ E-11）、（E-11）成立。 以下证明式（ E-9）注意到式（ E-2）和式（ E-4）即，标量对列向量求导x(x1 x2T ) xt(E-2)行向量对列向量求导 VXv1 v2vn1 2Xn )(1XXx2x2(E-4)xtxt2v1 vx12vn vxnx1x1式( E-9) 左 ( vi2 )x2v1 vi2vn vnxtxtv1vnx

13、1x1v1vnv12 VT Vxxtxtvn类似地，可以证得式（ E-10）成立。 再证明式（ E-13）注意到V T AV 是关于 x的标量函数，由式（ E-2）知，只需证明(VTAV ) 2 V T AVxixia11 v1 a12 v2a1n vn由于 VT (AxVi ) (v1v2 vn)xixixivann vinxian1 v1an2 v2xixiv1vna11v1a1nv1xian1 v1 vnaxixivnnnvnaxi( vi v2vx1 xiv1xiv1an1vnxivna11v1a1nv1xivn annvnxna11v1a1nv1n)an1vnannvnxiVT AV

14、xi所以式（ E-13）左 V AV +VxiT (AV) 2 VT AV 右xixi2. 正规方程设线性测量方程组与基残差方程组分别为AX = YL - AX = V(E-15) (E-16)式中A为n t阶常数矩阵，X 为t 阶待求向量，L是已知的n阶的测量向量， 注意 l1,l2, ln均是已测量所得），V 是 n阶残差向量。由最小二乘原理VTV = (L- AX)T(L- AX) Min(VTV ) 2 V VXX矩阵性质 (E-9) 式)2( L(X A )(L AX )XX注意到式( E-7)即常数阵的导数为零矩阵。LTX注意到式( E-11)即 (X TA) A，故X(X XA

15、) AT所以TT(VTV) 2AT(L - A)X2(ATL- AT AX)令(V TV ) 0得正规方程组的矩阵形式XATA X = ATL(E-18)当 AT A 满秩的情形，可求出T 1 TX (ATA) 1ATL(E-19)般地，可从式(E-15)出发，用稳定的数值解法，计算 A 的广义逆阵 A 1得(E-20)要进一步去研究此问题，可参阅有关近代矩阵分析及其数值方法的专著3待求量 X 的协方差矩阵已知测量向量 L 协方差矩阵DL E( L EL)( L EL)T =式中， Dl ii为lii 的方差：Dlij 为li 与l j的协的方差：Dl 11Dl12 Dl1nDl 21Dl 2

16、2 Dl 2nDl n1Dl n2 Dl nn22Dlii E(li E(li )2i2Dlijij i j这里，假设 l1,l2, ,ln 为等精度、独立测量的结果，有DL2I利用式( E-19)待求量 X的协方差E(X EX )(X EX)TT 1 T T 1 T T 1 T T 1 T T E( AT A) 1ATL E(ATA)1ATL)(ATA)1ATL E(ATA) 1ATL)TT 1 T T T 1 T T(ATA) 1(ATE(L- EL)(L- EL)T(ATA) 1ATTT 1 T T 1 T(ATA) 1ATDLA(AT A) 1)TT1 T T1(ATA) 1 AT D

17、LA( AT A) 1T1 T 2 T1=(ATA) 1AT 2IA(ATA) 1所以T 1 2DX (AT A) 1 2(E-21)4. 最小二乘法解的最佳性可以证明，在等精度、 独立和无系统误差的测量条件下， 最小二乘法的解具 有唯一性、无偏性、有效性和充分性。证明：1 . 唯一性因测量方程相互独立，且 nt, 则 AT A满秩，式（ E-18）有唯一解2 . 无偏性对 X 的估计式（ E-19 ）求数学期望。E(ATA)-1ATY) (ATA) 1ATEY (ATA) 1ATAx x 3 . 有效性 设另有 X 的无偏估计X* = GY则有EX GE Y GAX X故 G A I又DX*

18、 D(GY) GGTDY GGT I 2 GGT 2而DX (AT A) 1 2引入单位向量0i10其中第 i 行为 1，其它为 0 Xi* 与 Xi 的方差分别为i*2iGGT iT 2i*2i (AT A) iT 2以下证 i* 2i2T T 1 Ti(GGT (AT A) 1) iTi (ATA) 1 (AT A) 1ATGT GA( AT A) 1 GGT ) iTi (AT A) 1AT G)( A( AT A) 1 GT) iTCC T 0其中 第 一等 式利 用 了 GA I ， C i ( AT A) 1AT G) 是一 常数 ，故C CT C2 。最后得证 X 的方差最小，即

19、 X 的有效性成立。4 . 充分性y 取到了测量样本 y1, y2, ,yn 中的所有信息，故按( E-18)式求得 x的估计 量，显然也是充分的。正是由于最小二乘法的解具有最佳性， 所以，最小二乘法在精密测量的各个 领域获得广泛应用。三、精度估计对测量数据的最小二乘法处理，其最终结果不仅要给出待求量的最可信赖 值，还要确定其可信赖程度，即估计其精度。具体内容包含有两方面：一是估计 直接测量结果 y1, y2, , yn 的精度；二是估计待求量 x1,x2, ,xt 的精度。1直接测量结果的精度估计对 t 个未知量的线性测量方程组 AX Y 进行 n 次独立的等精度测量，得 l1,l2, ,l

20、n其残余误差 v1,v2, ,vn标准偏差 。如果 vi服从正态分布，那么vv 2服 从 2分布，其自由度 n-t ，有 2变量的数学期望 E vv/ 2 n t，以 S代 。即有令 t=1 ，由上式又导出了 Bessel 公式2待求量的精度估计按照误差传播的观点，估计量x1,x2, ,xt 的精度取决于直接测量数据l1,l2, ,ln 的精度以及建立它们之间联系的测量方程组。可求待求量的协方差（见二 3）DX(AT A) 1 2矩阵d11d12d1tT1d21d22d2t( AT A )dt1dt2dtt各元素 dij可由矩阵 ATA 求逆得，也可由下列各方程组分别解得。a1a1d 11 a

21、1a2d12a1atd1t 1a2a1d11 a2a2d12a2atd1t 0ata1d11 ata2d12atatd1t 0a1a1d 21 a1a2d22a1atd2t 1 a2a1d21 a2a2d22a2at d2t 0at a1d21 ata2d22atat d2t 0a1a1d t1 a1a2dt2a1at dtt 1(5-51)(5-52)a2a1dt1 a2a2dt 2a2at dtt 0ata1dt1 ata2dt2atat dtt 0是直接测量数据的标准差，可按S vv 估计nt待求量 xj 的方差22x2j d jj 2 (j 1,2, ,t)矩阵(AT A) 1中对角元

22、素 d jj就是误差传播系数。待求量 xi 与 xj 的相关系数。dijij (i, j 1,2, ,t)diid jj现在，可以解决本节开始提出的测量问题例 5-1 为精密测定 1 号、2 号和 3 号电容器的电容 x1, x2, x3 ，进行了等权独 立、无系统误差的测量。测得 1号电容值 C1 =0.3 ，2 号电容值 C2=-0.4 ，1 号和 3 号并联电容值 C3 =0.5 ， 2 号和 3 号并联电容值 C4 =-0.3 。试用最小二乘法求x1, x2, x3及其标准差。解：列出残差方程组v1 0.3 x1v2 0.4 x2v3 0.5 (x1 x2)v4 0.3 (x2 x3)

23、为计算方便，将数据列表如下：iai1ai2ai3yiai1ai1ai1ai 2ai1ai3ai1yiai2ai2ai2ai 3ai2yiai3ai3ai3yi11000.31000.3000002010-0.4000010-0.40031010.51010.500010.54011-0.3000011-0.31-0.32010.821-0.720.2按上表计算正规方程组各系数和常数项后，列出正规方程组2x1 0x2 x3 0.80x1 2x2 x3 0.7x1 x2 2x3 0.2解出 x1 =0.325, x2 =-0.425, x3 =0.150 代入残差方程组，计算v1 v2v3 v40

24、.025vv v12 v22 v32 v42 0.00250.00250.05 按式( 5-51 )，求出d11 =0.75, d22 =0.75, d33 =1 按式( 5-52 )，求出 xj d jjx1 0.0433， x2 0.0433 ， x3 0.050 最后得 1 号、2号和 3 号电容器的精密电容值x1 0.325 3 x1, x20.425 3 x2 ， x30.150 3 x3也可以用矩阵形式，这里显然：10A0100110.30.4Y0.50.3这样可求得 AT ,AT A,( AT A) 1求逆阵：则 X (ATA) 1ATY由 Y AX V 求得V vi (i 14

25、)由 D (ATA) 1，djj 可得d11,d22,d33 xjd jj 写出结果 4 非线性参数的最小二乘法在例 5-1 中，除了进行 4 次测量外，又对 1 号和 2 号电容器的串联电容 x1x2 /(x1 x2) 进行测量，测得 y5 ，方差仍为 2 ，那么如何处理呢？简单的办法 是把它线性化。所谓线性化，就是在未知量的附近，按泰勒级数展开取一次项，然后按线性参数最小二乘法进行迭代求解。线性化的具体步骤如下：设测量残差方程组yii (x1,x2, ,xt ) vi(4-1)(4-2)取xj 的初始近似值 Cj记j xj Cj则有yii (x1,x2, ,xt ) vii (C 1,C2

26、, ,Ct)1 x1 C xtt viCyi yii (C1,C2, ,Ct )(4-3)ai1C ，ai2xiC ， aitxt(4-4)于是得线性化残差方程组yiai1 1ai2 2ait t vi(4-5)作法：按线性参数最小二乘法解得 j，以至 xj，将此 xj作为新的 Cj，按式4-2 ），式（ 4-3 ），式（ 4-4 ）和式（ 4-5 ）进行反复迭代求解，直至 j 符合精 度要求为止。例 5-2 在例 5-1 的基础上，再增加一次测量串联电容 x1x2/（x1 x2） ，测得 y5 =0.14 。试用最小二乘法求 x1, x2 , x3及其标准差解：先列出测量方程组x1 =0.3

27、x1+x3 =0.5x2 =-0.4x2 +x3 =-0.3x1x2(x1 x2)0.14对前4个线性测量方程组，按例 5-1 求出解，作为初次近似解C1(1) 0.325, C2(1)0.425, C3(1) 0.150在(0.325,-0.425,0.150) 附近，取泰勒展开的一阶近似，y1 0.3 0.325 0.025, a111, a12 a13 0x1y2 0.4 0.425 0. 025, a21 a23 0, a 22 1y 3 0.5 (0.325 0.150) 0.025, a31 a33 1, a 32 0y4 0.3 ( 0.425 0.150) 0.025, a 4

28、1 0, a42 a43 1 0.325 0.425y5 0.14 ( ) 1.241255 0.325 0.4252a51x2 2 18.062551 (x1 x2 )22a5210.5625,x1(x1 x2 )2写出线性化残差方程组ALV1000.025v101010.025v210120.025v301130.025v418.062510.562501.24125v5整理得正规方程组328.254 190.785 1122.4201190.785 113.566 1213.11071 1 230解出1 0.0473, 2 0.0363, 3 0.0418取xj 的二次近似值C1(2)

29、C1(1) 1 0.2777, C2(2) C2(1) 2 0.4613, C3(2) 0.1918 重复上述过程再求出 1, 2 和 3 。依次迭代结果如表所示迭代次数123x1x2x300000.325-0.4250.1501-0.0473-0.03630.04180.278-0.4610.1922-0.0713-0.03730.05430.206-0.4990.2463-0.0472-0.05550.02640.159-0.5040.27340.001980.00105-0.006280.161-0.4940.2665-0.00113-0.001420.001270.160-0.4950

30、.26860.0003150.000419-0.0003670.160-0.4950.267可见，经 6 次迭代，精度已达 10-3，满足要求即可结束迭代 5 组合测量问题所谓组合测量， 是指直接或间接测量一组被测量的不同组合值， 从它们相互 组合所依赖的若干函数关系中， 确定出各被测量的最佳估计值。 组合测量的问题 常用最小二乘法， 以上两节所举精密测量电容值的问题就是一例。 本节再介绍几 个实例，以进一步说明组合测量方法的特点。例 4-3 如图所示，要求检定线纹尺 0,1,2,3 刻线间的距离 x1,x2,x3。已知用组 全测量法测得图所示刻线间隙的各种组合量。L1=1.01, L2=0.98, L3=1.02L4=2.02,L5=1.98,L6=3.03解：按前述方法，可以解得x1=1.028(0.011)，x2=0.983(0.011)，x3=1.013(0.011)这里，着重说明组合测量方法的优点。本例对刻度间隔 x1,x2与 x3分别测了 3 次，总共测量 6次若不采用组合测量， 按每刻度间隔重复测量 3次计，共需作 9 次测量，比组 合测量法多测 3 次