第三章 插值法和最小二乘法 3.7 数据拟合 (最小二乘法 )§ 3.7 数据拟合 (最小二乘法 )§ 实例 ： 考察某种纤维的强度与其拉伸倍数的关系 ,下表 是实际测定的 24个纤维样品的强度与相应的拉伸倍数 是记录 : 编 号 拉伸倍数 强 度 编 号 拉伸倍数 强 度 1 1.9 1.4 13 5 5.5 2 2 1.3 14 5.2 5 3 2.1 1.8 15 6 5.5 4 2.5 2.5 16 6.3 6.4 5 2.7 2.8 17 6.5 6 6 2.7 2.5 18 7.1 5.3 7 3.5 3 19 8 6.5 8 3.5 2.7 20 8 7 9 4 4 21 8.9 8.5 10 4 3.5 22 9 8 11 4.5 4.2 23 9.5 8.1 12 4.6 3.5 24 10 8.1 ii yx ii yx 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 纤维强度随拉伸 倍数增加而增加 系 要关系应是线性关 的主与拉伸倍数 因此可以认为强度 xy 并且 24个点大致分 布在一条直线附近 xxy 10)( bb += 为待定参数其中 10 , bb ---------(1) 越接近越好 样本点与所有的数据点我们希望 ),)(()( 10 ii yxxxy bb += 必须找到一种度量标准来衡量什么曲线最接近所有数据点 一 、 最小二乘法的基本概念 iii yxy ?= )(d令 一般使用 ∑ = = m i i 0 22 2 dd 在回归分析中称为残差 ∑ = ?= m i ii yxy 0 2))(( 准偏离程度大小的度量标与数据点作为衡量 ),()( ii yxxy 称为平方误差 在回归分析中称为残差平方和 从而确定 (1)中的待定系数 ∑ = = m i i 0 22 2 dd ∑ = ?= m i ii yxy 0 2))(( 注意 (1)式是一条直线 关系的关系并不一定是线性但 yx, 因此将问题一般化 )(, xSyyx =的关系为设 Φ来自函数类其中 )(xS 来自线性函数类中如 )()1( xy 为给定的一组数据设 ),,1,0)(,( miyx ii L= ),,1,0)(( nixi L=Φ j的基函数为设函数类 mn ≤一般要求 即生成的函数集是由也称 ,),,1,0)(( nixi L=Φ j )}(,),(),({ 10 xxxspan njjj L=Φ ∑ = = m i i 0 22 2 dd ∑ = ?= m i ii yxS 0 2))(( 仍然定义平方误差 ∑ = = n j jj xaxS 0 )()( j Φ∈ 我们选取的度量标准是 )(* xS中选取一个函数在函数类 Φ ∑ = = n j jj xaxS 0 * )()(* j )(*)(*)(* 1100 xaxaxa nn jjj +++= L 2 2*d ∑ = ?= m i ii yxS 0 2))(*( ∑ =Φ∈ ?= m i iixS yxS 0 2 )( ))((min22 )( min d Φ∈ = xS 中的任意函数为其中 Φ= ∑ = m j jj xaxS 0 )()( j ---------(2) ---------(3) 数据拟合的最小二乘法 的方法为的求函数称满足条件 ∑ = = n j jj xaxS 0 * )()(*)3( j 为最小二乘解∑ = = n j jj xaxS 0 * )()(* j 为拟合系数为拟合函数 ),,1,0(,)()( 0 njaxaxS j n j jj L== ∑ = j ),,1,0(,)( njaxS j L=如何求拟合系数后在确定了拟合函数 呢 ？满足拟合条件使得 )3()()(* 0 *∑ = = n j jj xaxS j 误差称为最小二乘解的平方22*d ∑ ∑ = = ?= m i i n j ijj yxa 0 2 0 ))(( j∑ = ?= m i ii yxS 0 2))(( 二 、 法方程组 2 2d ∑ = = n j jj xaxS 0 )()( j由 的函数为拟合系数 ),,1,0( nja j L= 可知 因此可假设 ),,,( 10 naaa Ly ∑ ∑ = = ?= m i i n j ijj yxa 0 2 0 ))(( j 因此求最小二乘解转化为 二次函数 的问题点极小值的最小值求 *,*,*,)(),,,( 1010 nn aaaaaa LLy 由多元函数取极值的必要条件 0),,,( 10 =?? k n a aaa Ly nk ,,1,0 L= )]())((2[ 0 0 ik m i i n j ijj xyxa jj∑ ∑ = = ?= ka? ?y 0= 得 即 ∑∑∑ == = = m i iki m i ik n j ijj xyxxa 00 0 )()()( jjj 0)]()()([ 0 0 =?∑ ∑ = = ik m i i n j ikijj xyxxa jjj ∑∑∑ == = = m i iki m i ik n j ijj xyxxa 00 0 )()()( jjj ∑∑ ∑ == = = m i iki n j jik m i ij xyaxx 00 0 )()]()([ jjj nk ,,1,0 L= ---------(4) ∑ ∑∑∑ = === = +++ m i iki ik m i innik m i iik m i i xy xxaxxaxxa 0 00 11 0 00 )( )()()()()()( j jjjjjj L nk ,,1,0 L= 即 元线性方程组的是一个关于显然 1,,,)4( 10 +naaa nL 引入记号 ))(,),(),(( 10 mrrr xxx jjj L=rj ),,,( 10 myyy L=f )()(),( 0 ij m i ikjk xx jjjj ∑ = = 则由内积的概念可知 i m i ikk yxf ∑ = = 0 )(),( jj ---------(5) ---------(6) ),( jk jj ),( kj jj=显然内积满足交换律 方程组 (4)便可化为 ),(),(),(),( 1100 faaa knknkk jjjjjjj =+++ L nk ,,1,0 L= ---------(7) 的线性方程组常数项为这是一个系数为 ),(),,( fkjk jjj 将其表示成矩阵形式 ?? ? ? ? ? ? ?? ? ? ? ? ? na a a M 1 0 ?? ? ? ? ? ? ?? ? ? ? ? ? = ),( ),( ),( 1 0 f f f nj j j M ?? ? ? ? ? ?? ? ? ? ? ),(),(),( 01000 njjjjjj L ),(),(),( 11101 njjjjjj L ),(),(),( 10 nnnn jjjjjj L MMM -----(8) baGn =简单记为 上的法方程组在点 式为函数序列称 m n xxx xxx ,,, )(,),(),()8( 10 10 L L jjj 的基为函数类由于 Φ)(,),(),( 10 xxx njjj L 必然线性无关因此 )(,),(),( 10 xxx njjj L 并且其系数矩阵为对称阵 所以法方程组的系数矩阵非奇异 ,即 0])),det[(()det( )1()1( ≠= +×+ nnjinG jj 根据 Cramer法则 ,法方程组有唯一解 *,*,*, 1100 nn aaaaaa === L *),*,*,( 10 naaa Ly ∑ ∑ = = ?= m i i n j ijj yxa 0 2 0 ))(( j),,,( 10 naaa Ly 即 是 的最小值 2 2*d ∑ = ? m i ii yxS 0 2))(*( ∑ =Φ∈ ?= m i iixS yxS 0 2 )( ))((min 2 2)(min dΦ∈= xS 所以 ∑ ∑ = = ? m i i n j ijj yxa 0 2 0 ))(*( j ∑ ∑ = =Φ∈ ?= m i i n j ijjxS yxa 0 2 0)( ))((min j ∑ ∑ = = ?= m i i n j ijj yxa 0 2 0 ))(*( j 为最小二乘解∑ = = n j jj xaxS 0 * )()(* j因此 的拟合函数作为常使用多项式 ),,1,0)(,()()( miyxxPxS iin L== 作为一种简单的情况 , 的基函数为拟合函数 )()( xPxS n= ,1)(0 =xj ,)(1 xx =j ,,)(, LL kk xx =j nn xx =)(j 基函数之间的内积为 )()(),( 0 ij m i ikjk xx jjjj ∑ = = ∑ = = m i j i k i xx 0 ∑ = += m i jk ix 0 i m i ikk yxf ∑ = = 0 )(),( jj ∑ = = m i i k i yx 0 =22*d平方误差 ∑ = ? m i ii yxS 0 2))(*( ∑ = ?= n j jj faff 0 ),(*),( j 例 1. 回到本节开始的实例 ， 从散点图可以看出 纤维强度和拉伸倍数之间近似与线性关系 xaaxy 10)( += 故可选取线性函数 为拟合函数 ， 其基函数为 1)(0 =xj xx =)(1j 建立法方程组 根据内积公式 ， 可得 24),( 00 =jj 5.127),( 10 =jj 61.829),( 11 =jj 1.113),( 0 =fj 6.731),( 1 =fj 法方程组为 ?????? 61.8295.127 5.12724 ?????? 1 0 a a ? ? ?? ? ?= 6.731 1.113 1505.00 =a 即为所求的最小二乘解xxy 8587.01505.0)(* += 8587.01 =a解得 6615.5* 22 =d平方误差为 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 拟合曲线与散点 的关系如右图 : 例 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 cos之间具有三角函数关系与 xexy 系之间还具有指数函数关与 xxy ln系之间还具有对数函数关与 因此假设拟合函数与基函数分别为 xcexbxaxS ++= cosln)( xex =)( 2jxx ln)(0 =j xx cos)(1 =j 0 0.5 1 1.5 2 2.5 3 -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 -49.0086 1002.5 ? ? ? ? ? ?? ? ? ? 1.6163 -2.3827 26.7728?? ? ? ? ?? ? ? ? 通过计算 ,得法方程组的系数矩阵及常数项矩阵为 y Go! 用 Gauss列 主元消去法 ,得 =??? ? ? ?? ? ? ? c b a -1.0410 -1.2613 0.030735?? ? ? ? ?? ? ? ? xexxxS 030735.0cos2613.1ln0410.1)(* +??= 的最小二乘解是关于 xy 2 2*d 2 0 ))(*(∑ = ?= m i ii yxS 2 0 )030735.0cos2613.1ln0410.1(∑ = ?+??= m i i x ii yexx i 92557.0= 拟合的平方误差为 图象如图 例 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 具有图示的图形的曲线很多 ， 本题特提供两种形式 tbaey =指数函数形式 bat ty +=双曲线形式 都是待定系数其中 ba, tbay 1lnln += tbay 11 += tbaey =指数函数形式).1( tbay 1lnln +=两边取对数 ,得 aattyy ln,1,ln =′=′=′设 tbay ′+′=′得 即为拟合函数 基函数为 ,1)(0 =′tj tt ′=′)(1j 0567.1,427.2 ?==′ ba解法方程组得 325.11=a tey 0567.1325.11 ?=最小二乘解为 11631.0* 221 =d平方误差为 bat ty +=双曲线形式).2( tbay 11 += 16272.0080174.0 == ba 用最小二乘法得 即 16272.0080174.0 += t ty 5621.1* 222 =d 无论从图形还是从平方误差考虑 在本例中指数函数拟合比双曲线拟合要好 0 2 4 6 8 10 12 14 16 4 5 6 7 8 9 10 11 t yyyy 平方误差为 三 、 加权最小二乘法 ),,1,0)(,( miyx ii L=对于一组给定的数据点 中在拟合的数据点 ),,1,0)(,( miyx ii L= 各点的重要性可能是不一样的 的重度表示数据点假设 ),( iii yxw 重度 : 即权重或者密度 ， 统称为权系数 mk ,,1,0 L= 定义加权 平方误差为 ∑== m i ii 0 22 2 dwd ∑ = ?= m i iii yxy 0 2))((w -----(9) Φ来自函数类设拟合函数 )(xS ),,1,0)(( nixi L=Φ j的基函数为函数类 )}(,),(),({ 10 xxxspan njjj L=Φ ∑ = ?= m i iii yxS 0 2))(*(w )(xS )()()( 1100 xaxaxa nnjjj +++= L 为拟合系数),,1,0( nja j L= ),,1,0(* nja j L=组拟合的目标仍然为找一 2 2*d ∑ =Φ∈ ?= m i iiixS yxS 0 2 )( ))((min w22 )( min d Φ∈ = xS 使得 ),,,( 10 naaa Ly求 ∑ ∑ = = ?= m i i n j ijji yxa 0 2 0 ))(( jw 的问题点极小值的最小值 *,*,*,)( 10 naaa L 由多元函数取极值的必要条件 0),,,( 10 =?? k n a aaa Ly nk ,,1,0 L= )]())((2[ 0 0 ik m i i n j ijji xyxa jjw∑ ∑ = = ?= ka? ?y 0=得 即 ∑∑ ∑ == = = m i ikii m i ik n j ijji xyxxa 00 0 )()()( jwjjw 0)]()()([ 0 0 =?∑ ∑ = = ik m i ii n j ikijji xyxxa jwjjw ∑∑ ∑ == = = m i ikii m i ik n j ijii xyxxa 00 0 )()()( jwjjw ∑∑ ∑ == = = m i ikii n j jik m i iji xyaxx 00 0 )()]()([ jwjjw nk ,,1,0 L= 元线性方程组的是一个关于显然 1,,,)10( 10 +naaa nL 引入记号 ))(,),(),(( 10 mrrr xxx jjj L=rj ),,,( 10 myyy L=f 定义加权内积 -----(10) )()(),( 0 ij m i ikijk xx jjwjj ∑ = = i m i ikik yxf ∑ = = 0 )(),( jwj ),(),(),(),( 1100 faaa knknkk jjjjjjj =+++ L nk ,,1,0 L= 矩阵形式 (法方程组 )为 ?? ? ? ? ? ? ?? ? ? ? ? ? na a a M 1 0 ?? ? ? ? ? ? ?? ? ? ? ? ? = ),( ),( ),( 1 0 f f f nj j j M ?? ? ? ? ? ?? ? ? ? ? ),(),(),( 01000 njjjjjj L ),(),(),( 11101 njjjjjj L ),(),(),( 10 nnnn jjjjjj L MMM 方程组 (10)式化为 -----(11) ---(12) 平方误差为 ∑ = ?= m i iii yxS 0 2))(*(w2 2*d 作为特殊情形 ,用多项式作拟合函数的法方程组为 ?? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? = ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ∑ ∑ ∑ ∑∑∑ ∑∑∑ ∑∑∑ = = = = + == + === === 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 w w w www www www MM L MMM L L -----(13) 四 、 用正交多项式作最小二乘拟合 * 选为基底的基函数若拟合函数 )()(xS ),()( 00 xPx =j ,L),()( 11 xPx =j )()( xPx nn =j 为正交多项式且 )(,),(),( 10 xPxPxP nL ),,1,0)(,( miyx ii L=对于一组给定的数据点 ),( jk PP ?? ?= jk ≠0 jk =kA ∑ = = m i ijiki xPxP 1 )()(w即 0>kA其中 正交多项式如何选取呢 -----(14) 线性无关显然 )(,),(),( 10 xPxPxP nL 线性表示次多项式均可由任意且 )(,),(),( 10 xPxPxPk kL 1)(,),(),( 10 时令其首项均为选取正交多项式 xPxPxP nL )(1 xxPk k次多项式考虑 + 线性表示显然其可由 )(,),(),( 110 xPxPxP k +L 1,,,,, 110 kk abbb ?L即存在系数 )()()( 1 1 0 xPxPxP kkk k j jj + ? = ++?= ∑ ab 使得 )(xxPk )()()( 1 1 0 xPxPxP kkk k j jj + ? = ++?= ∑ ab)(xxPk ),( sk PxP ))(),()()(( 1 1 0 xPxPxPxP skkk k j jj + ? = ++?= ∑ ab )( ks <),( sk PxP ),( sss PPb= ),( kk PxP ),( kkk PPa= 由 可知 因此 ),( ),( kk kk PP PxP= ka sb ),( ),( ss sk PP PxP= sb ),( ),( ss sk PP PxP=而 ),( ),( ss sk PP xPP= 次多项式为 1)( +sxxPs 线性表示可由正交多项式组 10)}({ +sj xP 时当 ks <+ 1 0),( =sk xPP时即 1?< ks 因此 sb ?? ??? ?= ?< = 1 10 ks ks ),( ),( 11 1 ?? ? kk kk PP xPP ),( ),( 11 1 ?? ? kk kk PP xPP ),( ),( 11 ?? = kk kk PP PP )()()()( 1 0 1 xPxPxxPxP kk k j jjkk ab ???= ∑ ? = + )()()( 11 xPxPx kkkk ????= ba 可知 最后可得正交多项式选取的方法 : 1)(0 =xP 01 )( a?= xxP ∑ =+ = m i iixm 0 0 1 1 wa ),( ),(? 00 00 PP PxP= )(1 xPk + )()()( 11 xPxPx kkkk ????= ba ),( ),( 11 1 ?? ? = kk kk k PP PPb ),( ),( kk kk PP PxP= ka -----(15) ni ,,2,1 L= )()()( 1 1 0 xPxPxP kkk k j jj + ? = ++?= ∑ ab)(xxPk由 ?? ? ? ? ? ? ?? ? ? ? ? ? na a a M 1 0 ?? ? ? ? ? ? ?? ? ? ? ? ? = ),( ),( ),( 1 0 f f f nj j j M ?? ? ? ? ? ?? ? ? ? ? ),(),(),( 01000 njjjjjj L ),(),(),( 11101 njjjjjj L ),(),(),( 10 nnnn jjjjjj L MMM 作拟合选择正交多项式 )(,),(),( 10 xPxPxP nL ),,1,0)(,( miyx iii L=的数据点对于一组给定的带权 w )()()()( 1100 xPaxPaxPaxS nn+++= L ∑ = ?= m i iii yxS 0 2))(*(w2 2*d ∑ =Φ∈ ?= m i iiixS yxS 0 2 )( ))((min w22)(min dΦ∈= xS 使得 由正交多项式的性质 ,法方程组 ),(),( fPaPP kkkk = ni ,,2,1,0 L= -----(16) ),( ),(* kk k k PP fPa = -----(17) ?? ? ? ? ? ? ?? ? ? ? ? ? na a a M 1 0 ?? ? ? ? ? ? ?? ? ? ? ? ? = ),( ),( ),( 1 0 fP fP fP n M ?? ? ? ? ? ?? ? ? ? ? 00),( 00 LPP 0),(0 11 LPP ),(00 nn PPL MMM 可化为 即 得 )(*)(*)(*)(* 1100 xPaxPaxPaxS nn+++= L即 为利用正交多项式的最小二乘解 ∑ = ?= m i iii yxS 0 2))(*(w2 2*d 平方误差为 ))(*,)(*( 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 ),(*),( 例 4. 如下及权重给定数据点 iii yx w),( 1111111 0.371.244.219.296.175.11 0.19.08.07.06.05.00 i i i y x w 是用最小二乘法求拟合这组数据的多项式 解 : 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3 从散点图可知 数据和二次多项式拟合较好 因此选用二次多项式作 这组数据的拟合函数 0.6 0.8 1 设拟合函数 为基函数)(),(),( 210 xPxPxP )()()()( 221100 xPaxPaxPaxS ++= 取 1)(0 =xP 01 )( a?= xxP 0a ),( ),( 00 00 PP PxP= 7 5.4= 642857.0?= x 0a ),( ),( 00 0 PP fP= 15.2 7 05.15 == 1a ),( ),( 11 1 PP fP= 978260.1 657143.0 300000.1 == 1a ),( ),( 11 11 PP PxP= 335403.0 657143.0 220408.0 == ),( ),( 00 11 0 PP PP=b 093878.0 7 657143.0 == )(2 xP )()()( 0011 xPxPx ba ??= 093878.0)642857.0)(335403.0( ???= xx 2a ),( ),( 22 2 PP fP= 999942.0 068660.0 068656.0 == 121738.0978260.02 +?= xx )()()()( 221100 xPaxPaxPaxS ++= 因此拟合多项式为 999993.0000057.1999942.0 2 ++= xx 平方误差为 2 2*d ∑ = ?= n k kkk PPaff 0 2 ),(*),( 1010313856.2 ?×= 五 、 利用正交多项式作最小二乘法的算法设计