最小二乘曲线拟合及其MATLAB实现

1、 现代测量数据处理方法学生课题论文论 文 题 目 ：最小二乘曲线拟合及其MATLAB实现学 院 ：土木工程学院年级专业班 ：2013级测绘工程一班学 生 姓 名： 学 生 学 号：指 导 老 师提 交 时 间：2016年1月成 绩教师签名目 录0 引 言31 曲线拟合与最小二乘法概述41.1 曲线拟合简介41.2 最小二乘法简介52 曲线拟合的最小二乘法原理62.1 原理的阐述及理论公式推导62.2 结合实例分析与理解82.3 总结归纳求解步骤113 基于MATLAB的最小二乘曲线拟合123.1 MATLAB软件介绍123.2 求解的基本理论阐述133.3 结合实例进行MATLAB解算144

2、最小二乘曲线拟合案例分析与解算164.1 案例叙述164.2 数据输入与分析174.3 进行拟合求解184.3.1 手工解算184.3.1 基于MATLAB的解算194.4 拟合函数的精度检测214.5拟合函数在实际运用中的优势225 结 论23参考文献24最小二乘曲线拟合及其MATLAB实现陈涛1（1. 重庆交通大学土木工程学院，重庆；）摘 要随着人类认识能力的不断进步以及计算技术的快速发展，对于变量之间的未知关系，应用曲线拟合的方法对揭示其内在规律具有重要的理论与现实意义。在科学实验数据的处理、分析时，实验数据拟合是经常采用的一种方法。本文将采用最小二乘法对给定的实验数据进行拟合并得到拟合

3、曲线，加深大家对最小二乘曲线拟合原理的理解。同时将根据最小二乘拟合理论，并利用MATLAB数值分析软件进行编程，解决最小二乘曲线拟合在塔机起重量监测系统中的应用问题，实现相应案例数据的曲线拟合，获得了曲线模型对相应数据的拟合曲线，很好地解决了该工程案例的曲线拟合问题。关键词：曲线拟合，最小二乘法，MATLAB0 引 言在科学实验的统计方法研究中，往往要从一组数据中，寻找自变量x与因变量y之间的函数关系。由于观测数据往往不准确，因此不要求经过所有带点，而只要求在给定点上的误差按某种标准最小。若记，就是要求向量的范数最小。如果用最大范数，计算上困难较大，通常采用欧式范数作为误差度量的标准。的函数类

4、型往往与实验的物理背景以及数据的实际分布有关，它一般含有某些待定参数。如果是所有待定参数的线性函数，那么相应的问题称为线性最小二乘问题，否则称为非线性最小二乘问题。最小二乘法还是实验数据参数估计的重要工具。这是因为这种方法比其他方法更容易理解，即使在其他方法失效的情况下，用最小二乘法还能提供解答，而且从统计学的观点分析，用该方法求得各项估计具有最优统计特征，因此这一方法也是系统识别的重要基础。用最小二乘法求拟合曲线时，首先要确定的形式，然后利用最小二乘曲线拟合去构造一个近似解析式。利用该方法“拟合”出的函数曲线虽然不能保证通过所有的样本点，但是很好地“逼近”了它们，充分反映了已知数据间内在的数

5、量关系。因此，这种方法在生产实践和科学实验中具有广泛的应用前景。本文针对最小二乘曲线拟合的有关理论和应用问题以及相应的MATLAB实现进行探讨。1 曲线拟合与最小二乘法概述1.1 曲线拟合简介实际工作中，变量间未必都有线性关系，如服药后血药浓度与时间的关系；疾病疗效与疗程长短的关系；毒物剂量与致死率的关系等常呈曲线关系。曲线拟合是指选择适当的曲线类型来拟合观测数据，并用拟合的曲线方程分析两变量间的关系。曲线拟合是用连续曲线近似地刻画或比拟平面上离散点组所表示的坐标之间的函数关系的一种数据处理方法。用解析表达式逼近离散数据的一种方法。在科学实验或社会活动中，通过实验或观测得到量x与y的一组数据对

6、（i1，2，m），其中各是彼此不同的 。人们希望用一类与数据的背景材料规律相适应的解析表达式，来反映量x与y之间的依赖关系，即在一定意义下“最佳”地逼近或拟合已知数据。常称作拟合模型 ，在式中是一些待定参数。当c在中线性出现时，称为线性模型，否则称为非线性模型。有许多衡量拟合优度的标准，最常用的一种做法是选择参数c使得拟合模型与实际观测值在各点的残差(或离差) 的加权平方和达到最小，此时所求曲线称作在加权最小二乘意义下对数据的拟合曲线。有许多求解拟合曲线的成功方法，对于线性模型一般通过建立和求解方程组来确定参数，从而求得拟合曲线。至于非线性模型，则要借助求解非线性方程组或用最优化方法求得所需参

