最小二乘法综述及算例一最小二乘法的历史简介1801年,意大利天文学家朱赛普·皮亚齐发现了第一颗小行星谷神星。
经过40天的跟踪观测后,由于谷神星运行至太阳背后,使得皮亚齐失去了谷神星的位置。
随后全世界的科学家利用皮亚齐的观测数据开始寻找谷神星,但是根据大多数人计算的结果来寻找谷神星都没有结果。
时年24岁的高斯也计算了谷神星的轨道。
奥地利天文学家海因里希·奥尔伯斯根据高斯计算出来的轨道重新发现了谷神星。
高斯使用的最小二乘法的方法发表于1809年他的著作《天体运动论》中。
经过两百余年后,最小二乘法已广泛应用与科学实验和工程技术中,随着现代电子计算机的普及与发展,这个方法更加显示出其强大的生命力。
二最小二乘法原理最小二乘法的基本原理是:成对等精度测得的一组数据),...,2,1(,n i y x i i =,是找出一条最佳的拟合曲线,似的这条曲线上的个点的值与测量值的差的平方和在所有拟合曲线中最小。
设物理量y 与1个变量l x x x ,...,2,1间的依赖关系式为:)(,...,1,0;,...,2,1n l a a a x x x f y =。
其中n a a a ,...,1,0是n +l 个待定参数,记()21∑=-=mi i i y vs 其中 是测量值, 是由己求得的n a a a ,...,1,0以及实验点),...,2,1)(,...,(;,2,1m i v x x x i il i i =得出的函数值)(,...,1,0;,...,2,1n il i i a a a x x x f y =。
在设计实验时, 为了减小误差, 常进行多点测量, 使方程式个数大于待定参数的个数, 此时构成的方程组称为矛盾方程组。
通过最小二乘法转化后的方程组称为正规方程组(此时方程式的个数与待定参数的个数相等) 。
我们可以通过正规方程组求出a最小二乘法又称曲线拟合, 所谓“ 拟合” 即不要求所作的曲线完全通过所有的数据点, 只要求所得的曲线能反映数据的基本趋势。
三曲线拟合曲线拟合的几何解释: 求一条曲线, 使数据点均在离此曲线的上方或下方不远处。
(1)一元线性拟合设变量y 与x 成线性关系x a a y 10+=,先已知m 个实验点),...,2,1(,m i v x i i =,求两个未知参数1,0a a 。
令()2110∑=--=mi i i x a a y s ,则1,0a a 应满足1,0,0==∂∂i a si。
即i v i v化简得从中解出∑∑∑∑====-==⎪⎪⎪⎭⎫⎝⎛=-=-=∑∑∑m i m i mi mi ii mi m i i i mi i i i i x m a y m a xxm y x y x m a 112111011111(2)多元线性拟合设变量y 与n 个变量)1(,...,2,1≥n n xx x 的内在联系是线性的,即有下式∑==+nj j o x a a y 11设j x 的第i 次测量值为ij x ,对应的函数值为),...,2,1(m i i y ==,则偏差平方和()()∑∑=--=-==mi i mi i i x a a y y y s 11021'为使s 去得最小值的方程组⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧=⎪⎭⎫ ⎝⎛∑-=∂∂=⎪⎭⎫ ⎝⎛∑-=∂∂=⎪⎭⎫ ⎝⎛∑-=∂∂∑∑∑=--=--=--===m i in ij j i n m i i ij j i m i ij i x x a a y a s x x a a y a s x a a y a s n j n j n j 101101110002....................................................0202111 即⎪⎪⎩⎪⎪⎨⎧=⎪⎭⎫ ⎝⎛+=⎪⎭⎫⎝⎛+∑∑∑∑∑∑∑=======n j mi ik j m i ik ij mi ik mi ij n j m i ij y x a x x a x y a x ma 111101110n k ,...,2,1=。
(4) 将实验数据()i ij y x ,代入(4)式,即得n a a a ,...,1,0。
∑∑=--=--=-=∂∂=-=∂∂m i i i mi i i i x a a y a sx a a y a s11011000)(20)(2∑∑∑∑====+=+mi m i m i m i ii i ii y x a x a y m x m a a 111110101(3)多项式拟合科学实验后得到一组数据时,常会遇到因变量y 与自变量x 之间根本不存在线性关系。
此可以考虑用一个n 次多项式来拟合y 与x 之间的函数关系。
对于n 次多项式∑==ni ii xa y 0,令),...,1,0(n x x j ij ==,则可将其化为线性形式:∑=+=nj j j x a a y 1对于i=1,2,...,m 个实验点有j i ij x x =,代入(3)式有n k y x a x x a x y a x ma mi n j mi i ik j m i ik ij ik mi ij n j m i ij ,...,2,1111101110=⎪⎪⎩⎪⎪⎨⎧=⎪⎭⎫ ⎝⎛+=⎪⎭⎫⎝⎛+∑∑∑∑∑∑∑======= 从而得出多项式的最小二乘法拟合的方程n k y x a x i mi k i i ni m i k j i ,...,1,0111==⎪⎭⎫⎝⎛∑∑∑===+写成矩阵的形式即为⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫⎝⎛∑∑∑∑∑∑∑∑∑∑∑∑∑∑====+=+==+======mi i n i m i i i m i i n m i n i mi n imi n im i n i m i n i mi imi im i i mi mi nimi iiy x y x y a a a x xxx x xxx x xxm 1110121211111131211112.............................. 从中可以解出n a a a ,...,1,0。
(4)指数函数拟合此时拟合函数具有形式bxae y =(a ,b 为待定系数)。
两端取自然对数有(*)ln ln bxa y +=令a b yY ln ln 0==则(*)式化为线性形式 bx b Y +=0再利用(1)式和(2)式,即可求出b b ,0。
从而有ob e a =。
故bxb o e y +=。
四最小二乘法应用举例例:已知某铜棒的电阻与温度关系为:t R R t ⋅+=α0。
实验测得7组数据(见表1)如下:试用最小二乘法求出参量R 0、α 以及确定它们的误差。
此例中只有两个待定的参量R 0和α,为得到它们的最佳系数,所需要的数据有n 、∑i x 、∑iy 、∑2ix、∑2iy和∑iiyx 六个累加数,为此在没有常用的科学型计算器时,通过列表计算的方式来进行,这对提高计算速度将会有极大的帮助(参见表2),并使工作有条理与不易出错。
其中表内双线右边的计算是为了确定R 0和α的误差项用的。
根据表2中所求得的数据,代入公式(12))则可得: C k 02/28788.035.51156.1472)5.245(8.9340700.5665.2458.200607Ω==-⨯⨯-⨯==αΩ=⋅-==76078.7075.24528788.0700.5660b R 把测量数据代入式(13)和(15)中可求出相关系数)]7)00.566(45826[(]7)5.245(8.9340[700.5665.2458.20060])(1[])(1[1222222-⨯-⨯-=-⋅--=∑∑∑∑∑∑∑i i i i ii i i y n y x n x y x n y x γ99757.07)00.566(458267)5.245(8.934028788.0)(1)(1222222=--⨯=--⨯=∑∑∑∑i i i iy n y x n x k说明:电阻R t 与温度t 的线性关系良好,所以取R 0的有效数字与R 对齐,即R 0=70.76Ω;又因为t 7-t 1 = 31.0 ℃,R 7-R 1 = 8.80Ω,取k 有效数字为以上两个差值中较少的位数3位,则k = 0.288Ω/︒C 。
由此可以得到电阻与温度的相关关系为:t R t 288.076.70+=按补充资料中的公式计算k 和b 的不确定度,可得)(239.027102845242Ω=-⨯=-==-∑n S S iR y t ν)C /(0088.003699.0239.07)5.245(8.9340239.0)(222︒Ω=⨯=-=-==∑∑nx x S S S i i yk α)(33.078.93400088.020Ω=⨯===∑n xS S S ikR b故 Ω±=Ω±=)3.08.70()33.076.70(0R ,C /)009.0288.0(C /)009.02879.0(︒Ω±=︒Ω±=α 则 t R t 288.08.70+=参考文献:1.《最小二乘法与测量平差》 郭禄光,樊功瑜著 同济大学出版社 19852.《近代最小二乘法》 测绘出版社 19803.《最小二乘法的拟合及其应用》 邓亮章 兰州教育学院学报 2012.114.《最小二乘法的创立及其思想方法》 贾小勇,徐传胜,白欣 西北大学学报2006.6。