最小二乘法曲線擬合原理及matlab實(shí)現(xiàn)

1、曲線擬合( curve-fitting )： 工程實(shí)踐中，用測(cè)量到的一些離散的數(shù)據(jù)( X, yj,i 0,1,2,.m求一個(gè)近似的函數(shù)(x)來(lái)擬合這組數(shù)據(jù)，要求所得的擬合曲 線能最好的反映數(shù)據(jù)的基本趨勢(shì)(即使 (x)最好地逼近f x，而不必滿足插值 原則。因此沒(méi)必要取 (X)=yi，只要使i (Xi) yi盡可能地小)。原理：給定數(shù)據(jù)點(diǎn) ( xi,yi),i0 ,1 , 2, . . . m 。求近似曲線 ( x) 。并且使得近似曲線與 f x 的偏差最小。近似曲線在該點(diǎn)處的偏差 i(xi ) yi ， i=1,2,.,m 。常見(jiàn)的曲線擬合方法：1. 使偏差絕對(duì)值之和最小2. 使偏差絕對(duì)值最大

2、的最小3. 使偏差平方和最小最小二乘法：按偏差平方和最小的原則選取擬合曲線， 并且采取二項(xiàng)式方程為擬合曲線的方法 ,稱(chēng)為最小 二乘法。推導(dǎo)過(guò)程：1. 設(shè)擬合多項(xiàng)式為：2. 各點(diǎn)到這條曲線的距離之和，即偏差平方和如下：3問(wèn)題轉(zhuǎn)化為求待定系數(shù) a0.ak對(duì)等式右邊求 q偏導(dǎo)數(shù)，因而我們得到了：4、把這些等式化簡(jiǎn)并表示成矩陣的形式，就可以得到下面的矩陣：5. 將這個(gè)范德蒙得矩陣化簡(jiǎn)后可得到 ：6. 也就是說(shuō) X*A=Y ，那么 A = (X*X)-1*X*Y ，便得到了系數(shù)矩陣 A ，同時(shí)，我們也就得 到了擬合曲線。MATLAB 實(shí)現(xiàn)：MATLAB 提供了 polyfit ()函數(shù)命令進(jìn)行最小二乘曲

3、線擬合。調(diào)用格式： p=polyfit(x,y,n)p,s= polyfit(x,y,n)p,s,mu=polyfit(x,y,n)x,y 為數(shù)據(jù)點(diǎn)， n 為多項(xiàng)式階數(shù)，返回 p 為冪次從高到低的多項(xiàng)式系數(shù)向量 p。 x 必須是單調(diào)的。矩陣s包括R (對(duì)x進(jìn)行QR分解的三角元素)、df(自由度)、 normr(殘差)用于生成預(yù)測(cè)值的誤差估計(jì)。p,s,mu=polyfit(x,y,n) 在擬合過(guò)程中，首先對(duì) x 進(jìn)行數(shù)據(jù)標(biāo)準(zhǔn)化處理，以在擬合 中消除量綱等影響，mu包含標(biāo)準(zhǔn)化處理過(guò)程中使用的x的均值和標(biāo)準(zhǔn)差。 polyval( ) 為多項(xiàng)式曲線求值函數(shù)，調(diào)用格式 ： y=polyval(p,x)y

4、,DELTA=polyval(p,x,s) y=polyval(p,x)為返回對(duì)應(yīng)自變量x在給定系數(shù)P的多項(xiàng)式的值。y,DELTA=polyval(p,x,s) 使用polyfit函數(shù)的選項(xiàng)輸出s得出誤差估計(jì)丫 DELTA。它假設(shè)polyfit函數(shù)數(shù)據(jù)輸入的誤差是獨(dú)立正態(tài)的，并且方差為常數(shù)。 則 丫 DELTA 將至少包含 50%的預(yù)測(cè)值。如下給定數(shù)據(jù)的擬合曲線： x=0.5,1.0,1.5,2.0,2.5,3.0， y=1.75,2.45,3.81,4.80,7.00,8.60。 解： MATLAB 程序如下： x=0.5,1.0,1.5,2.0,2.5,3.0;y=1.75,2.45,3.

5、81,4.80,7.00,8.60;p=polyfit(x,y,2) x1=0.5:0.05:3.0;y1=polyval(p,x1); plot(x,y,*r,x1,y1,-b) 運(yùn)行結(jié)果如圖 1 計(jì)算結(jié)果為： p =0.5614 0.8287 1.1560 即所得多項(xiàng)式為 y=0.5614xA2+0.08287x+1.15560圖1 最小二乘法曲線擬合示例 對(duì)比檢驗(yàn)擬合的有效性： 例：在0, n區(qū)間上對(duì)正弦函數(shù)進(jìn)行擬合，然后在0,2 n區(qū)間畫(huà)出圖形，比較擬合 區(qū)間和非擬合區(qū)間的圖形，考察擬合的有效性。在 MATLAB 中輸入如下代碼：clearx=0:0.1:pi; y=sin(x);p,

6、mu=polyfit(x,y,9) x1=0:0.1:2*pi;y1=sin(x1);% 實(shí)際曲線y2=polyval(p,x1);%根據(jù)由區(qū)間0到pi上進(jìn)行擬合得到的多項(xiàng)式計(jì)算 0到2pi上 的函數(shù)值，%需要注意的是polyval ()返回的函數(shù)值在pi到2pi上并 沒(méi)有進(jìn)行擬合plot(x1,y2,k*,x1,y1,k-)運(yùn)行結(jié)果：p =0.0000 0.0000 -0.0003 0.0002 0.0080 0.0002 -0.1668 0.0000 1.0000 0.0000mu =R: 10x10 doubledf: 22 normr: 1.6178e-07MATLAB 的最優(yōu)化工具箱

7、還提供了 lsqcurvefit ()函數(shù)命令進(jìn)行最小二乘曲線 擬合(Solve non li near curve-fitti ng (data-fitti ng) problems in least-squares sense)。調(diào)用格式：x = lsqcurvefit(fun,x0,xdata,ydata)x = lsqcurvefit(fun,x0,xdata,ydata,lb,ub)x = lsqcurvefit(fun,x0,xdata,ydata,lb,ub,options)x = lsqcurvefit(problem)x,resnorm = lsqcurvefit(.)x,r

8、esnorm,residual = lsqcurvefit(.)x,resnorm,residual,exitflag = lsqcurvefit(.)x,resnorm,residual,exitflag,output = lsqcurvefit(.) x,resnorm,residual,exitflag,output,lambda = lsqcurvefit(.) x,resnorm,residual,exitflag,output,lambda,jacobian = x0 為初始解向量 ;xdata ， ydata 為滿足關(guān)系 ydata=F(x, xdata) 的數(shù)據(jù) ;lb、ub為

9、解向量的下界和上界，若沒(méi)有指定界，則lb= , ub=;options 為指定的優(yōu)化參數(shù) ;fun 為擬合函數(shù)，其定義方式為 :x = lsqcurvefit(myfun,x0,xdata,ydata) ， 其中 myfun 已定義為 function F = myfun(x,xdata)F =%計(jì)算x處擬合函數(shù)值fun的用法與前面相同；resnorm=sum (fun(x,xdata)-ydataF2)，即在 x 處殘差的平方和； residual=fun(x,xdata)-ydata ，即在 x 處的殘差 ;exitflag 為終止迭代的條件 ;output 為輸出的優(yōu)化信息 ;lambd

10、a 為解 x 處的 Lagrange 乘子 ;jacobian 為解 x 處擬合函數(shù) fun 的 jacobian 矩陣。例： lsqcurvefit() 優(yōu)化程序Data = .0.00005.89550.10003.56390.20002.51730.30001.97900.40001.89900.50001.39380.60001.13590.70001.00960.80001.03430.90000.84351.00000.68561.10000.61001.20000.53921.30000.39461.40000.39031.50000.54741.60000.34591.7000

11、0.13701.80000.22111.90000.17042.00000.2636;t = Data(:,1);y = Data(:,2);% axis(0 2 -0.5 6)plot(t,y,ro)title(Data points)%We would like to fit the function y = c(1)*exp(-lam(1)*t) + c(2)*exp(-lam(2)*t) to the data %The lsqcurvefit function solves this type of problem easily.%To begin, define the param

12、eters in terms of one variable x:%x(1) = c(1)%x(2) = lam(1)%x(3) = c(2)%x(4) = lam(2)%Then define the curve as a function of the parameters x and the data t:F = (x,xdata)x(1)*exp(-x(2)*xdata) + x(3)*exp(-x(4)*xdata);x0 = 1 1 1 0;x,resnorm,exitflag,output = lsqcurvefit(F,x0,t,y)hold onplot(t,F(x,t)ho

13、ld offFsumsquares = (x)sum(F(x,t) - y).A2);opts = optimset(LargeScale,off);xunc,ressquared,eflag,outputu = .fminunc(Fsumsquares,x0,opts)fprintf(There were %d iterations using fminunc, . and %d using lsqcurvefit.n, .outputu.iterations,output.iterations)fprintf(There were %d function evaluations using

14、 fminunc, . and %d using lsqcurvefit., .outputu.funcCount,output.funcCount)type fitvectorx02 = 1 0;F2 = (x,t) fitvector(x,t,y);x2,resnorm2,exitflag2,output2 = lsqcurvefit(F2,x02,t,y)fprintf(There were %d function evaluations using the 2-d .formulation, and %d using the 4-d formulation., . output2.fu

15、ncCount,output.funcCount)x0bad = 5 1 1 0;xbad,resnormbad,exitflagbad,outputbad = .lsqcurvefit(F,x0bad,t,y)hold onplot(t,F(xbad,t),g)legend(Data,Global fit,Bad local fit,Location,NE)hold offfpri ntf(The residual norm at the good ending point is %f, .and the residual norm at the bad ending point is %f., .res no rm,res no rmbad)displayE ndOfDemoMessage(mfile name)擬合效果如下：直線的最小二乘擬合：y= a+bx式中有兩個(gè)待定參數(shù)，a代表截距，b代表斜率。對(duì)于等精度測(cè)量所得到的 N組數(shù)據(jù)(xi，yi)，i = 1, 2,N, xi值被認(rèn)為是準(zhǔn)確的，所有的誤差只聯(lián) 系著yi。下面利用最小二乘法把觀測(cè)數(shù)據(jù)擬合為直線。用最小二乘法估計(jì)參數(shù)時(shí)，要求觀測(cè)值yi的偏差的加權(quán)平方和為最小。對(duì)于等精度觀測(cè)值的直線擬合來(lái)說(shuō)，可使下式的值最?。荷鲜椒謩e對(duì)a、b求偏導(dǎo)得：整理后得到方程組：解上述方程組便可求得直線參數(shù)a