7、数才能得到拟合曲线，有时称之为非线性最小二乘拟合。1.2 最小二乘法简介最小二乘法是法国大数学家A.M.Legendre 最先于1805年发表的，其动机是为处理一类从天文学和测地学中提出的数据分析问题。它通过最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以简便地求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方和为最小。最小二乘法还可用于曲线拟合，工程施工中，我们会经常取得一些相关的数据，这些数据往往来自与施工密切相关的测量或实验中，我们可以通过作图或多段插值取得变量之间的联系，但作图和插值查图往往误差较大。这时可采用最小二乘法先拟合出一个多项式，再根据此多项式求解任一自

8、变量所对应的因变量较精确的结果，据此绘图可得到较精确、较合理的曲线。1801年，意大利天文学家朱赛普皮亚齐发现了第一颗小行星谷神星2,经过40天的跟踪观测后，由于谷神星运行至太阳背后，使得皮亚齐失去了谷神星的位置。随后全世界的科学家利用皮亚齐的观测数据开始寻找谷神星，但是根据大多数人计算的结果来寻找谷神星都没有结果。时年24岁的高斯也计算了谷神星的轨道。奥地利天文学家海因里希奥尔伯斯根据高斯计算出来的轨道重新发现了谷神星。高斯使用的最小二乘法的方法发表于1809年他的著作天体运动论中。 法国科学家勒让德于1806年独立发现“最小二乘法”，但因不为世人所知而默默无闻。 勒让德曾与高斯为谁最早创立

9、最小二乘法原理发生争执。 1829年，高斯提供了最小二乘法的优化效果强于其他方法的证明，因此被称为高斯-莫卡夫定理。最小二乘法（又称最小平方法）是一种数学优化技术。它通过最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以简便地求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方和为最小。2 曲线拟合的最小二乘法原理2.1 原理的阐述及理论公式推导给定数据，设拟合函数形式为 2.11其中为已知的线性无关函数（如果存在不全为零的常熟，使得，则称函数线性相关，否则称为线性无关）。求系数，使得 2.12最小，若 2.13则称相应的为最小二乘拟合曲线。特别的，若则称为次最小二乘拟合多项

10、式。下面用求多元函数极值的方法来求最小点。将（2.12）式两边对求偏导。并令 化简得 2.14为了进一步化简，可以引入内积符号。在线性代数中，中两个向量及的内积定义为，将它加以推广，得到下面结论：设与是两个已知函数，记，令利用内积的定义，式（2.14）可以写为 2.15其中， ， 2.16方程组（2.15）称为正规方程组或法方程组，其中系数矩阵是对称的。可以证明，当函数线性无关时，方程组（2.15）是对称正定的，因此有唯一解。求出方程组（2.15）的解后，代入式（2.11）即可得最小二乘拟合函数。另外，对带权的最小二乘拟合函数有如如下的定义：设，给定在个节点上的函数值及一组权系数，若有函数，满

11、足则称为在个节点上关于权系数的最小二乘拟合函数。2.2 结合实例分析与理解Intel公司董事长Moore在上个世纪的60年代就观察到一个很有趣的现象：集成电路上可容纳的单晶体数量每隔一年半左右并会增长一倍，从而使集成电路的性能也能提高一倍。据此他提出了轰动世界的Moore定律，预测这种增长趋势会一直延续下去。下面给出Moore数据，如表 1所示：（年）19591962196319641965（增长倍数）13456表 1 Moore数据画出相应的散点图如图 1所示： 图 1 Moore数据散点图表 1中第二行数据为芯片上晶体数目在不同年代与1959年时的数目比较的倍数，通过观察k与t中间大致呈线

12、性关系，如图 1所示。据此导出了著名的Moore定律。通过以上的分析，可设 2.21将表 1中的数据代入式（2.21）的超定方程组，其中，t表示时间，k表示增长倍数，a,b为待定系数。若将表 1中的数据代入式（2.21），得线性方程组 2.22方程组（2.22）是一个朝顶方程组，在这五个线性方程中，任意两个联立求解可得到十组不同的解。即是说该方程组不存在通常意义上的解。现将线性方程组（2.22）写出矩阵形式,其中，此超定方程组五常义解，即是说不存在使得，但是该超定方程组存在最小二乘解，也就是说存在，使得达到最小，并且是线性方程组 2.23的解。我们称式（2.23）为法方程组，在本例中它是一个二

