数值计算方法37

1、ii i j jiji i l xlb x 1 1 nnnn n n aaa aaa aaa A 21 22221 11211 bAx ni, 3 ,2 第二章 插值与逼近 2.9 数据拟合数据拟合(最小二乘法最小二乘法) 华长生制作2 2.9 数据拟合数据拟合(最小二乘法最小二乘法) 实例：考察某种纤维的强度与其拉伸倍数的关系,下表 是实际测定的24个纤维样品的强度与相应的拉伸倍数 是记录: 编 号 拉伸倍数 强 度编 号 拉伸倍数 强 度 11.91.41355.5 221.3145.25 32.11.81565.5 42.52.5166.36.4 52.72.8176.56 62.72.

2、5187.15.3 73.531986.5 83.52.72087 944218.98.5 1043.52298 114.54.2239.58.1 124.63.524108.1 ii yx ii yx 华长生制作3 12345678910 1 2 3 4 5 6 7 8 9 12345678910 1 2 3 4 5 6 7 8 9 纤维强度随拉伸 倍数增加而增加 系 要关系应是线性关 的主与拉伸倍数 因此可以认为强度 xy 并且24个点大致分 布在一条直线附近 xxy 10 )( 为待定参数其中 10 , -(1) 华长生制作4 越接近越好 样本点与所有的数据点我们希望),)()( 10i

3、i yxxxy 必须找到一种度量标准来衡量什么曲线最接近所有数据点 一、最小二乘法的基本概念 iii yxy)(令 一般使用 m i i 0 2 2 2 在回归分析中称为残差 m i ii yxy 0 2 )( 准偏离程度大小的度量标与数据点作为衡量),()( ii yxxy 称为平方误差 华长生制作5 在回归分析中称为残差平方和 从而确定(1)中的待定系数 m i i 0 2 2 2 m i ii yxy 0 2 )( 注意(1)式是一条直线 关系的关系并不一定是线性但yx, 因此将问题一般化 华长生制作6 )(,xSyyx的关系为设 来自函数类其中)(xS来自线性函数类中如)()1(xy

4、为给定的一组数据设), 1 , 0)(,(miyx ii ), 1 , 0)(nix i 的基函数为设函数类 mn 一般要求 即生成的函数集是由也称,), 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 )()( 华长生制作7 我们选取的度量标准是 )(* xS中选取一个函数在函数类 n j jj xaxS 0 * )()(* )(*)(*)(* 1100 xaxaxa nn 2 2 * m i ii yxS 0 2 )(*( m i ii xS yxS 0 2

5、 )( )(min 2 2 )( min xS 中的任意函数为其中 m j jj xaxS 0 )()( -(2) -(3) 华长生制作8 数据拟合的最小二乘法 的方法为的求函数称满足条件 n j jj xaxS 0 * )()(*)3( 为最小二乘解 n j jj xaxS 0 * )()(* 为拟合系数为拟合函数), 1 ,0(,)()( 0 njaxaxS j n j jj ), 1 ,0(,)(njaxS j 如何求拟合系数后在确定了拟合函数 呢？满足拟合条件使得)3()()(* 0 * n j jj xaxS 误差称为最小二乘解的平方 2 2 * 华长生制作9 m i i n j i

6、jj 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 )( 因此求最小二乘解转化为 二次函数 华长生制作10 的问题点极小值的最小值求*,*,*,)(),( 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

7、 ijj xyxxa 000 )()()( 0)()()( 00 ik m i i n j ikijj xyxxa 华长生制作11 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 -(4) m i iki ik m i innik m i iik m i i xy xxaxxaxxa 0 00 11 0 00 )( )()()()()()( nk, 1 , 0 即 华长生制作12 元线性方程组的是一个关于显然1,)4( 10 naaa n 引入记号)(,),(),

8、( 10mrrr xxx r ),( 10m yyyf )()(),( 0 ij m i ikjk xx 则由内积的概念可知 i m i ikk yxf 0 )(),( -(5) -(6) ),( jk ),( kj 显然内积满足交换律 华长生制作13 方程组(4)便可化为 ),(),(),(),( 1100 faaa knknkk nk, 1 , 0 -(7) 的线性方程组常数项为这是一个系数为),(),(f kjk 将其表示成矩阵形式 n a a a 1 0 ),( ),( ),( 1 0 f f f n ),(),(),( 01000n ),(),(),( 11101n ),(),(),

9、( 10nnnn -(8) 华长生制作14 上的法方程组在点 式为函数序列称 m n xxx xxx , )(,),(),()8( 10 10 的基为函数类由于)(,),(),( 10 xxx n 必然线性无关因此)(,),(),( 10 xxx n 并且其系数矩阵为对称阵 所以法方程组的系数矩阵非奇异,即 0),det( nnji 根据Cramer法则,法方程组有唯一解 *,*,*, 1100nn aaaaaa 华长生制作15 *),*,*,( 10n aaa m i i n j ijj yxa 0 2 0 )(),( 10n aaa 即 是 的最小值 2 2 * m i ii yxS 0

10、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 * )()(*因此 华长生制作16 的拟合函数作为常使用多项式), 1 , 0)(,()()(miyxxPxS iin 作为一种简单的情况, 的基函数为拟合函数)()(xPxS n , 1)( 0 x,)( 1 xx ,)(, k k xx n n xx )( 基函数

11、之间的内积为 )()(),( 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 ),(*),( 华长生制作17 例1. 回到本节开始的实例，从散点图可以看出 纤维强度和拉伸倍数之间近似与线性关系 xaaxy 10 )( 故可选取线性函数 为拟合函数，其基函数为 1)( 0 xxx )( 1 建立法方程组 根据内积公式，可得 华长生制作18 24),( 00 5 .127),( 10 61.

12、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 平方误差为 华长生制作19 12345678910 1 2 3 4 5 6 7 8 9 12345678910 1 2 3 4 5 6 7 8 9 拟合曲线与散点 的关系如右图: 华长生制作20 例2. 求拟合下列数据的最小二乘解 x=.24 .65 .95 1.24 1.73 2.

13、01 2.23 2.52 2.77 2.99 y=.23 -.26 -1.10 -.45 .27 .10 -.29 .24.56 1 解: 从数据的散点图可以看出 xxycos之间具有三角函数关系与 x exy系之间还具有指数函数关与 xxyln系之间还具有对数函数关与 因此假设拟合函数与基函数分别为 x cexbxaxScosln)( x ex )( 2 xxln)( 0 xxcos)( 1 华长生制作21 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

14、 -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! 华长生制作22 用Gauss列主元消去法,得 c b a -1.0410 -1.2613 0.030735 x exxxS030735. 0cos2613. 1ln0410. 1)(* 的最小二乘解是关于xy 2 2 * 2 0 )(*( m i ii yxS 2 0 )030735. 0cos2613. 1ln0410. 1( m i i x ii yexx i 92557.

15、0 拟合的平方误差为 图象如图 华长生制作23 例3.在某化学反应里,测得生成物浓度y%与时间t的 数据如下，试建立y关于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,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

16、 11 华长生制作24 t b aey 指数函数形式).1( t bay 1 lnln两边取对数,得 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 平方误差为 华长生制作25 bat t y 双曲线形式).2( t ba y 11 16272. 0080174. 0ba 用最小二乘法得 即 16272. 0080174. 0 t t y 5621. 1* 2 2 2 无论从图形还是从

17、平方误差考虑 在本例中指数函数拟合比双曲线拟合要好 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 平方误差为 华长生制作26 三、加权最小二乘法 ), 1 , 0)(,(miyx ii 对于一组给定的数据点 中在拟合的数据点), 1 , 0)(,(miyx ii 各点的重要性可能是不一样的 的重度表示数据点假设),( iii yx 重度: 即权重或者密度，

18、统称为权系数 mk, 1 ,0 定义加权 平方误差为 m i ii 0 2 2 2 m i iii yxy 0 2 )(-(9) 华长生制作27 来自函数类设拟合函数)(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 i iii xS yxS 0 2 )( )(min 2 2 )( min xS 使得 华长生制作28 ),( 10n aaa求

19、 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 000 )()()( 0)()()( 00 ik m i ii n j ikijji xyxxa 华长生制作29 m i ikii m i ik n j ijii xyxxa 000 )()()( m i ikii n j jik m i iji xya

20、xx 000 )()()( nk, 1 , 0 元线性方程组的是一个关于显然1,)10( 10 naaa n 引入记号)(,),(),( 10mrrr xxx r ),( 10m yyyf 定义加权内积 -(10) 华长生制作30 )()(),( 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 ),(),(),( 10n

21、nnn 方程组(10)式化为 -(11) -(12) 华长生制作31 平方误差为 m i iii yxS 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 -(13) 华长生制作32 四、用正交多项式作最小二乘拟合* 选为基底的基函数若拟合函数)()(x

22、S ),()( 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) 华长生制作33 线性无关显然)(,),(),( 10 xPxPxP n 线性表示次多项式均可由任意且)(,),(),( 10 xPxPxPk k 1)(,),(),( 10 时令其首项均为选取正交多项式xPxPxP n )(1xxPk k 次多项式考虑 线性

23、表示显然其可由)(,),(),( 110 xPxPxP k 1 , 110kk 即存在系数 )()()( 1 1 0 xPxPxP kkk k j jj 使得 )(xxP k 华长生制作34 )()()( 1 1 0 xPxPxP kkk k j jj )(xxP k ),( sk PxP )(),()()( 1 1 0 xPxPxPxP skkk k j jj )(ks ),( sk PxP),( sss PP ),( kk PxP),( kkk PP 由 可知 因此 ),( ),( kk kk PP PxP k s ),( ),( ss sk PP PxP 华长生制作35 s ),( ),

24、( ss sk PP PxP 而 ),( ),( ss sk PP xPP 次多项式为1)(sxxP s 线性表示可由正交多项式组 1 0 )( s j xP 时当ks10),( sk xPP时即1 ks 因此 s 1 10 ks ks ),( ),( 11 1 kk kk PP xPP ),( ),( 11 1 kk kk PP xPP ),( ),( 11 kk kk PP PP 华长生制作36 )()()()( 1 0 1 xPxPxxPxP kk k j jjkk )()()( 11 xPxPx kkkk 可知 最后可得正交多项式选取的方法: 1)( 0 xP 01 )( xxP m

25、i iix m 0 0 1 1 ),( ),( 00 00 PP PxP )( 1 xP k )()()( 11 xPxPx kkkk ),( ),( 11 1 kk kk k PP PP ),( ),( kk kk PP PxP k -(15) ni,2 , 1 )()()( 1 1 0 xPxPxP kkk k j jj )(xxP k 由 华长生制作37 n a a a 1 0 ),( ),( ),( 1 0 f f f n ),(),(),( 01000n ),(),(),( 11101n ),(),(),( 10nnnn 作拟合选择正交多项式)(,),(),( 10 xPxPxP n

26、 ), 1 , 0)(,(miyx iii 的数据点对于一组给定的带权 )()()()( 1100 xPaxPaxPaxS nn m i iii yxS 0 2 )(*( 2 2 * m i iii xS yxS 0 2 )( )(min 2 2 )( min xS 使得 由正交多项式的性质,法方程组 华长生制作38 ),(),(fPaPP kkkk ni,2 , 1 , 0 -(16) ),( ),( * kk k k PP fP a -(17) n a a a 1 0 ),( ),( ),( 1 0 fP fP fP n 00),( 00 PP 0),(0 11 PP ),(00 nn PP 可化为 即 得 )(*)(*)(*)(* 1100 xPaxPaxPaxS nn 即 为利用正交多项式的最小二乘解 华长生制作39 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 n k kkk PPaff 0 2 ),(*),( 华长生制作40 例4. 如下及权重给定数据点 iii yx),( 1111111 0 . 371.