最小二乘法原理
函數插值是內插補點函數p(x)與被插函數f(x)在節點處函數值相同,即p( )=f( ) (i=0,1,2,3……,n),而曲線拟合函數 不要求嚴格地通過所有資料點( ),也就是說拟合函數 在 處的偏差
=
不都嚴格地等于零。但是,為了使近似曲線能盡量反應所給資料點的變化趨勢,要求| |按某種度量标準最小。
即
=
為最小。這種要求誤差平方和最小的拟合稱為曲線拟合的最小二乘法。
(一)線性最小二乘拟合
根據線性最小二乘拟合理論,我們得知關于系數矩陣A的解法為A=R\Y。
例題 假設測出了一組,由下面的表格給出,且已知函數原型為
y(x)=c1+c2*e^(-3*x)+c3*cos(-2*x)*exp(-4*x)+c4*x^2
x | 0.2 | 0.4 | 0.7 | 0.9 | 0.92 | 0.99 | 1.2 | 1.4 | 1.48 | 1.5 | |
y | 2.88 | 2.2576 | 1.9683 | 1.9258 | 2.0862 | 2.109 | 2.1979 | 2.5409 | 2.9627 | 3.155 | 3.2052 |
試用已知資料求出待定系數的值。
在Matlab中輸入以下程式
x=[0,0.2,0.4,0.7,0.9,0.92,0.99,1.2,1.4,1.48,1.5]';
y=[2.88;2.2576;1.9683;1.9258;2.0862;2.109;
2.1979;2.5409;2.9627;3.155;3.2052];
A=[ones(size(x)) exp(-3*x),cos(-2*x).*exp(-4*x) x.^2];
c=A\y;
c'
運作結果為
ans =
1.2200 2.3397 -0.6797 0.8700
下面畫出由拟合得到的曲線及已知的資料散點圖
x1=[0:0.01:1.5]';
A1=[ones(size(x1)) exp(-3*x1),cos(-2*x1).*exp(-4*x1) x1.^2];
y1=A1*c;
plot(x1,y1,x,y,'o')
事實上,上面給出的資料就是由已知曲線
y(x)= 0.8700-0.6797*e^(-3*x)+ 2.3397*cos(-2*x)*exp(-4*x)+ 1.2200*x^2
産生的,由上圖可見拟合效果較好。
多項式最小二乘拟合
在Matlab的線性最小二乘拟合中,用得較多的是多項式拟合,其指令是
A=polyfit(x,y,m)
其中 表示函數中的自變量矩陣, 表示因變量矩陣,是輸出的系數矩陣,即多項式的系數。
多項式在自變量x處的函數值y可用以下指令計算:
y=polyval(A,x)
例題 對下面一組資料作二次多項式拟合,即要求出二次多項式 中的 ,使最小。
x | 0.1 | 0.2 | 0.3 | 0.4 | 0.5 | 0.6 | 0.7 | 0.8 | 0.9 | 1 | |
y | -0.447 | 1.978 | 3.28 | 6.16 | 7.08 | 7.34 | 7.66 | 9.56 | 9.48 | 9.30 | 11.2 |
在Matlab中輸入以下指令
x=0:.1:1;
y=[-0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2];
a=polyfit(x,y,2)
運作結果為
a =
-9.8108 20.1293 -0.0317
f=vpa(poly2sym(a),5)
%vpa(polyval2sym(),n)隻适用于關于多項式函數的拟合。因為此函數對于自變量統一規定為“x”,将由polyfit()所得出的系數按自變量幂次升降放在相應的位置。
運作結果為
f =
-9.8108*x^2+20.129*x-.31671e-1
下面畫出由拟合得到的曲線及已知的資料散點圖
y1=polyval(a,x);
plot(x,y,'o',x,y1)
(二)非線性最小二乘拟合
(1)lsqcurvefit( )
lsqcurvefit( )是非線性最小二乘拟合函數,其本質上是求解
最優化問題。其使用格式為
x=lsqcurvefit(fun,x0,xdata,ydata)
其中,fun是要拟合的非線性函數,x0是初始參數,xdata,ydata是拟合點的資料,該函數最終傳回系數矩陣。
例題 假設已知
并已知該函數滿足原型為 ,其中 為待定系數。
在Matlab中輸入以下指令
x=0:.1:10;
y=0.12*exp(-0.213*x)+0.54*exp(-0.17*x).*sin(1.23*x);
f=inline('a(1)*exp(-a(2)*x)+a(3)*exp(-a(4)*x).*sin(a(5)*x)','a','x');
%建立函數原型,則可以根據他來進行下面的求取系數的計算
[a,res]=lsqcurvefit(f,[1,1,1,1,1],x,y);
a',res
運作結果為
ans =
0.1197
0.2125
0.5404
0.1702
1.2300
res =
7.1637e-007
所求得的系數與原式中的系數相近。
如果向進一步提高精度,則需修改最優化的選項,函數的調用格式也将随之改變。
在Matlab中輸入以下指令
ff=optimset;ff.TolFun=1e-20;ff.TolX=1e-15;%修改精度,即是修改其限制條件
[a,res]=lsqcurvefit(f,[1,1,1,1,1],x,y,[],[],ff);%兩個空矩陣表示系數向量的上下限
a',res
運作結果為
ans =
0.1200
0.2130
0.5400
0.1700
1.2300
res =
9.5035e-021
下面繪圖
x1=0:0.01:10;
y1=f(a,x1);
plot(x1,y1,x,y,'o')
(2)lsqnonlin( )
lsqnonlin( )函數是另一種求最小二乘拟合的函數,其本質上是求解最優化問題
最優化解。它的應用格式為
x=lsqnonlin(fun,x0)
其中,fun的定義與lsqnonlin( )函數中fun的定義有差别, x0仍為初始參數向量,将輸出的系數結果放在變量 中。
說明lsqnonlin()函數的使用方法。
首先編寫目标函數 (curve_fun.m)
function y=curve_fun(p)%非線性最小二乘拟合函數
x=[0.02 0.02 0.06 0.06 0.11 0.11 0.22 0.22 0.56 0.56 1.10 1.10];
y=[76 47 97 107 123 139 159 152 191 201 207 200];
y=p(1)*x./(p(2)+x)-y;
再用lsqnonlin()函數求解,輸入
[p,resnorm,residual]=lsqnonlin(@curve_fun,[200,0.1])
運作結果為
p =
212.6836 0.0641
resnorm =
1.1954e+003
residual =
Columns 1 through 11
-25.4339 3.5661 5.8111 -4.1889 11.3617 -4.6383 5.6847 12.6847 -0.1671 -10.1671 -6.0313
Column 12
0.9687
上面的兩種方法都可以作非線性最小二乘曲線拟合
(3)非線性函數的線性化
在進行非線性拟合時,以往由于計算機和相關軟體水準有限,常常先把非線性的曲線拟合線性化,然後再進行拟合。下面比較一下先線性化再進行拟合和直接進行非線性拟合的差異。
例題 已知資料
t | 0.25 | 0.5 | 1 | 1.5 | 2 | 3 | 4 | 6 | 8 |
c | 19.21 | 18.15 | 15.36 | 14.10 | 12.89 | 9.32 | 7.45 | 5.24 | 3.01 |
滿足曲線 通過資料拟合求出參數和 。
方法一:先将非線性函數轉化為線性函數
編寫Matlab程式如下
d=300;
t=[0.25 0.5 1 1.5 2 3 4 6 8];
c=[19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01];
y=log(c);
a=polyfit(t,y,1)
運作結果為
a =
-0.2347 2.9943
k=-a(1)
k =
0.2347
v=d/exp(a(2))
v =
15.0219
由此也可以求出相關系數。
方法二:應用非線性拟合直接求解系數
建立m檔案:
function f=curvefun3(x,tdata)
d=300
f=(x(1)\d)*exp(-x(2)*tdata)%x(1)=v;x(2)=k
運作程式
tdata=[0.25 0.5 1 1.5 2 3 4 6 8];
cdata=[19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01];
x0=[10 0.5];
x=lsqcurvefit('curvefun3',x0,tdata,cdata)
運作結果為
x =
14.8212 0.2420
下面繪圖
f=curvefun3(x,tdata);
plot(tdata,cdata,'o',tdata,f)
我們發現兩種求法求出的系數很接近。
(三)線性拟合和非線性拟合差別與聯系
在許多實際問題中,變量之間内在的關系并不想前面說的那樣簡單。呈線性關系,但有些非線性拟合曲線可以通過适當的變量替換轉化為線性曲線,進而用線性拟合進行處理。對于一個實際的曲線拟合問題,一般先根據觀測值在直角坐标平面上描出散點圖,看一看散點的分布同哪類曲線圖形接近,讓後選用相接近的曲線拟合方程,再通過适當的變量替換轉化為線性拟合問題,按線性拟合解出後再還原為原變量所表示的曲線拟合方程。
表1.1
線性拟合方程 | 變量變換 | 變換後線性拟合方程 |
Y= | , | |
Y=a | Y=a | |
Y= | , | |
Y= | ||
Y= |
例題
測出一組實際資料 見下表 是對其進行函數拟合。
X | 1.1052 | 1.2214 | 1.3499 | 1.4918 | 1.6478 | 3.6693 |
Y | 0.6795 | 0.6006 | 0.5309 | 0.4693 | 0.4148 | 0.1546 |
X | 1.8221 | 2.0138 | 2.2255 | 2.4596 | 2.7183 | |
Y | 0.3666 | 0.3241 | 0.2864 | 0.2532 | 0.2238 |
>>x=[1.1052,1.2214,1.3499,1.4918,1.6478,1.8221,2.0138,2.2255,2.4596,2.7183,3.6693];
>>y=[0.6795,0.6006,0.5309,0.4693,0.4148,0.3666,0.3241,0.2864,0.2532,0.2238,0.1546];
>>plot(x,y,x,y,'*')
見下圖
由上圖可以看出是非線性曲線y=a 由表1.1第一行可知
Y= | , |
分别對x,y進行對數變換
>>x1=log(x);y1=log(y);plot(x1,y1)
用線性函數拟合的方法可以求出現行參數,使得ln y=aln x+b,即y= ,求解系數a,b及 。
>> A=[x1',ones(size(x1'))];c=[A\y1']
c =
-1.2338
-0.2631
>> exp(c(2))
ans =
0.7686
拟合函數為y(x)=0.768 703 388 199 24