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

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

2、即偏差平方和如下： 3. 問題轉(zhuǎn)化為求待定系數(shù).對等式右邊求偏導(dǎo)數(shù)，因而我們得到了： . 4、 把這些等式化簡并表示成矩陣的形式，就可以得到下面的矩陣： 5. 將這個范德蒙得矩陣化簡后可得到: 6. 也就是說X*A=Y，那么A = (X*X)-1*X*Y，便得到了系數(shù)矩陣A，同時，我們也就得到了擬合曲線。MATLAB實現(xiàn)：MATLAB提供了polyfit（）函數(shù)命令進行最小二乘曲線擬合。調(diào)用格式：p=polyfit(x,y,n) p,s= polyfit(x,y,n) p,s,mu=polyfit(x,y,n)x,y為數(shù)據(jù)點，n為多項式階數(shù)，返回p為冪次從高到低的多項式系數(shù)向量p。x必須是單調(diào)

3、的。矩陣s包括R（對x進行QR分解的三角元素）、df(自由度）、normr(殘差）用于生成預(yù)測值的誤差估計。 p,s,mu=polyfit(x,y,n)在擬合過程中，首先對x進行數(shù)據(jù)標準化處理，以在擬合中消除量綱等影響，mu包含標準化處理過程中使用的x的均值和標準差。polyval( )為多項式曲線求值函數(shù)，調(diào)用格式： y=polyval(p,x) y,DELTA=polyval(p,x,s)y=polyval(p,x)為返回對應(yīng)自變量x在給定系數(shù)P的多項式的值。y,DELTA=polyval(p,x,s) 使用polyfit函數(shù)的選項輸出s得出誤差估計Y DELTA。它假設(shè)polyfit函數(shù)

4、數(shù)據(jù)輸入的誤差是獨立正態(tài)的，并且方差為常數(shù)。則Y DELTA將至少包含50%的預(yù)測值。如下給定數(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)運行結(jié)果如圖1計算結(jié)果為：p =0.5614 0.8287 1.1560即所得多項式為y=0

5、.5614x2+0.08287x+1.15560 圖1 最小二乘法曲線擬合示例對比檢驗擬合的有效性：例：在0,區(qū)間上對正弦函數(shù)進行擬合，然后在0,2區(qū)間畫出圖形，比較擬合區(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);%實際曲線y2=polyval(p,x1);%根據(jù)由區(qū)間0到pi上進行擬合得到的多項式計算0到2pi上的函數(shù)值， %需要注意的是polyval（）返回的函數(shù)值在pi到2pi上并沒有進行擬合plot(x1,y2,k*,x1

6、,y1,k-)運行結(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ù)命令進行最小二乘曲線擬合(Solve nonlinear curve-fitting (data-fitting) problems in least-squares sense)。調(diào)用格式：x = lsqcurvefit(fun,x0,xdata,ydata)x = l

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

8、lambda = lsqcurvefit(.)x,resnorm,residual,exitflag,output,lambda,jacobian = x0為初始解向量;xdata，ydata為滿足關(guān)系ydata=F(x, xdata)的數(shù)據(jù);lb、ub為解向量的下界和上界 ，若沒有指定界，則lb= ，ub= ;options為指定的優(yōu)化參數(shù);fun為擬合函數(shù)，其定義方式為:x = lsqcurvefit(myfun,x0,xdata,ydata)，其中myfun已定義為 function F = myfun(x,xdata)F = % 計算x處擬合函數(shù)值fun的用法與前面相同;resnorm

9、=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.4000 1.8990 0.5000 1.3938 0.6000 1.1359 0.7000 1.0096

10、 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)title(Data points)%We would like to fit the function y = c

11、(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 define the curve as a function of the parameters x and the data t:F =

12、 (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);xunc,ressquared,eflag,outputu = .fminunc(Fsumsquares,x0,opts)fprintf(There were %d ite

13、rations using fminunc, . and %d using lsqcurvefit.n, . outputu.iterations,output.iterations)fprintf(There were %d function evaluations using fminunc, . and %d using lsqcurvefit., . outputu.funcCount,output.funcCount)type fitvectorx02 = 1 0;F2 = (x,t) fitvector(x,t,y);x2,resnorm2,exitflag2,output2 =

14、lsqcurvefit(F2,x02,t,y)fprintf(There were %d function evaluations using the 2-d . formulation, and %d using the 4-d formulation., . output2.funcCount,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,Ba

15、d local fit,Location,NE)hold offfprintf(The residual norm at the good ending point is %f, . and the residual norm at the bad ending point is %f., . resnorm,resnormbad)displayEndOfDemoMessage(mfilename)擬合效果如下：直線的最小二乘擬合：ya+bx式中有兩個待定參數(shù)，a代表截距，b代表斜率。對于等精度測量所得到的N組數(shù)據(jù)（xi，yi），i1，2，N，xi值被認為是準確的，所有的誤差只聯(lián)系著yi。下面

16、利用最小二乘法把觀測數(shù)據(jù)擬合為直線。 用最小二乘法估計參數(shù)時，要求觀測值yi的偏差的加權(quán)平方和為最小。對于等精度觀測值的直線擬合來說，可使下式的值最?。荷鲜椒謩e對a、b求偏導(dǎo)得：整理后得到方程組：解上述方程組便可求得直線參數(shù)a和b的最佳估計值。1、 可看成是一階多項式擬合，跟前面曲線擬合方法一樣。2、 利用linefit()函數(shù)進行最小二乘的直線擬合 使用： clear x=0.5 1 1.5 2 2.5 3; y= 1.75 2.45 3.81 4.8 8 8.6; k,b=linefit(x,y) %得到斜率k和常數(shù)b y1=polyval(k,b,x); plot(x,y1,k-,x,y,k*)MATLAB一元到多元線性回歸方程的計算和