曲线拟合的最小二乘法

1、1 2009, Henan Polytechnic University 1 实例：考察某种纤维的强度与其拉伸倍数的关系实例：考察某种纤维的强度与其拉伸倍数的关系, ,下下 表是实际测定的表是实际测定的2424个纤维样品的强度与相应的拉伸倍个纤维样品的强度与相应的拉伸倍 数是记录数是记录: : 编 号 拉伸倍数 强 度编 号 拉伸倍数 强 度 11.91.41355.5 221.3145.25 32.11.81565.5 42.52.5166.36.4 52.72.8176.56 62.72.5187.15.3 73.531986.5 83.52.72087 944218.98.5 1043.

2、52298 114.54.2239.58.1 124.63.524108.1 ii yx ii yx 2 2009, Henan Polytechnic University 2 12345678910 1 2 3 4 5 6 7 8 9 12345678910 1 2 3 4 5 6 7 8 9 纤维强度随拉伸 倍数增加而增加 系系 要要关关系系应应是是线线性性关关 的的主主与与拉拉伸伸倍倍数数 因因此此可可以以认认为为强强度度 xy 并且24个点大致分 布在一条直线附近 为待定参数为待定参数其中其中 10, xxy 10 )( 越接近越好越接近越好 样本点样本点与所有的数据点与所有的数据点

3、我们希望我们希望),)()( 10ii yxxxy 必须找到一种度量标准来衡量什么曲线最接近所有数据点。必须找到一种度量标准来衡量什么曲线最接近所有数据点。 3 2009, Henan Polytechnic University 3 iii yxy )( 令令 一般使用一般使用 在回归分析中称为在回归分析中称为 m i i 0 2 2 2 m i ii yxy 0 2 )( 准准偏离程度大小的度量标偏离程度大小的度量标与数据点与数据点作为衡量作为衡量),()( ii yxxy 在回归分析中称为残差平方和在回归分析中称为残差平方和. . 从而确定从而确定(1)(1)中的待定系数中的待定系数:

4、: m i i 0 2 2 2 m i ii yxy 0 2 )( 注意注意(1)(1)式是一条直线式是一条直线, ,关关系系的的关关系系并并不不一一定定是是线线性性但但yx, 因此将问题一般化为因此将问题一般化为: : 4 2009, Henan Polytechnic University 4 )(,xSyyx 的的关关系系为为设设 ,)( 来自函数类来自函数类其中其中xS来自线性函数类来自线性函数类中中如如)()1(xy 为给定的一组数据为给定的一组数据设设), 1 , 0)(,(miyx ii ), 1 , 0)(nix i 的基函数为的基函数为设函数类设函数类 mn 一一般般要要求求

5、 即即生成的函数集生成的函数集是由是由也称也称,), 1 , 0)(nix i )(,),(),( 10 xxxspan n m i i 0 2 2 2 m i ii yxS 0 2 )( 仍然定义平方误差仍然定义平方误差 n j jj xaxS 0 )()( 5 2009, Henan Polytechnic University 5 我们选取的度量标准是我们选取的度量标准是 )(* xS中中选选取取一一个个函函数数在在函函数数类类 n j jj xaxS 0 * )()(* )(*)(*)(* 1100 xaxaxa nn 2 2 * m i ii yxS 0 2 )(*( m i ii

6、xS yxS 0 2 )( )(min 2 2 )( min xS 中中的的任任意意函函数数为为其其中中 m j jj xaxS 0 )()( 6 2009, Henan Polytechnic University 6 数据拟合的最小二乘法数据拟合的最小二乘法 的方法为的方法为的求函数的求函数称满足条件称满足条件 n j jj xaxS 0 * )()(*)3( 为为最最小小二二乘乘解解 n j jj xaxS 0 * )()(* 为为拟拟合合系系数数为为拟拟合合函函数数), 1 , 0(,)()( 0 njaxaxS j n j jj ), 1 , 0(,)(njaxS j 如如何何求求拟

7、拟合合系系数数后后在在确确定定了了拟拟合合函函数数 满满足足拟拟合合条条件件呢呢？使使得得 n j jj xaxS 0 * )()(* 误误差差称称为为最最小小二二乘乘解解的的平平方方 2 2 * 7 2009, Henan Polytechnic University 7 m i i n j ijj yxa 0 2 0 )( m i ii yxS 0 2 )( 2 2 n j jj xaxS 0 )()(由由 的的函函数数为为拟拟合合系系数数), 1 , 0(nja j 可知可知 因此可假设因此可假设 ),( 10n aaa m i i n j ijj yxa 0 2 0 )( 因此求最小二

8、乘解转化为因此求最小二乘解转化为 二次函数 8 2009, Henan Polytechnic University 8 的问题的问题点点极小值极小值的最小值的最小值求求*,*,*,)(),( 1010nn aaaaaa 由多元函数取极值的必要条件由多元函数取极值的必要条件 0 ),( 10 k n a aaa nk, 1 , 0 )()(2 00 ik m i i n j ijj xyxa k a 0 得得 即即 m i iki m i ik n j ijj xyxxa 000 )()()( 0)()()( 00 ik m i i n j ikijj xyxxa 9 2009, Henan

9、Polytechnic University 9 m i iki m i ik n j ijj xyxxa 000 )()()( m i iki n j jik m i ij xyaxx 000 )()()( nk, 1 , 0 m i iki ik m i innik m i iik m i i xy xxaxxaxxa 0 00 11 0 00 )( )()()()()()( nk, 1 , 0 即即 10 2009, Henan Polytechnic University 10 元线性方程组元线性方程组的的显然上式是一个关于显然上式是一个关于1, 10 naaa n 引入记号引入记号)

10、(,),(),( 10mrrr xxx r ),( 10m yyyf )()(),( 0 ij m i ikjk xx 则由内积的概念可知则由内积的概念可知 i m i ikk yxf 0 )(),( ),( jk ),( kj 显然内积满足交换律显然内积满足交换律 11 2009, Henan Polytechnic University 11 方程组便可化为方程组便可化为 ),(),(),(),( 1100 faaa knknkk nk, 1 , 0 的的线线性性方方程程组组常常数数项项为为这这是是一一个个系系数数为为),(),(f kjk 将其表示成矩阵形式将其表示成矩阵形式 n a a

11、 a 1 0 ),( ),( ),( 1 0 f f f n ),(),(),( 01000n ),(),(),( 11101n ),(),(),( 10nnnn baGn简单记为 12 2009, Henan Polytechnic University 12 上上的的法法方方程程组组在在点点 列列称称上上述述方方程程组组为为函函数数序序 m n xxx xxx , )(,),(),( 10 10 的基的基为函数类为函数类由于由于 )(,),(),( 10 xxx n 必然线性无关必然线性无关因此因此)(,),(),( 10 xxx n 并且其系数矩阵为对称阵并且其系数矩阵为对称阵 所以法方

12、程组的系数矩阵非奇异所以法方程组的系数矩阵非奇异, ,即即 0),det()det( )1()1( nnjin G 根据根据Cramer法则法则, ,法方程组有唯一解法方程组有唯一解 *,*,*, 1100nn aaaaaa 13 2009, Henan Polytechnic University 13 *),*,*,( 10n aaa m i i n j ijj yxa 0 2 0 )(),( 10n aaa 即即 是是 的最小值 2 2 * m i ii yxS 0 2 )(*( m i ii xS yxS 0 2 )( )(min 2 2 )( min xS 所以所以 m i i n

13、j ijj yxa 0 2 0 )(*( m i i n j ijj xS yxa 0 2 0 )( )(min m i i n j ijj yxa 0 2 0 )(*( 为为最最小小二二乘乘解解 n j jj xaxS 0 * )()(* 因此因此 14 2009, Henan Polytechnic University 14 的拟合函数的拟合函数作为作为常使用多项式常使用多项式), 1 , 0)(,()()(miyxxPxS iin 作为一种简单的情况, 的基函数为的基函数为拟合函数拟合函数)()(xPxS n , 1)( 0 x,)( 1 xx ,)(, k k xx n n xx )

14、( 基函数之间的内积为基函数之间的内积为 )()(),( 0 ij m i ikjk xx m i j i k i xx 0 m i jk i x 0 i m i ikk yxf 0 )(),( m i i k i yx 0 2 2 * 平方误差 m i ii yxS 0 2 )(*( n j jj faff 0 ),(*),( 15 2009, Henan Polytechnic University 15 例例1 1. 回到本节开始的实例，从散点图可以看出 纤维强度和拉伸倍数之间近似与线性关系 xaaxy 10 )( 故可选取线性函数 为拟合函数，其基函数为 1)( 0 xxx )( 1

15、建立法方程组 根据内积公式，可得 16 2009, Henan Polytechnic University 16 24),( 00 5 .127),( 10 61.829),( 11 1 .113),( 0 f 6 .731),( 1 f 法方程组为法方程组为 61.8295 .127 5 .12724 1 0 a a 6 .731 1 .113 1505. 0 0 a 即即为为所所求求的的最最小小二二乘乘解解xxy8587. 01505. 0)(* 8587. 0 1 a解得解得 6615. 5* 2 2 平方误差为平方误差为 17 2009, Henan Polytechnic Univ

16、ersity 17 12345678910 1 2 3 4 5 6 7 8 9 12345678910 1 2 3 4 5 6 7 8 9 拟合曲线与散点拟合曲线与散点 的关系如右图的关系如右图: : 18 2009, Henan Polytechnic University 18 例例2 2. 求拟合下列数据的最小二乘解求拟合下列数据的最小二乘解 x=.24 .65 .95 1.24 1.73 2.01 2.23 2.52 2.77 2.99 y=.23 -.26 -1.10 -.45 .27 .10 -.29 .24.56 1 解解: :从数据的散点图可以看出从数据的散点图可以看出 xxy

17、cos之之间间具具有有三三角角函函数数关关系系与与 x exy系系之间还具有指数函数关之间还具有指数函数关与与 xxyln系系之之间间还还具具有有对对数数函函数数关关与与 因此假设拟合函数与基函数分别为因此假设拟合函数与基函数分别为 x cexbxaxScosln)( x ex )( 2 xxln)( 0 xxcos)( 1 19 2009, Henan Polytechnic University 19 00.511.522.53 -1.5 -1 -0.5 0 0.5 1 x y 6.7941 -5.3475 63.2589 -5.3475 5.1084 -49.0086 63.2589 -

18、49.0086 1002.5 1.6163 -2.3827 26.7728 通过计算通过计算, ,得法方程组的系数矩阵及常数项矩阵为得法方程组的系数矩阵及常数项矩阵为 00.511.522.53 -1.5 -1 -0.5 0 0.5 1 x y Go!Go! 20 2009, Henan Polytechnic University 20 用用Gauss列主元消去法列主元消去法, ,得得 c b a -1.0410 -1.2613 0.030735 x exxxS030735. 0cos2613. 1ln0410. 1)(* 的的最最小小二二乘乘解解是是关关于于xy 2 2 * 2 0 )(*

19、( m i ii yxS 2 0 )030735. 0cos2613. 1ln0410. 1( m i i x ii yexx i 92557. 0 拟合的平方误差为拟合的平方误差为 图象如图图象如图 21 2009, Henan Polytechnic University 21 例例3.3. 在某化学反应里在某化学反应里, ,测得生成物浓度测得生成物浓度y%y%与时间与时间t t的的 数据如下，试建立数据如下，试建立y y关于关于t t的经验公式的经验公式 t=1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16 y=4.00,6.40,8.00,8.80,9.22

20、,9.50,9.70,9.86,10.00, 10.20,10.32,10.42,10.50,10.55,10.58,10.60 解解: : 的的散散点点图图与与浓浓度度画画出出时时间间yt 具有图示的图形的曲线很多，本题特提供两种形式具有图示的图形的曲线很多，本题特提供两种形式 t b aey 指数函数形式指数函数形式 bat t y 双双曲曲线线形形式式 都都是是待待定定系系数数其其中中ba, t bay 1 lnln t ba y 11 22 2009, Henan Polytechnic University 22 t b aey 指数函数形式指数函数形式).1( t bay 1 ln

21、ln两边取对数两边取对数, ,得得 aa t tyyln, 1 ,ln 设设 t bay 得得即为拟合函数即为拟合函数 基函数为基函数为 , 1)( 0 t tt ) ( 1 0567. 1,427. 2ba解法方程组得解法方程组得325.11a t ey 0567. 1 325.11 最小二乘解为最小二乘解为 11631. 0* 2 2 1 平方误差为平方误差为 23 2009, Henan Polytechnic University 23 bat t y 双双曲曲线线形形式式).2( t ba y 11 16272. 0080174. 0ba 用最小二乘法得用最小二乘法得 即即 1627

22、2. 0080174. 0 t t y 5621. 1* 2 2 2 无论从图形还是从平方误差考虑无论从图形还是从平方误差考虑 在本例中指数函数拟合比双曲线拟合要好在本例中指数函数拟合比双曲线拟合要好 0246810121416 4 5 6 7 8 9 10 11 t y 0246810121416 4 5 6 7 8 9 10 11 t y 0246810121416 4 5 6 7 8 9 10 11 t y 0246810121416 4 5 6 7 8 9 10 11 t y 平方误差为平方误差为 24 2009, Henan Polytechnic University 24 从本例

23、看到，拟合曲线的数学模型并不从本例看到，拟合曲线的数学模型并不 是一开始就能选好的，往往要通过分析是一开始就能选好的，往往要通过分析 确定若干模型之后，再经过实际计算，确定若干模型之后，再经过实际计算， 才能选到较好的模型。才能选到较好的模型。 25 2009, Henan Polytechnic University 25 ), 1 , 0)(,(miyx ii 对于一组给定的数据点对于一组给定的数据点 中中在拟合的数据点在拟合的数据点), 1 , 0)(,(miyx ii 各点的重要性可能是不一样的 的重度的重度表示数据点表示数据点假设假设),( iii yx 重度: 即权重或者密度，统称

24、为权系数 mk, 1 ,0 定义加权 平方误差为: m i ii 0 2 2 2 m i iii yxy 0 2 )( 26 2009, Henan Polytechnic University 26 来自函数类来自函数类设拟合函数设拟合函数)(xS ), 1 , 0)(nix i 的基函数为的基函数为函数类函数类 )(,),(),( 10 xxxspan n m i iii yxS 0 2 )(*( )(xS)()()( 1100 xaxaxa nn 为为拟拟合合系系数数), 1 , 0(nja j ), 1 , 0(*nja j 组组拟拟合合的的目目标标仍仍然然为为找找一一 2 2 * m

25、 i iii xS yxS 0 2 )( )(min 2 2 )( min xS 使得使得 27 2009, Henan Polytechnic University 27 ),( 10n aaa 求求 m i i n j ijji yxa 0 2 0 )( 的问题的问题点点极小值极小值的最小值的最小值*,*,*,)( 10n aaa 由多元函数取极值的必要条件由多元函数取极值的必要条件 0 ),( 10 k n a aaa nk, 1 , 0 )()(2 00 ik m i i n j ijji xyxa k a 0得得 即即 m i ikii m i ik n j ijji xyxxa 0

26、00 )()()( 0)()()( 00 ik m i ii n j ikijji xyxxa 28 2009, Henan Polytechnic University 28 m i ikii m i ik n j ijii xyxxa 000 )()()( m i ikii n j jik m i iji xyaxx 000 )()()( nk, 1 , 0 元线性方程组元线性方程组的的是一个关于是一个关于显然显然1,)10( 10 naaa n 引入记号引入记号)(,),(),( 10mrrr xxx r ),( 10m yyyf 定义加权内积定义加权内积 29 2009, Henan

27、Polytechnic University 29 )()(),( 0 ij m i ikijk xx i m i ikik yxf 0 )(),( ),(),(),(),( 1100 faaa knknkk nk, 1 , 0 矩阵形式矩阵形式( (法方程组法方程组) )为为 n a a a 1 0 ),( ),( ),( 1 0 f f f n ),(),(),( 01000n ),(),(),( 11101n ),(),(),( 10nnnn 方程组化为方程组化为 30 2009, Henan Polytechnic University 30 平方误差为平方误差为 m i iii yx

28、S 0 2 )(*( 2 2 * 作为特殊情形作为特殊情形, ,用多项式作拟合函数的法方程组为用多项式作拟合函数的法方程组为 m i i n ii m i iii m i ii n i m i i n i m i i n i m i i n i m i ii m i ii m i i n i m i ii m i i m i i yx yx y a a a xxx xxx xx 0 0 0 1 0 2 0 1 00 1 0 2 00 000 31 2009, Henan Polytechnic University 31 选为选为基底基底的基函数的基函数若拟合函数若拟合函数)()(xS ),(

29、)( 00 xPx ,),()( 11 xPx )()(xPx nn 为正交多项式为正交多项式且且)(,),(),( 10 xPxPxP n ), 1 , 0)(,(miyx ii 对于一组给定的数据点对于一组给定的数据点 ),( jk PP jk 0 jk k A m i ijiki xPxP 1 )()( 即即 0 k A其中其中 正交多项式如何选取呢正交多项式如何选取呢 -(14) 32 2009, Henan Polytechnic University 32 0( ) 1P x 110 ( )()( )P xxP x 0 (,)()()() m kjikiji i xxx 1( )k

30、 Px 11 ()( )( ) kkkk xP xPx 11 (,) (,) kk k kk P P PP (,) (,) kk kk xP P P P 1k 1,2,1kn 0,1,2,1kn 33 2009, Henan Polytechnic University 33 n a a a 1 0 ),( ),( ),( 1 0 f f f n ),(),(),( 01000n ),(),(),( 11101n ),(),(),( 10nnnn 作拟合作拟合选择正交多项式选择正交多项式)(,),(),( 10 xPxPxP n ), 1 , 0)(,(miyx iii 的数据点的数据点对于一

31、组给定的带权对于一组给定的带权 )()()()( 1100 xPaxPaxPaxS nn m i iii yxS 0 2 )(*( 2 2 * m i iii xS yxS 0 2 )( )(min 2 2 )( min xS 使得使得 由正交多项式的性质由正交多项式的性质, ,法方程组法方程组 34 2009, Henan Polytechnic University 34 ),(),(fPaPP kkkk ni,2 , 1 , 0 ),( ),( * kk k k PP fP a n a a a 1 0 ),( ),( ),( 1 0 fP fP fP n 00),( 00 PP 0),(

32、0 11 PP ),(00 nn PP 可化为可化为 即即 得得 )(*)(*)(*)(* 1100 xPaxPaxPaxS nn 即即 为利用正交多项式的最小二乘解为利用正交多项式的最小二乘解 35 2009, Henan Polytechnic University 35 m i iii yxS 0 2 )(*( 2 2 * 平方误差为平方误差为 )(*,)(*(fxSfxS ),()*,(2*)*,(fffSSS ),(),(*2),(* 00 2 fffPaPPa n k kk n k kkk ),(),(*2),(* 0 2 0 2 ffPPaPPa n k kkk n k kkk

33、n k kkk PPaff 0 2 ),(*),( 36 2009, Henan Polytechnic University 36 例例4.4.如下如下及权重及权重给定数据点给定数据点 iii yx ),( 1111111 0 . 371. 244. 219. 296. 175. 11 0 . 19 . 08 . 07 . 06 . 05 . 00 i i i y x 是用最小二乘法求拟合这组数据的多项式是用最小二乘法求拟合这组数据的多项式 解解: : 00.10.20.30.40.50.60.70.80.91 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3 从散

34、点图可知从散点图可知 数据和二次多项式拟合较好数据和二次多项式拟合较好 因此选用二次多项式作因此选用二次多项式作 这组数据的拟合函数这组数据的拟合函数 00.10.20.30.40.50.60.70.80.91 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3 37 2009, Henan Polytechnic University 37 设拟合函数设拟合函数 为基函数为基函数)(),(),( 210 xPxPxP )()()()( 221100 xPaxPaxPaxS 取取1)( 0 xP 01 )( xxP 0 ),( ),( 00 00 PP PxP 7 5

35、. 4 642857. 0 x 0 a ),( ),( 00 0 PP fP 15. 2 7 05.15 1 a ),( ),( 11 1 PP fP 978260. 1 657143. 0 300000. 1 38 2009, Henan Polytechnic University 38 1 ),( ),( 11 11 PP PxP 335403. 0 657143. 0 220408. 0 ),( ),( 00 11 0 PP PP 093878. 0 7 657143. 0 )( 2 xP)()()( 0011 xPxPx 093878. 0)642857. 0)(335403. 0(

36、xx 2 a ),( ),( 22 2 PP fP 999942. 0 068660. 0 068656. 0 121738. 0978260. 0 2 xx 39 2009, Henan Polytechnic University 39 )()()()( 221100 xPaxPaxPaxS 因此拟合多项式为因此拟合多项式为 999993. 0000057. 1999942. 0 2 xx 平方误差为平方误差为 2 2 * n k kkk PPaff 0 2 ),(*),( 10 10313856. 2 40 2009, Henan Polytechnic University 40 *)

37、,*,*,( 10n aaa m i i n j ijj yxa 0 2 0 )(),( 10n aaa 即即 是是 的最小值 2 2 * m i ii yxS 0 2 )(*( m i ii xS yxS 0 2 )( )(min 2 2 )( min xS 所以所以 m i i n j ijj yxa 0 2 0 )(*( m i i n j ijj xS yxa 0 2 0 )( )(min m i i n j ijj yxa 0 2 0 )(*( 为为最最小小二二乘乘解解 n j jj xaxS 0 * )()(* 因此因此 41 2009, Henan Polytechnic University 41 00.511.522.53 -1.5 -1 -0.5 0 0.5 1 x y 6.7941 -5.3475 63.2589 -5.3475 5.1084 -49.0086 63.25