《[理学]数值计算第6章曲线拟合的最小二乘法.doc》由会员分享，可在线阅读，更多相关《[理学]数值计算第6章曲线拟合的最小二乘法.doc（30页珍藏版）》请在三一办公上搜索。
1、第6章 曲线拟合的最小二乘法6.1 拟合曲线通过观察或测量得到一组离散数据序列，当所得数据比较准确时，可构造插值函数逼近客观存在的函数，构造的原则是要求插值函数通过这些数据点，即。此时，序列与是相等的。如果数据序列，含有不可避免的误差（或称“噪音”），如图6.1所示；如果数据序列无法同时满足某特定函数，如图6.2所示，那么，只能要求所做逼近函数最优地靠近样点，即向量与的误差或距离最小。按与之间误差最小原则作为“最优”标准构造的逼近函数，称为拟合函数。图6.1 含有“噪声”的数据图6.2 一条直线公路与多个景点插值和拟合是构造逼近函数的两种方法。插值的目标是要插值函数尽量靠近离散点；拟合的目标是
2、要离散点尽量靠近拟合函数。向量与之间的误差或距离有各种不同的定义方法。例如：用各点误差绝对值的和表示：用各点误差按模的最大值表示：用各点误差的平方和表示：或 （6.1）其中称为均方误差，由于计算均方误差的最小值的方法容易实现而被广泛采用。按均方误差达到极小构造拟合曲线的方法称为最小二乘法。本章主要讲述用最小二乘法构造拟合曲线的方法。在运筹学、统计学、逼近论和控制论中，最小二乘法都是很重要的求解方法。例如，它是统计学中估计回归参数的最基本方法。关于最小二乘法的发明权，在数学史的研究中尚未定论。有材料表明高斯和勒让德分别独立地提出这种方法。勒让德是在1805年第一次公开发表关于最小二乘法的论文，这
3、时高斯指出，他早在1795年之前就使用了这种方法。但数学史研究者只找到了高斯约在1803年之前使用了这种方法的证据。在实际问题中，怎样由测量的数据设计和确定“最贴近”的拟合曲线？关键在选择适当的拟合曲线类型，有时根据专业知识和工作经验即可确定拟合曲线类型；在对拟合曲线一无所知的情况下，不妨先绘制数据的粗略图形，或许从中观测出拟合曲线的类型；更一般地，对数据进行多种曲线类型的拟合，并计算均方误差，用数学实验的方法找出在最小二乘法意义下的误差最小的拟合函数。例如，某风景区要在已有的景点之间修一条规格较高的主干路，景点与主干路之间由各具特色的支路联接。设景点的坐标为点列；设主干路为一条直线，即拟合函
4、数是一条直线。通过计算均方误差最小值而确定直线方程（见图6.2）。6.2线性拟合和二次拟合函数线性拟合给定一组数据，做拟合直线，均方误差为 （6.2）是二元函数，的极小值要满足整理得到拟合曲线满足的方程： （6.3）或 称式（6.3）为拟合曲线的法方程。用消元法或克莱姆法则解出方程：a=例6.1 下表为P. Sale及R. Dybdall在某处作的鱼类抽样调查，表中为鱼的数量，为鱼的种类。请用线性函数拟合鱼的数量和种类的函数关系。13151621222325293031361110111212131312141617404255606264707210013013142214212124172
5、334解：设拟合直线，并计算得下表：编号xyxyx2123452113151621221309561110111212343441431501762522644420189131692252564414841690061640将数据代入法方程组（6.3）中，得到：解方程得：= 8.2084，= 0.1795拟合直线为：= 8.2084 + 0.1795二次拟合函数给定数据序列，用二次多项式函数拟合这组数据。设，作出拟合函数与数据序列的均方误差：（6.4）由多元函数的极值原理，的极小值满足整理得二次多项式函数拟合的法方程：（6.5）解此方程得到在均方误差最小意义下的拟合函数。方程组（6.5）称为
6、多项式拟合的法方程，法方程的系数矩阵是对称的。当拟保多项式阶时，法方程的系数矩阵是病态的，在计算中要用双精度或一些特殊算法以保护解的准确性。例6.2 给定一组数据，如下表。用二次多项式函数拟合的这组数据。32101234230125解：设，由计算得下表：32101234230125112430141539941014928368301845727810182708116101168196将数据代入式（6.5），相应的法方程为：解方程得：=0.66667，=1.39286，=0.13095= 0.666671.392860.13095拟合曲线的均方误差：=3.09524结果见图6.3。图6.3拟
7、合曲线与数据序6.3解矛盾方程组在6.2节中用最小二乘法构造拟合函数，本节中用最小二乘法求解线性矛盾方程组的方法构造拟合函数。给定数据序列，作拟合多项式，如果要求过这些点，那么有矛盾方程组： 即： (6.6)我们作一辅助函数这是自变量为的多元函数，要使 达到最小值，采用多元函数求权值的方法， 对每一个自变量的偏导数为0。整理成以 为未知数的线性方程组整理成矩阵形式，其中：这是一个的对称方程组，称为法方程。只要非奇异，就可以得出唯一解。这就是矛盾方程组的最小二乘解。有什么快捷的方法来求法方程的解？把矛盾方程组(6.6)写成矩阵形式，其中容易验证 ，法方程就是。例如，拟合直线的矛盾方程组的形式如下
8、：化简得到与（6.3）相同的法方程：在线性代数中，我们知道，关于方程组，若秩秩，则方程组无解，这时方程组称为矛盾方程组。在数值代数中对矛盾方程计算的是在均方误差极小意义下的解，也就是在最小二乘法意义下的矛盾方程的解。定理6.1将证明，方程组的解就是矛盾方程组在最小二乘法意义下的解。定理6.11为行列的矩陈，为列向量，秩，称为矛盾方程的的法方程，法方程恒有解。2是的解，当且仅当满足，即是法方程的解。证明：1）对作行初等变换，使， 秩而秩秩秩 方程组有解而且解惟一存在。2）设满足，任取，则由于是任取的，故法方程组的解为极小问题的解。事实上，对离散数据作次多项式曲线拟合，要计算的极小问题。这与解矛盾
9、方程组或求的极小问题是一回事。在这里故对离散数据所作的次拟合曲线，可通过解下列方程组求得：例6.3 给定一组数据，见下表，用二次多项式函数拟合这组数据。32101234230125解：记二次拟合曲线为 ，形成法方程： 而=得到：解方程得：= 0.66667，= 1.39286，= 0.13095 = 0.666671.39286 0.13095例6.4 给出一组数据，见下表。用最小二乘法求形如的经验公式。x32124y14.38.34.78.322.7解：列出法方程： =法方程为：解方程得到：=10.675，= 0.137拟合曲线为：10.678+0.137 有些非线性函数经过转换以后可化为线
10、性函数计算。例如，令，则化拟合曲线为：。例6.5 求一个形如的经验函数公式，使它能够和下列数据相拟合。1234567815.320.527.436.649.165.687.8117.6解：化经验公式为线性形式，对经验公式的两边取自然对数有由矛盾方程组得到法方程组即 解方程组得：拟合曲线的均方误差为：计算结果见图6.4。图6.4 拟合曲线图示例6.6 解矛盾方程组解：写出法方程组 ，即得解法方程得：1.5917，0.5899，0.7572。6.4权有的实际问题中，已知数据不一定都是一次观测的结果，对于不同的可能观测次数不同，在矛盾方程组中，一组确定一个方程，而最小二乘解对每个方程来说都还存在误差
11、，这样，把每对都同等看待就不太合理，希望观测次数多的（即可靠性大些的）数据在矛盾方程组中占的比重大些。为了使问题的提法更有一般性，通常把最小二乘法中偏差平方和改为加权平方和，使最小，其中称作权。同样，可由误差函数求极值，得到法方程，其中，事实上，只需对矛盾方程组作一些修改，每个方程乘上，得到，其中这样，就使，仍然可用前面的方法，对，左乘用求得最小二乘解。6.5用正交多项式作最小二乘拟合我们不仅可用多项式来拟合函数，还可用一般的函数来拟合。定义6.1 若，如果当且仅当时成立，则称在上线性无关，称为上线性无关族。最小二乘法的一般提法为：对给定的一组数据，要求在函数类中找一个函数，使加权平地均和，其
12、中这里是线性无关的函数族，是上的权函数，点处的表示该点数据的重要程度。求误差函数的极小值点，由多元函数极值的必要条件得：这是个未知数个方程的方程组，称为法方程式。定义6.2：称为与关于点集的内积。这样，法方程式可简写为，记为，其中称为克莱姆行列式，记作。定理6.2 的充要条件是线性无关。证明：若存在使对此式两边分别取与的内积得：这是一个以为未知数的齐次方程组，有非零解的充要条件是系数矩阵行列式等于零，于是的充要条件是方程有全零解，即全为0，所以线性无关。证毕。由于法方程有惟一解的充要条件是，因而线性无关也是法方程有惟一解的充要条件。特别当取为时，由于是在中的线性无关函数族，因而必有最小二乘解。
13、用上的多项式拟合，需要解一个的线性方程组，且当取得大一此时，法方程的系数矩阵会出现病态。从系数矩阵B的形式看，里面的元素都是些内积，是否能取某些函数族，使对非对角元素全变为0？如果有这样的函数族，那么方程容易解，病态也得到改善。定义6.3 函数族如果在点集上满足称为点集带权的正交函数族。例6.7 三角函数族上是正交函数族（权）.实际上，而如果拟合函数在上取，且是正交函数族，则法方程式成为：直接可得到，不用解线性方程组了。且容易估算，是否病态也容易判断。完成以上工作的关键在于如何构造正交函数族。正交多项式是最简单的正交函数族。常用的正交多项式如：Chebyshev（切比雪夫）多项式、Legend
14、re（勒让德）多项式等。现在我们根据给定结点及权函数，造出带权正交的多项式族，用递推公式表示如下：其中这样给出的是正交的，这一点可以用归纳法证明。例6.8 已知函数表，利用正交多项式求拟合多项。123441018261111解：设所以：所以：所以：得：6.6程序示例程序6.1 用线性函数 拟合给定数据。算法描述输入值，及。解方程组：输出。程序源码/purpose：（x_i，y_i）线性拟合函数 /# include stdio.h# define MAX_N 25 /（x_i，y_i）的最大维数typedef struct tagPOINT / 点的结构 double x；double y；
15、POINT；int main ( ) int m； int i； POINT points MAX_N ； static double u11，u12，u21，u22，c1，c2； double a，b，tmp； printf (“ nInput n value：”)；/ 输入点数m scanf (“% d”，&m)；if (mMAX_N ) printf (“The input n is larger then MAX_N，please redefine the MAX_N. n”)； return 1；if (m= 0) printf (“Please input a number bet
16、ween 1 and % d. n”，MAX_N )；return 1；printf (“Now input the（x_i，y_i），i=0，% d: n”，m1)；/ 输入点for (i=0；i= m；i+) scanf (“% lf”，&tmp)；printf i.x=tmp；scanf (“%lf ”，& tmp)； printf i.y=tmp；/ scanf (“%lf %lf ”，& printf i.x，& printf i.y)；for (i=0；i= n；i+) / 列出方程U(a，b)T = cu21+ = points i.x；u22+ = points i.xpoin
17、ts i.x； c1+ = points i.y； c2+ = points i.ypoints i.y；u12 = u21；u11 = m； / 求解a = ( clu22c2u12) / (u11 u22u12u21)；b = ( cl *u21c2u11) / (u21 u12u22u11)；printf(“Solve：p(x) = %f +%f x n”，a，b)； / 输出return 0； / End of File计算实例用线性函数拟合下列数据（例6.1）。13151621222325293031361110111212131312141617404255606264707210
18、013013142214212124172334 程序输入输出input m value：21Now input the (x_i，y_i)，i=0，20：13 11 15 10 16 11 21 12 22 12 23 13 25 13 29 1230 14 31 16 36 17 40 13 42 14 55 22 60 14 62 2164 21 70 24 72 17 100 23 130 34Solve：p(x) = 8.208408 + 0.179522 x程序6.2 用形如的函数拟合定数据。算法描述输入值，及。解方程组：输出所求拟合函数 。程序源码/purpose：（x_i，y_
19、i）形如y=aExp(bx)的拟合函数 /# include stdio.h# include math.h# define MAX_N 25 /（x_i，y_i）的最大维数typedef struct tagPOINT / 点的结构 double x；double y； POINT；int main ( ) int m；int i；POINT points MAX_N ； static double u11，u12，u21，u22，c1，c2； double A，B，tmp； printf (“ nInput n value：”)； / 输入点数mscanf (“% d”，&m)；if (m
20、MAX_N ) printf (“The input n is larger then MAX_N，please redefine the MAX_N. n”)； return 1；if (m= 0） printf (“Please input a number between 1 and % d. n”，MAX_N )；return 1；printf (“Now input the（x_i，y_i），i=0，% d: n”，m1)；/ 输入点for (i=0；i= m；i+） scanf (“% lf ”，&tmp)； printf i.x = tmp；scanf (“%lf ”，& tmp
21、)； printf i.y = tmp；/ scanf (“%lf %lf ”，& printf i.x，& printf i.y)；for (i=0；i= n；i+) / 列出方程X(a，b)T = Y / 其中Y为对原方程取对数后的数据u21+ = points i.x；u22+ = points i.xpoints i.x； c1+ = log (points i.)y； c2+ = points i.ylog ( points i.y)； u12 = u21；u11 = m； / 求解a = ( clu22c2u12) / (u11 u22u12u21)；b = ( clu21c2u1
22、1) / (u21 u12u22u11)；printf(“Solve：p(x) = % lfexp (%lf x) n”，exp(A)，B)； / 输出return 0； / End of File计算实例求一个形如的经验函数公式，使它能够和下列数据相拟合（例6.5）。1234567815.320.527.436.649.165.687.8117.6程序输入输出input m value：8Now input the（x_i，y_i），i=0，7：1 15.3 2 20.5 3 27.4 4 36.6 5 49.1 6 65.6 7 87.8 8 117.6Solve：p(x) = 11.437069exp(0.291215x)
