数据拟合法

第四章 数据拟合法在科学实验和生产实践中，有许多函数关系仅能用由实验或观测得到的一组数据表(,)(0,1,,)i i x y i m =来表示，例如某种物质的化学反应，能够测得生成物的浓度与时间关系的一组数据表.而它们的解析表达式)(t f y =是不知道的。

但是为了要知道化学反应速度，必须要利用已知数据给出它的近似表达式，有了近似表达式，通过求导数便可知道化学反应速度。

可见已知一组数据求它的近似表达式是非常有意义的.如何求它的近似表达式呢？第二章介绍的插值方法是一种有效的方法.但是由于数据(,)(0,1,,)i i x y i m =是由测量或观测得到的，它本身就有误差，作插值时一定要通过型值点),(i i y x 似乎没有必要；其次当m 很大时，采用插值（特别是多项式插值）很不理想（会出现龙格现象），非多项式插值计算又很复杂。

为此，本章介绍一种“整体”近似的方法，即对于给定的数据(,),0,1,,i i x y i n =，选一个线性无关函数系)(,),(),(10x x x n ϕϕϕ ，以它们为基底构成的线性空间为{}0span (),,()n x x ϕϕ=Φ.在此空间内选择函数()()nj j j x x ϕαϕ==∑其中(0,1,,)j j n α=为待定常数。

要求它逼近真实函数)(x f y =的误差尽可能小，这就是数据拟合问题.§1 最小二乘法一、最小二乘法设有数据(,),0,1,,i i x y i m =，令()(),0,1,,ni i i i j j i j r y x y x i m ϕαϕ==-=-=∑.并称Tm r r r r ),,,(10 =为残向量，用)(x ϕ去拟合)(x f y =的好坏问题变成残量的大小问题。

判断残量大小的标准，常用的有下面几种：(1) 确定参数(0,1,,)j j n α=，使残量绝对值中最大的一个达到最小，即i mi r ≤≤0max 为最小。

(2) 确定参数(0,1,,)j j n α=，使残量绝对值之和达到最小，即∑=mi ir为最小。

(3) 确定参数(0,1,,)j j n α=，使残量的平方和达到最小，即2mT ii rr r ==∑ 最小(1)和(2)两个标准很直观，但因为有绝对值，所以实际应用很不方便；而标准(3)既直观，使用又很方便。

按标准(3)确定待定参数，得到近似函数的方法，通常称为最小二乘法。

在实际问题中如何选择基函数()(0,1,,)j x j n ϕ=是一个复杂的问题，一般要根据问题本身的性质来决定。

如果从问题本身得不到这方面的信息，那么通常可取的基函数有多项式、三角函数、指数函数、样条函数等。

下面重点介绍多项式的情况。

设基函数取为()(0,1,,)j j x x j n ϕ==. 已知列表函数()(0,1,,)i i y f x i m ==，且n m . 用多项式01()n n n p x a a x a x =+++ (1.1)去近似)(x f ，问题是应该如何选择n a a a ,,,10 使)(x p n 能较好地近似列表函数)(x f . 按最小二乘法，应选择n a a a ,,,10 使得 2010(,,,)[()()]mn i n i i s a a a f x p x ==-∑ (1.2)取最小。

注意到s 是非负的，且是n a a a ,,,10 的二次多项式，它必有最小值。

求s 对na a a ,,,10 的偏导数，并令其等于零，得到10[()]0,0,1,,mn k ii n i i i y aa x a x x k n =-+++==∑进一步将上式写成如下方程组010002101000012010000(1)()(),()()(),()()(),m m mni i n i i i i m m m mn i i i n i i i i i i m m m mn n n n i i i n i i i i i i m a x a x a y x a x a x a x y x a x a x a x y ===+====+====⎧++++=⎪⎪⎪+++=⎪⎨⎪⎪⎪+++=⎪⎩∑∑∑∑∑∑∑∑∑∑∑ 再将方程组写成矩阵形式20000231000012200001 m m mn i i i i i i m m m m n i i i i i i i i m m m mn n n n i i i i i i i i m x x x a x x x x x x x x ===+====++====⎡⎤+⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦∑∑∑∑∑∑∑∑∑∑∑0100m i i m i i i n m n i ii y x y a a x y ===⎡⎤⎢⎥⎢⎥⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎢⎥⎢⎥⎣⎦∑∑∑. (1.3) 若令200000211111211,,1n n n n m mmm a y x x x a y x x x A Y a y x x x α⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥===⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦则(1.3)可简单地表示为.Y A A A T T =α定义1 方程组(1.4)称为法方程组（也叫正规方程组或正则方程组），而Y A =α（1n +个未知量，1m +个方程式） (1.5) 称为超定方程组（也叫做矛盾方程组）.可以证明α为超定方程组(1.4)的最小二乘解的充分必要条件是α满足(1.3).定理1 法方程组(1.4)有唯一一组解。

定理2 设01,,,m a a a 是法方程组(1.4)的解，则多项式0()mi m i i p x a x ==∑是问题的解。

正规方程组方按下表来构造：试按最小二乘法求)(x f 的二次近似多项式.法方程组为.857.5185.7942.9826.1 090.2 503.2090.2 530.2 250.3503.2 250.35210⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡a a a解得0121.036,0.751,0.928.a a a ===故22() 1.0360.7510.928p x x x =++下表给出了)(在节点处的误差：在利用最小二乘法建立和式(1.2)时，所有点i x 都起到了同样的作用，但是有时依据某种理由认为∑中某些项的作用大些，而另外一些作用小些（例如，一些i y 是由精度高的仪器或由操作上比较熟练的人员获得的，自然应该予比较大的信任），在数学上常表现为用2(()())miinii f x p x ρ=-∑ (1.6)替代(1.2)取最小值，此处诸,0>i ρ且∑==mi i1ρ，并称i ρ为权，而(1.6)称为加权和, 并称)(x p n 为)(x f y =在点集},,{0m x x 上关于权函数}{i ρ的最小二乘逼近多项式。

二、内积表示作(),()f x g x 关于权函数()x ρ及01,,,m x x x 的内积(,)()()()mi i i i f g x f x g x ρ==∑ (1.7)其中权函数()x ρ满足()0,0,1,2,,i x i m ρ>=，以4,2m n ==为例，方程组(1.4)化为000110220000111122110021122222(,)(,)(,)(,),(,)(,)(,)(,),(,)(,)(,)(,)a a a f a a a f a a a f ϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕ++=⎧⎪++=⎨⎪++=⎩ (1.8)其中44(,),,0,1,2;(,),0,1,2j k k j k i ik i i i i x x j k f y x k ϕϕϕ======∑∑，20122,4,()1,()1,(),(),()(0,1,2,3,4)i i n m x x x x x x y f x i ρϕϕϕ==≡=====.用矩阵表示为001020000111201102122022(,)(,)(,)(,)(,)(,)(,)(,)(,)(,)(,)(,)a f a f a f ϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕ⎛⎫⎛⎫⎛⎫ ⎪⎪ ⎪= ⎪⎪ ⎪ ⎪ ⎪⎪⎝⎭⎝⎭⎝⎭ (1.12)例2 已知函数()y f x =的数据为0.20.50.70.8511.221 1.6492.014 2.340 2.718i i x y试用最小二乘法求()f x 的二次近似多项式22012()p x a a x a x =++.解 根据题意，得20122,4,()1,()1,(),(),()(0,1,2,3,4)i i m N x x x x x x y f x i ρϕϕϕ==≡=====01234012340.2,0.5,0.7,0.85,1,1;1.221,1.649,2.014,2.430, 2.718, 2.x x x x x N y y y y y m ============44420010200004442011121000442202122200(,)115,(,)1 3.250,(,)1 2.503,(,)1 3.250,(,) 2.503,(,) 2.090,(,)1 2.503,(,) 2.090,(,i i i i i i i i i i i i i i i i i i x x x x x x x x x x ϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕ=========⨯==⨯==⨯==⨯==⨯==⨯==⨯==⨯=∑∑∑∑∑∑∑∑422044420120) 1.826,(,)19.942,(,)7.185,(,) 5.857,i i i i i i i i i i i x x f y f y x f y x ϕϕϕ=====⨯==⨯==⨯==⨯=∑∑∑∑得法方程组01253.250 2.5039.9423.250 2.503 2.0907.1852.503 2.090 1.826 5.858a a a ⎛⎫⎛⎫⎛⎫ ⎪⎪ ⎪= ⎪⎪ ⎪ ⎪⎪ ⎪⎝⎭⎝⎭⎝⎭解得0121.036, 1.036,0.928.a a a === 于是，所求多项式为22() 1.036 1.0360.928p x x x =++.[注] 1) 实际计算表明：当m 较大时，法方程组(1.8)往往是病态的。

因此提高拟合多项式的次数不一定能改善逼近效果。

实际计算中常采用不同的低次多项式去拟合不同的分段，这种方法称为分段拟合。

2) 如何找到更符合实际情况的数据拟合，一方面要根据专业知识和经验来确定经验曲线的近似公式；另一方面要根据散点图的分布形状及特点来选择适当的曲线取拟合这些数据。