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

1、曲線擬合（curve-fitting）：工程實(shí)踐中，用測(cè)量到的一些離散的數(shù)據(jù)求一個(gè)近似的函數(shù)來(lái)擬合這組數(shù)據(jù)，要求所得的擬合曲線能最好的反映數(shù)據(jù)的基本趨勢(shì)（即使最好地逼近，而不必滿足插值原則。因此沒(méi)必要取=，只要使盡可能地?。Ｔ恚航o定數(shù)據(jù)點(diǎn)。求近似曲線。并且使得近似曲線與的偏差最小。近似曲線在該點(diǎn)處的偏差，i=1,2,.,m。 常見(jiàn)的曲線擬合方法: 1.使偏差絕對(duì)值之和最小 2.使偏差絕對(duì)值最大的最小 3.使偏差平方和最小 最小二乘法： 按偏差平方和最小的原則選取擬合曲線，并且采取二項(xiàng)式方程為擬合曲線的方法,稱(chēng)為最小二乘法。推導(dǎo)過(guò)程： 1. 設(shè)擬合多項(xiàng)式為： 2. 各點(diǎn)到這條曲線的距離之和，

2、即偏差平方和如下： 3. 問(wèn)題轉(zhuǎn)化為求待定系數(shù).對(duì)等式右邊求偏導(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)行最小二乘曲線擬合。調(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ù)

3、向量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,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ì)Y DELTA。它假

4、設(shè)polyfit函數(shù)數(shù)據(jù)輸入的誤差是獨(dú)立正態(tài)的，并且方差為常數(shù)。則Y 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.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 =

5、0.5614 0.8287 1.1560即所得多項(xiàng)式為y=0.5614x2+0.08287x+1.15560 圖1 最小二乘法曲線擬合示例對(duì)比檢驗(yàn)擬合的有效性：例：在0,區(qū)間上對(duì)正弦函數(shù)進(jìn)行擬合，然后在0,2區(qū)間畫(huà)出圖形，比較擬合區(qū)間和非擬合區(qū)間的圖形，考察擬合的有效性。在MATLAB中輸入如下代碼：clearx=0:0.1:pi;y=sin(x);p,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ù)值在

6、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 double df: 22 normr: 1.6178e-07MATLAB的最優(yōu)化工具箱還提供了lsqcurvefit（）函數(shù)命令進(jìn)行最小二乘曲線擬合(Solve nonlinear curve-fitting (data-fitting) problems in least-squares se

7、nse)。調(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,resnorm,residual = lsqcurvefit(.)x,resnorm,residual,exitflag = lsqcurvefit(.)x,resnorm,residual,exitflag,output =

8、 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為解向量的下界和上界 ，若沒(méi)有指定界，則lb= ，ub= ;options為指定的優(yōu)化參數(shù);fun為擬合函數(shù)，其定義方式為:x = lsqcurvefit(myfun,x0,xdata,ydata)，其中myfun已定義為 function F

9、= myfun(x,xdata)F = % 計(jì)算x處擬合函數(shù)值fun的用法與前面相同;resnorm=sum (fun(x,xdata)-ydata).2)，即在x處殘差的平方和;residual=fun(x,xdata)-ydata，即在x處的殘差;exitflag為終止迭代的條件;output為輸出的優(yōu)化信息;lambda為解x處的Lagrange乘子;jacobian為解x處擬合函數(shù)fun的jacobian矩陣。例：lsqcurvefit()優(yōu)化程序Data = . 0.0000 5.8955 0.1000 3.5639 0.2000 2.5173 0.3000 1.9790 0.400

10、0 1.8990 0.5000 1.3938 0.6000 1.1359 0.7000 1.0096 0.8000 1.0343 0.9000 0.8435 1.0000 0.6856 1.1000 0.6100 1.2000 0.5392 1.3000 0.3946 1.4000 0.3903 1.5000 0.5474 1.6000 0.3459 1.7000 0.1370 1.8000 0.2211 1.9000 0.1704 2.0000 0.2636;t = Data(:,1);y = Data(:,2);% axis(0 2 -0.5 6) plot(t,y,'ro'

11、;)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 parameters in terms of one variable x:%x(1) = c(1)%x(2) = lam(1)%x(3) = c(2)%x(4) = lam(2)%Then

12、 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)hold offFsumsquares = (x)sum(F(x,t) - y).2);opts = optimset('LargeScale','off

13、9;);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 fminunc,' . ' and %d using lsqcurve

14、fit.', . 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.funcCou

15、nt,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 offfprintf('The residual norm at the good ending point

16、 is %f,' . ' and the residual norm at the bad ending point is %f.', . resnorm,resnormbad)displayEndOfDemoMessage(mfilename)擬合效果如下：直線的最小二乘擬合：ya+bx式中有兩個(gè)待定參數(shù)，a代表截距，b代表斜率。對(duì)于等精度測(cè)量所得到的N組數(shù)據(jù)（xi，yi），i1，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和b的最