最小二乘法原理

函式插值是差值函式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

試用已知資料求出待定係數的值。

在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)

例題 對下面一組資料作二次多項式擬合，即要求出二次多項式 中的 ，使最小。

在Matlab中輸入以下命令

x=0:.1:1;

y=[-0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.489.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.560.56 1.10 1.10];

y=[76 47 97 107 123 139 159 152 191 201 207200];

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）非線性函式的線性化

在進行非線性擬合時，以往由於計算機和相關軟體水平有限，常常先把非線性的曲線擬合線性化，然後再進行擬合。下面比較一下先線性化再進行擬合和直接進行非線性擬合的差異。

例題 已知資料

滿足曲線 通過資料擬合求出引數和 。

方法一：先將非線性函式轉化為線性函式

編寫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.243.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.243.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

例題

測出一組實際資料 見下表 是對其進行函式擬合。

>>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第一行可知

分別對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