参数的最小二乘法估计讲解

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

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

3、 1 1 v PPi1 n exp 1 ( i )2 (dx)n ii ( 2 )n2 i i 显然，最可信赖值应使出现的概率 P 为最大，即使上式中页指数中的因子达 最小，即 2 vi i2 Min i i2 权因子： wi o2 即权因子 wi 12 ， 则 ii2 2 wvvwiv i Min 再用微分法，得最可信赖值 x n wi xi 即加权算术平均值 i1 xn wi i1 这里为了与概率符号区别，以 i 表示权因子。 特别是等权测量条件下，有： vvvi2Min 以上最可信赖值是在残差平方和或加权残差平方和为最小的意义下求得的， 称之为最小二乘法原理。它是以最小二乘方而得名。 为

4、从一组测量数据中求得最佳结果，还可使用其它原理 例如 (1) 最小绝对残差和法：vi Min (2) 最小最大残差法： max vi Min (3) 最小广义权差法： maxvi min vi Min 以上方法随着电子计算机的应用才逐渐引起注意，但最小二乘法便于解析， 至今仍用得最广泛 3. 线性参数最小二乘法 先举一个实际遇到的测量问题，为精密测定三个电容值 : x1,x2,x3 采用的测 量方案是，分别等权、独立测得 x1,x2,x1 x3,x2 x3, 列出待解的数学模型。 x1=0.3 x2 =-0.4 x1 + x3 =0.5 x2 + x3 =-0.3 这是一个超定方程组，即方程个

5、数多于待求量个数，不存在唯一的确定解， 事实上，考虑到测量有误差，记它们的测量误差分别为 v1,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 这组解称之为原超定方程组 的最小二乘解。 以下，一般地讨论线性参数测量方程组的最小二乘解及其精度估计。 一、正规方程组

6、 设线性测量方程组的一般形式为： yiai1x1 ai2 x2 ait xt (i 1,2, ,n) y1 a11 x1 a12x2a1t xt y 2 a 21 x1 a22x2a2t xt yn an1x1 an2x2ant xt 式中，有 n个直接测得值 y1,y2, ,yn，t 个待求量 x1,x2, ,xt 。nt, 各 yi等 权，无系统误差和粗大误差 固 yi 含有测量误差， 每个测量方程都不严格成立， 故有相应的测量残差方程 t vi yiaij xj (i 1,2, , n) j1 yi 实测值 x j 待估计量，最佳估计值，最可信赖值 t aij x j 最可信赖的“ y”

7、值。 j1 按最小二乘法原理，待求的 x j 应满足 n n t vvvi2 yiaij xj 2 Min i 1 i 1 j 1 上式分别对 xj 求偏导数，且令其等于零，经推导得 a1a1x1a1a2 x2a1at xta1 y a2a1x1a2a2x2a2 at xta2 y ata1x1ata2x2atat xtat y 式中， aj， y分别为如下列向量 aj a1 j a2 j y y1 y2 an j yn alak和aj y分别为如下两列向量的内积： al ak =a1l a1k a2la2kanl ank aj y=a1jy1 a2jy2anjyn 正规方程组有如下特点： （

8、1）主对角线系数是测量方程组各列系数的平方和，全为正数。 （2）其它系数关于主对角线对称 3）方程个数等于待求量个数，有唯一解 由此可见，线性测量方程组的最小二乘解归结为对线性正规方程组的求解。 为了便于进一步讨论问题，下面借助矩阵工具给出正规方程组的矩阵形式。 记列向量 Y yy1 y2 X x1 x2 v1 L= ll12 2 V v2 yn xt vn ln 和 n t 阶矩阵 a11a12a1t A a21a22a2 t an1an 2ant 则测量方程组可记为： AX = Y 般意义下的方程组 测量残差方程组记为 当估计出的 xj 已经是最可信赖的值，则 AX 是 yi 的最佳结果