13、阶线性方程组，即解这个方程组得.由此得到Moore公式.需要说明的是，对于，显然，但是根据曲线拟合的最小二乘原理，从整体趋势上使偏差达到最小，此处的偏差，这个值已经很小了、满足要求。2.3 总结归纳求解步骤下面我们就以上摩尔（Moore）预测公式实例总结利用最小二乘曲线拟合原理求解实际问题的步骤：（1）分析数据，根据散点图设定拟合函数（2）代入数据得到超定方程组，该超定方程组的矩阵形式为，其中，.（3）如表 2所示，建立法方程组.1959119591962358861963478521964598201965611790表 2 据表 2中计算结果得, 其中m为实测数据组数。（4）解法方程组得拟

14、合参数向量并据此得到拟合曲线函数（5）通过将所得的拟合函数曲线与原始数据散点图进行同坐标对比或计算总体趋势上的偏差值检验拟合函数的精度。3 基于MATLAB的最小二乘曲线拟合3.1 MATLAB软件介绍MATLAB是matrix和laboratory两个词的组合，意为矩阵工厂（矩阵实验室）。是由美国mathworks公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上

15、摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式，代表了当今国际科学计算软件的先进水平。MATLAB和Mathematica、Maple并称为三大数学软件。它在数学类科技应用软件中在数值计算方面首屈一指。MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。MATLAB的基本数据单位是矩阵，它的指令表达式与数学、工程中常用的形式十分相似，故用MATLAB来解算问题要比用C，FORTRAN等语言完成相同的事情简捷得多，并且MATLAB也吸收了像Ma

16、ple等软件的优点，使MATLAB成为一个强大的数学软件。在新的版本中也加入了对C，FORTRAN，C+，JAVA的支持。可以直接调用,用户也可以将自己编写的实用程序导入到MATLAB函数库中方便自己以后调用，此外许多的MATLAB爱好者都编写了一些经典的程序，用户可以直接进行下载就可以用。 3.2 求解的基本理论阐述假设有一组数据，且已知这组数据满足某一函数原型，其中为待定系数向量，则最小二乘曲线拟合的目标就是求出这一组待定系数值，使得目标函数最小。MATLAB的统计工具箱提供了Isqcurvefit函数，可以解决最小二乘曲线拟合问题。 该函数的调用格式如下：其中，Fun为原始函数的MATL

17、AB表示，可以是M-函数或inline函数；为最优化的初值；x，y为原始输入输出向量。调用该函数，将返回待定系数向量，以及在此系数下的目标函数的值Jm.。3.3 结合实例进行MATLAB解算此处我们就结合上面的Moore实例进行分析，通过对散点图的分析我们已经假设出了初始函数.其实现的MATLAB程序如下：t=1959 1962 1963 1964 1965;k=1 3 4 5 6;令，这样，原函数就可以写出，可以用MATLAB程序代码写出：function k = K( a,t )k=a(1)+a(2)*t; end %定义原型函数kformat long%小数精度定义为小数点后15位a,J

18、M=lsqcurvefit(K,1;1,t,k); %调用Isqcurvefit函数求系数和偏差值结果如下：aa = 1.0e+03 * -1.1662 0.2693JMJM = 0.2990由结果可知a=a(1)=-1625.,b=a(2)=0.,即由MATLAB解算出来的拟合函数为这与上面手工解算的结果基本一致。另外，我们可以将原始数据的散点图和得到的拟合曲线画在一个坐标画面上以检测拟合函数的精度，编写程序如下：ti=1959 1962 1963 1964 1965;ki=-1625.+0.*ti;plot(t,k,o,ti,ki)结果如图 2所示，最小二乘法曲线拟合的结果是找到符合经验公

19、式的最优曲线，但这一经验公式是否有效还需要事后检验，一般就是从图像上做出判断。定量的方法也是有的，一般是计算残差平方和，再进行统计检验，对此就不做多余的讲解了。通过对图 2的分析可知拟合曲线与原始数据是比较稳合的，满足要求。图 2 拟合函数对原始数据的逼近4 最小二乘曲线拟合案例分析与解算4.1 案例叙述起重量限制器是用来保护塔机的重要装置之一，是用于防止因超重而引起起升电机、传动机构、钢丝绳的损坏。但是它只能在极限的状态下保护塔机起升机构不会受到损坏，不能够显示起重量值，因此，司机在操作过程中不了解塔机每次起吊重量的具体状况。为了进一步提高塔机的安全性能和工作效率，增加塔机起重量在线监测装置