9、最小二乘原理记为 (L- AX)T(L- AX) Min 利用矩阵的导数及其性质有 (V TV ) 2 V V xx 2( LxT(XTxAT)(L- AX) xx 2 AT ( L - AX ) 2ATL- 2ATAX 令 （V TV ） 0 ，得正规方程组的矩阵形式。 x AT AX = AL 展开系数矩阵 AT A和列向量 ATL ，可得代数形式的正规方程组。 上述和矩阵的导数有关，因此，我们来分析“矩阵最小二乘法” 、矩阵最小二乘法 1. 矩阵的导数 设n t 阶矩阵 a11a12a1t Ai a21a22a2t ani an2ant (aij ) (A1A2 At) n 阶列向量（

10、n+1 阶矩阵） V 和 t 阶列向量 X v1 x1 V v2 X x2 vn xt 定义如下几个导数 （1）矩阵对标量 x 的导数 矩阵内A元素aij是x的函数，对矩阵 AX的导数，定义为各元素对 x的导数， 构成新的导数矩阵。 若 aij 是变量 x 的函数，则定义 dAdaij dx ( dx ) (E-1) （2）标量函数对向量的导数 标量函数 ，对列向量 X 的导数，等于标量函数 xi（i 1 t）的导数组成的列向量（行向量的转置） 对向量 X 的组成元素 yx ( xy xyxy)T xx1 x2xt (E-2) 标量函数 ，对行向量 XT 的导数，等于标量函数 对向量 X 的组

11、成元素 xi（i 1t） 的导数组成的行向量 y xT ( y yy x1 x2xt ) ( yx)T x (E-3) （3）行（列）向量对列（行）向量的导数 V T （v1v2 vn） 行向量V T对列向量 X的导数等于行向量各组成元素，对列向量各组成元素 分别求得 Vv1 T(T xv2T x xvTn )T x vx11 v1 xt v2 xt V T T X v2 x1 X x vn vn xi xt (E-5) v1 x2 vn x1 vi v2 xx vn x vn x2 (E-4) vn 关于矩阵的导数有如下性质： 1）矩阵 A和B乘积对标量 x 的导数 d(AB) AdB dA

12、 B dx dx dx (E-6) (E-7) (E-8) (VTV) x =2 VTV x (E-9) (VXTTV) 2VT V XT (E-10) 2）常数阵的导数为零矩阵。 dA 0 dx 3）向量关于自身转置向量的导数为单位方阵 X T = dX X = dX T 4）向量与向量转置乘积的导数 5）关于常数矩阵与向量乘积的导数 T (X TA) A (E-11) X T (ATX)= AT (E-12) XT T (VTAV)= 2 V AV XX (E-13) T (V T AV ) = 2V T A T (E-14) XX 利用（E-1）、（E-4）、和（ E-5）三个定义式，容

13、易证明式（ E-6）、（E-7）、 E-8）、和（ E-11）、（E-11）成立。 以下证明式（ E-9） 注意到式（ E-2）和式（ E-4）即， 标量对列向量求导 x( x1 x2 T ) xt (E-2) 行向量对列向量求导 V X v1 v2vn 1 2Xn ) (1 XX x2 x2 (E-4) xt 式( E-9) 左 ( vi2 ) x 2v1 vi2vn vn xt xt v1vn x1x1 v1vn v1 2 VT V x xtxt vn 类似地，可以证得式（ E-10）成立。 再证明式（ E-13） 注意到V T AV是关于x的标量函数，由式（ E-2）知，只需证明 (VT

14、AV ) 2 V T AV xixi a11 v1 a12 v2a1n vn 由于 VT (AxVi ) (v1v2 vn) xixixi v ann vin xi an1 v1an2 v2 xixi v1vn a11v1a1nv1 xi an1 v1 vna xi xi vn nnvn axi ( vi v2v x1 xi v1 xi v1 an1vn xi vn a11v1a1nv1 xi vn annvn xn a11v1a1nv1 n) an1vnannvn xi VT AV xi 所以式（ E-13）左 V AV +V xi T (AV) 2 VT AV 右 xi xi 2. 正规方

15、程 设线性测量方程组与基残差方程组分别为 AX = Y L - 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 V XX 矩阵性质 (E-9) 式) 2( L(X A )(L AX ) XX 注意到式( E-7)即常数阵的导数为零矩阵。 LT X 注意到式( E-11)即 (X TA) A，故 X (X XA ) AT 所以 TT (VTV) 2AT(L - A) X 2(

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

17、ii E(li E(li )2i2 Dlij 为li 与l j的协的方差： Dlij ij i j 这里，假设 l1,l2, ,ln 为等精度、独立测量的结果，有 DL2I 利用式( E-19)待求量 X的协方差 E(X EX )(X EX)T T 1 T T 1 T T 1 T T 1 T T E( AT A) 1ATL E(ATA)1ATL)(ATA)1ATL E(ATA) 1ATL)T T 1 T T T 1 T T (ATA) 1(ATE(L- EL)(L- EL)T(ATA) 1ATT T 1 T T 1 T (ATA) 1ATDLA(AT A) 1)T T1 T T1 (ATA)

18、1 AT DLA( AT A) 1 T1 T 2 T1 =(ATA) 1AT 2IA(ATA) 1 所以 T 1 2 DX (AT A) 1 2 (E-21) 对 X 的估计式（ E-19 ）求数学期望。 E(ATA)-1ATY) (AT A) (ATA) 1 (AT A) 1ATGT GA( AT A) 1 GGT ) iT AT EY (ATA) 1ATAx x 3 . 有效性 设另有 X 的无偏估计 X*= GY 则有 EX GE Y GAX X 故 G A I 又 DX* D(GY) GGTDY GGT I 2 GGT 2 而DX (AT A) 1 2 引入单位向量 0 i1 0 其中