20、非常重要，实时准确地测量出起重量是在线监测的关键。在实时临测系统中，在原有起重量限制器的基础上加装了拉力传感器，传感器所测量的拉杆拉力Q与钢丝绳的张力F之间存在着一定的函数关系，起重量增加，拉杆拉力也相应增加，因此可通过间接测量拉杆拉力的方法先测出钢丝绳的张力，然后根据吊钩处的钢丝绳倍率关系计算出实际起重量，从而在拉杆拉力与塔机起重量之间建立起函数关系。4.2 数据输入与分析由于塔机起重量G与钢丝绳张力F之间有确定的函数关系，在实际应用中，以塔机起重量G代替钢丝绳张力F作为输出样本，以拉杆拉力Q作为输入样本。塔机QTZ63最大额定起重量为，分别以，为起重量，测量相应的拉力传感器拉力Q，以获取样

21、本表 3所示。样本拉杆拉力Q/kN起重量G/kg10020.4560030.94120041.44180052.10240062.61300073.36360084.27420095.164800106.055400117.336000表 3实测样本、估算值及相对误差 利用MATLAB画出其散点图，分析其函数模型，程序代码如下：G=0 600 1200 1800 2400 3000 3600 4200 4800 5400 6000;Q=0 0.45 0.94 1.44 2.10 2.61 3.36 4.27 5.16 6.05 7.33;plot(Q,G,o)结果如图 3所示，根据散点图的走势

22、我们可以设原函数为三次多项式函数模型： 图 3样本数据散点图4.3 进行拟合求解 由于该案例的运算量不是太大，所以在这里我们在采用MATLAB解算的同时也进行了一次手工解算，同时也可以就两种方法的对比体现出MATLAB解算的高效、简单与快捷的特点，但是在这个大数据时代，我们碰到的更多是手工解算所不能完成的大数据，到时候就只能利用MATLAB实现了。4.3.1 手工解算（1）由原始数据散点图得出函数模型如下（2）将表 3中的原始数据代入上式得超定方程组 该超定方程组的矩阵形式为,其中（3）建立法方程组(4)解法方程组，得所以所得的拟合函数为4.3.1 基于MATLAB的解算通过对散点图的分析我们

23、已经假设出了函数原型.求解过程实现的MATLAB程序如下：G=0 600 1200 1800 2400 3000 3600 4200 4800 5400 6000;Q=0 0.45 0.94 1.44 2.10 2.61 3.36 4.27 5.16 6.05 7.33;令，这样，原函数就可以写出，用MATLAB程序代码写出：function G = F( a,Q )G=a(1)*Q.3+a(2)*Q.2+a(3)*Q; end%定义原型函数Gformat long%小数精度定义为小数点后15位a,JM=lsqcurvefit(F,1;1;1,Q,G); %调用Isqcurvefit函数求系数

24、和偏差值结果如下：a = 1.0e+03 * 0.4153 -0.6377 1.2993 JM = 1.6404e+04由结果可知a=a(1)=4.6019,b=a(2)=-110.5342,c=a(3)=1383.5979。即由MATLAB解算出来的拟合函数为拟合函数与上面手工解算的结果是一致的。通过传统手工解算与MATLAB解算过程的对比，能够明显的看出MATLAB解算比手工解算要简单方便的多，这种对比在数据量更大的案例中会更加的显著，而且在那种大数据处理中手工解算是很容易出错的，MATLAB解算可以避免这种错误。4.4 拟合函数的精度检测由于手工解算和MATLAB解算的拟合结果是一样的，

25、我们采用画图法检测精度是只需要画一幅图如图 4所示：图 4拟合函数对原始数据的逼近观察散点图可知拟合函数曲线与原始数据的吻合度是非常高的，定性分析的整体精度是满足要求的。不同之处在于这里MATLAB解算的偏差值从表面上看起来是很大的，这主要是因为该案例采用的数据原本就大，我们可以求其相对误差，从另一方面定量来看一下其拟合精度，相对误差表如下表 4所示：样本拉杆拉力Q/kN起重量G/kg拟合函数所求起重量G/kg相对误差/%1000020.45600600.660.1530.9412001206.740.5741.4418001776.921.3052.1024002460.722.4762.6

26、130002940.042.0473.3636003575.570.6984.2742004250.881.2095.1648004828.570.59106.0554005344.011.05117.3360006015.270.25表 4 拟合曲线所求起重量相对理论值的相对误差表观察该表其单个样本的相对误差最大也就为2.47%，从事件发生的概率来讲这个概率的事件属于小概率事件，另外 GB5l4494规定：起重机应安装起重量限制器，对最大起重量大于6t的起重机如设有提示装置，则其数值误差不得大于指示值的5，因此定量分析结果的精度是满足要求的。4.5拟合函数在实际运用中的优势塔机起重量监测中存