19、第 i 行为 1，其它为 0 Xi* 与 Xi 的方差分别为 i*2iGGT iT 2 i*2i (AT A) iT 2 以下证 i* 2i2 T T 1 T i(GGT (AT A) 1) iT i (AT A) 1AT G)( A( AT A) 1 GT) iT CC T 0 其中第一等式利用了 GA I ， C i ( AT A) 1AT G) 是一常数，故 C CT C2 。最后得证 X 的方差最小，即 X 的有效性成立。 4 . 充分性 y 取到了测量样本 y1, y2, ,yn 中的所有信息，故按( E-18)式求得 x的估计 量，显然也是充分的。 正是由于最小二乘法的解具有最佳性

20、， 所以，最小二乘法在精密测量的各个 领域获得广泛应用。 三、精度估计 对测量数据的最小二乘法处理，其最终结果不仅要给出待求量的最可信赖 值，还要确定其可信赖程度，即估计其精度。具体内容包含有两方面：一是估计 直接测量结果 y1, y2, , yn 的精度；二是估计待求量 x1,x2, ,xt 的精度。 1直接测量结果的精度估计 对 t 个未知量的线性测量方程组 AX Y 进行 n 次独立的等精度测量，得 l1,l2, ,ln其残余误差 v1,v2, ,vn标准偏差 。如果 vi服从正态分布，那么vv 2服 从 2分布，其自由度 n-t ，有 2变量的数学期望 E vv/ 2 n t，以 S代

21、 。 即有 令 t=1 ，由上式又导出了 Bessel 公式 2待求量的精度估计 按照误差传播的观点，估计量 x1,x2, ,xt 的精度取决于直接测量数据 l1,l2, ,ln 的精度以及建立它们之间联系的测量方程组。 可求待求量的协方差（见二 3） DX (AT A) 1 2 矩阵 d11 d12 d1t T1 d21 d22 d2t ( AT A ) dt1 dt2 dtt 各元素 dij可由矩阵 ATA 求逆得，也可由下列各方程组分别解得。 a1a1d 11 a1a2d12a1atd1t 1 a2a1d11 a2a2d12a2atd1t 0 ata1d11 ata2d12atatd1t