27、在的非线性问题中，采用数据拟合理沦，建立了起重量 G和拉杆拉力Q之间的函数关系式，使塔机起重量监测在PLC中得以实现。在实际运用中，该方法具有如下优点：(1)计算结果惟一，计算量小，便于在PLC、单片机等硬件设备上实现；(2)可精确、方便地实现起重量的实时监测；(3)当钢丝绳倍率改变时，只需调整对应多项式的系数，不必改动其它硬件设施；(4)保留了原有起重量限制器中的超重预警开关和超重报警开关，能够实现起重量预警和报警的双重保护。5 结 论当今最小二乘法已经广泛的应用于各类学科，成为了不可缺少的重要工具。目前在物理学、地质勘探学、概率论、统计学等领域有着重要的应用。而最小二乘法曲线拟合的出现，又

28、使得图像呈现更加直观，程序代码简单，使用方便，已经成为研究人员开展科研工作的有效工具之一。在做完了这篇论文后，学习到了许多新的知识。对于最小二乘法有了深一步的认识，了解了它的计算原理以及对于现在的测量估算上的意义，并对MATLAB也有了重新的认识，感受到了MATLAB在现代数据处理中的重要地位。最小二乘法如果想将曲线拟合的比较完美，必须应用适当的模拟曲线，如果模拟曲线选择不够适当，那么用最小二乘法计算完后，会发现拟合曲线误差比较大，均方误差也比较大，而如果拟合曲线选择适当，那么效果较好，例如在本文中，不论是Moore定律案例还是塔机的起重量监测案例其实最初的曲线模型确定都是非常重要而且并不是那

29、么容易的，只是本文侧重于利用理论求解问题而忽视了对这一点的讲解，不同的曲线模型得出的结果是完全不一样的。因此，在实际应用中需要对已知点根据分布规律选取多个可能的近似拟合曲线，算出后比较误差与均方误差，得到最佳拟合曲线。但是如果已知点分布非常不规律，无法观察或是无法正确观察出其近似曲线，那么根本无法使用最小二乘法进行曲线拟合，我们只能使用其它方法进行逼近。 参考文献1 胡晓东，董辰辉.MATLAB从入门到精通M.一.北京:人民邮电出版社, 2010.xiao-dong hu, Dong Chen fai. MATLAB from entry to master M. Beijing: peopl

30、es posts and telecommunications publishing house, 2010. 2 薛莲.数值计算方法M.一.北京:电子工业出版社, 2007.Xue Lian. Numerical calculation method M. Beijing: electronic industry press, 2007. 3 张德丰.MATLAB概率与数理统计分析M.一.北京:机械工业出版社, 2010.de-feng zhang. MATLAB probability and mathematical statistics analysis M. Beijing: mec

31、hanical industry publishing house, 2010. 4 张光澄，张雷.实用数值分析M.二.四川:四川大学出版社, 2009.Zhang Guangcheng, zhang lei. Practical numerical analysis M. 2. Sichuan: sichuan university press, 2009. 5 武汉大学测绘学院测量平差学科组.误差理论与测量平差基础M.三.武汉:武汉大学出版社, 2014.wuhan university institute of surveying and mapping surveying adjust

32、ment discipline groups. Error theory and measurement adjustment M. 3. Wuhan: wuhan university press, 2014. 6 盛骤，谢式千，潘承毅.概率论与数理统计M.四.北京:高等教育出版社, 2013.fill flash, xie qian, Pan Chengyi. Probability and mathematical statistics M. Beijing: higher education press, 2013. 7 同济大学数学系.工程数学线性代数M.五.北京:高等教育出版社,

33、2007.tongji university mathematics. Engineering mathematics, linear algebra M. 5. Beijing: higher education press, 2007. 8 北京信息科技大学.最小二乘法原理，VC+实现及应用EB/OL.2015-12-15./view/f35289b202d276ae91.html?from=search.Beijing university of science and technology information. The least squa

34、res principle, application and vc + + realize EB/OL. 2015-12-15. /view/f35289b202d276ae91.html? The from = search. 9 燕山大学.数据拟合的几个应用实例 EB/OL.2015-12-15.http:/ /view/dce54b64caaedd3383c4d36f.html?from=searchSeveral applications in the yanshan university. The data fitting EB/OL.2015-12-15. /view/dce54b64caaedd3383c4d36f.html? The from = search English TitleThe principle