22、 0 a1a1d 21 a1a2d22a1atd2t 1 a2a1d21 a2a2d22a2at d2t 0 at a1d21 ata2d22atat d2t 0 a1a1d t1 a1a2dt2a1at dtt 1 (5-51) (5-52) a2a1dt1 a2a2dt 2a2at dtt 0 ata1dt1 ata2dt2atat dtt 0 是直接测量数据的标准差，可按 S vv 估计 nt 待求量 xj 的方差 22 x2j d jj 2 (j 1,2, ,t) 矩阵(AT A) 1中对角元素 d jj就是误差传播系数。 待求量 xi 与 xj 的相关系数。 dij ij (i, j

23、 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 x1 v2 0.4 x2 v3 0.5 (x1 x2) v4 0.3 (x2 x3) 为计算方便，将数据列表如下： i ai1 ai2 ai3 yi a

24、i1ai1 ai1ai 2 ai1ai3 ai1yi ai2ai2 ai2ai 3 ai2yi ai3ai3 ai3yi 1 1 0 0 0.3 1 0 0 0.3 0 0 0 0 0 2 0 1 0 -0.4 0 0 0 0 1 0 -0.4 0 0 3 1 0 1 0.5 1 0 1 0.5 0 0 0 1 0.5 4 0 1 1 -0.3 0 0 0 0 1 1 -0.3 1 -0.3 2 0 1 0.8 2 1 -0.7 2 0.2 按上表计算正规方程组各系数和常数项后，列出正规方程组 2x1 0 x2 x3 0.8 0 x1 2x2 x3 0.7 x1 x2 2x3 0.2 解出 x

25、1 =0.325, x2 =-0.425, x3 =0.150 代入残差方程组，计算 v1 v2v3 v40.025 vv v12 v22 v32 v42 0.0025 0.0025 0.05 按式( 5-51 )，求出 d11 =0.75, d22 =0.75, d33 =1 按式( 5-52 )，求出 xj d jj x1 0.0433， x2 0.0433 ， x3 0.050 最后得 1 号、2号和 3 号电容器的精密电容值 x1 0.325 3 x1, x20.425 3 x2 ， x30.150 3 x3 也可以用矩阵形式，这里显然： A 10 01 0 0 1 1 0.3 0.4

26、 Y 0.5 0.3 这样可求得 AT ,AT A,( AT A) 1求逆阵： 则X (ATA) 1ATY 由 Y AX V 求得V vi (i 14) 由D (ATA) 1，djj可得d11,d22,d33 xj d jj 写出结果 4 非线性参数的最小二乘法 在例 5-1 中，除了进行 4 次测量外，又对 1 号和 2 号电容器的串联电容 x1x2 /(x1 x2) 进行测量，测得 y5 ，方差仍为 2 ，那么如何处理呢？简单的办法 是把它线性化。所谓线性化，就是在未知量的附近，按泰勒级数展开取一次项，然后按线性参数最小二乘法进行迭代求解。 线性化的具体步骤如下： 设测量残差方程组 yii

27、 (x1,x2, ,xt ) vi (4-1) (4-2) 则有 yii (x1,x2, ,xt ) vi i (C 1,C2, ,Ct) 1 x1 C xt t vi C yi yii (C1,C2, ,Ct ) (4-3) ai1 C ，ai2xi C ， aitxt (4-4) 取xj 的初始近似值 Cj记 j xj Cj 于是得线性化残差方程组 yiai1 1ai2 2ait t vi (4-5) 作法：按线性参数最小二乘法解得j ，以至 xj，将此 xj 作为新的 Cj，按式 4-2 ），式（ 4-3 ），式（ 4-4 ）和式（ 4-5 ）进行反复迭代求解，直至 j 符合精 度要求为

28、止。 例 5-2 在例 5-1 的基础上，再增加一次测量串联电容 x1x2 /（x1 x2） ，测得 y5 =0.14 。试用最小二乘法求 x1, x2 , x3及其标准差 解：先列出测量方程组 x1 =0.3 x1+x3 =0.5 x2 =-0.4 x2 +x3 =-0.3 x1x2 (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 0

29、x1 y2 0.4 0.425 0. 025, a21 a23 0, a 22 1 y 3 0.5 (0.325 0.150) 0.025, a31 a33 1, a 32 0 y4 0.3 ( 0.425 0.150) 0.025, a 41 0, a42 a43 1 0.325 0.425 y5 0.14 ( ) 1.24125 5 0.325 0.425 2 a51x2 2 18.0625 51 (x1 x2 )2 a52 2 x1 (x1 x2 )2 10.5625, 写出线性化残差方程组 A L V 1 0 0 0.025 v1 0 1 0 1 0.025 v2 1 0 1 2 0.

30、025 v3 0 1 1 3 0.025 v4 18.0625 10.5625 0 1.24125 v5 整理得正规方程组 328.254 190.785 1 1 22.4201 190.785 113.566 1 2 13.1107 1 1 2 3 0 解出 1 0.0473, 2 0.0363, 3 0.0418 取xj 的二次近似值 C1(2) C1(1) 1 0.2777, C2(2) C2(1) 2 0.4613, C3(2) 0.1918 重复上述过程再求出 1, 2 和 3 。 依次迭代结果如表所示 迭代次数 1 2 3 x1 x2 x3 0 0 0 0 0.325 -0.425

31、 0.150 1 -0.0473 -0.0363 0.0418 0.278 -0.461 0.192 2 -0.0713 -0.0373 0.0543 0.206 -0.499 0.246 3 -0.0472 -0.0555 0.0264 0.159 -0.504 0.273 4 0.00198 0.00105 -0.00628 0.161 -0.494 0.266 5 -0.00113 -0.00142 0.00127 0.160 -0.495 0.268 6 0.000315 0.000419 -0.000367 0.160 -0.495 0.267 可见，经 6 次迭代，精度已达 10-

32、3，满足要求即可结束迭代 5 组合测量问题 所谓组合测量， 是指直接或间接测量一组被测量的不同组合值， 从它们相互 组合所依赖的若干函数关系中， 确定出各被测量的最佳估计值。 组合测量的问题 常用最小二乘法， 以上两节所举精密测量电容值的问题就是一例。 本节再介绍几 个实例，以进一步说明组合测量方法的特点。 例 4-3 如图所示，要求检定线纹尺 0,1,2,3 刻线间的距离 x1,x2,x3。已知用组 全测量法测得图所示刻线间隙的各种组合量。 L1=1.01, L2=0.98, L3=1.02 L4=2.02,L5=1.98,L6=3.03 解：按前述方法，可以解得 x1=1.028(0.01

33、1)，x2=0.983(0.011)，x3=1.013(0.011) 这里，着重说明组合测量方法的优点。 本例对刻度间隔 x1,x2与 x3分别测了 3 次，总共测量 6次 若不采用组合测量， 按每刻度间隔重复测量 3次计，共需作 9 次测量，比组 合测量法多测 3 次。如果待检定的刻度间隔远多于 3 个。那么可以类似分析得出， 采用组合测量法可以大大减少测量次数，提高测量的工作效率。 本 例 测 量 方 程 的 个 数 是 6 ， 待 求 量 的 个 数 是 3 。 假 设 v2 v1 v2 v3 v4 v5 v6 v 。按 S 有 S 2v 。如果测量方程减少为 4 nt 个，那么有S2 4v2。如果两种情形的误差传播系数 djj相近，那么按式 xjdjj