楊 洪，林大偉,張宏禮

(黑龍江八一農(nóng)墾大學(xué) 文理學(xué)院，黑龍江 大慶 163319)

0 引言

在科學(xué)數(shù)據(jù)的處理過程中，常常需要建立變量之間的數(shù)學(xué)模型。所建立的數(shù)學(xué)模型往往是非線性的，因?yàn)樵趯?shí)際應(yīng)用中非線性模型常常比線性模型更能反映變量之間的關(guān)系。數(shù)據(jù)擬合是一種重要的建立變量之間數(shù)學(xué)模型的方法，其求解方法通常是對(duì)給定的含有待定參數(shù)的數(shù)學(xué)模型應(yīng)用最小二乘法求解。隨著研究的不斷深入，非線性數(shù)據(jù)擬合問題越來越受到重視。

1 基于最小二乘法的數(shù)據(jù)擬合

在Matlab下有三個(gè)適于數(shù)據(jù)擬合的命令，它們分別是統(tǒng)計(jì)工具箱下的nlinfit()、優(yōu)化工具箱下的 lsqnonlin()、優(yōu)化工具箱下的lsqcurvefit()，這三個(gè)命令都依據(jù)最小二乘法求解[1]。

最小二乘法一般是求解最小二乘問題，最小二乘問題就是對(duì)由若干個(gè)函數(shù)的平方和構(gòu)成的目標(biāo)函數(shù)求極小值的問題。目標(biāo)函數(shù)一般可以寫成

其中，x=(x1,x2,…xn)是集合En中的點(diǎn)，一般m≥n。目標(biāo)函數(shù)極小化可以表示成

特別地，當(dāng)fi(x),i=1,2,…m中至少有一個(gè)是x的非線性函數(shù)時(shí)，稱為非線性最小二乘問題。

2 三個(gè)命令的算法原理、格式

2.1 nlinfit函數(shù)的算法原理

nlinfit()實(shí)際上是非線性回歸函數(shù)，一般用非線性最小二乘法確定回歸方程中的系數(shù)[3]，其基本算法是Guass-Newton法。

Guass-Newton法可以非常方便地求解最小二乘問題。在Guass-Newton法中需要求解牛頓方程

2f(x(k))d=f(x(k))

由Hesse矩陣

可知，這里需要計(jì)算ri(x)(i=1,2,…m)，計(jì)算量大。為了簡(jiǎn)化計(jì)算，省略Hesse矩陣的第二項(xiàng)，即求解方程組

J(x(k))TJ(x(k))d=-J(x(k))Tr(x(k))

得d(k)，然后置

x(k+1)=x(k)+d(k)

這樣就得到了Guass-Newton法[4]。

2.2 lsqnonlin函數(shù)的算法原理及格式

該函數(shù)主要的算法是Guass-Newton法和Levenberg-Marquardt法。

lsqnonlin函數(shù)求解非線性數(shù)據(jù)擬合問題，不僅僅要計(jì)算目標(biāo)函數(shù)f(x)(平方和)，lsqnonlin函數(shù)還需要用戶指定函數(shù)來計(jì)算向量值函數(shù)

然后，將優(yōu)化問題重新寫成向量的形式：

式中，x為一向量，F(xiàn)(x)為一返回向量值的函數(shù)。

用lsqnonlin函數(shù)求解數(shù)據(jù)擬合問題，其調(diào)用格式為：

x = lsqnonlin(fun,x0)

x = lsqnonlin(fun,x0,lb,ub)

x = lsqnonlin(fun,x0,lb,ub,options)

x = lsqnonlin(fun,x0,lb,ub,options,P1,P2,…)

[x,resnorm] = lsqnonlin(...)

[x,resnorm,residual] = lsqnonlin(...)

[x,resnorm,residual,exitflag] = lsqnonlin(...)

[x,resnorm,residual,exitflag,output] = lsqnonlin(...)

[x,resnorm,residual,exitflag,output,lambda] = lsqnonlin(...)

[x,resnorm,residual,exitflag,output,lambda,jacobian] = lsqnonlin(...)

2.3 lsqcurvefit函數(shù)的算法原理及格式

lsqcurvefit函數(shù)的基本算法與lsqnonlin函數(shù)基本相同。同樣是應(yīng)用了Guass-Newton法和Levenberg-Marquardt法。

lsqcurvefit解決非線性數(shù)據(jù)擬合問題的數(shù)學(xué)模型為

其中xdata和ydata為向量，F(xiàn)(x,xdata)為向量值函數(shù)。

用lsqcurvefit函數(shù)求最小二乘意義上的非線性數(shù)據(jù)擬合問題。即根據(jù)輸入數(shù)據(jù)xdata和得到的輸出數(shù)據(jù)ydata，找到與方程F(x,xdata)最佳的擬合系數(shù)。

lsqcurvefit求解非線性數(shù)據(jù)擬合問題。lsqcurvefit需要一個(gè)用戶定義函數(shù)來計(jì)算向量值函數(shù)F(x,xdata)。用戶定義的函數(shù)的向量的大小必須與ydata的大小相同。

該函數(shù)的調(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,resnorm] = lsqcurvefit(...)

[x,resnorm,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] = lsqcurvefit(...)

3 nlinfit、lsqnonlin和lsqcurvefit命令的實(shí)例比較

nilinfit()函數(shù)是應(yīng)用了Guass-Newton法進(jìn)行求解的，在Guass-Newton法中，當(dāng)J(x(k))的各列是線性相關(guān)，或接近線性相關(guān)，換句話說，矩陣J(x(k))TJ(x(k))奇異或是病態(tài)的，那么求解線性方程組就會(huì)出現(xiàn)一點(diǎn)困難，而lsqnonlin命令是基于Guass-Newton法和Levenberg-Marquardt法進(jìn)行求解。從而利用Levenberg-Marquardt法設(shè)置參數(shù)即可完成。已知數(shù)據(jù)如表1所示

表1 數(shù)據(jù)表

在Matlab中輸入

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];

f=inline('p(1)*x./(p(2)+x)','p','x');

[p,r]=nlinfit(x,y,f,[200,0.1])

運(yùn)行結(jié)果為

p = 212.6826 0.0641

r =25.4332 -3.5668 -5.8119 4.1881 -11.3623 4.6377 -5.6848 -12.6848 0.1676 10.1676 6.0319 -0.9681

下面畫出由擬合得到的數(shù)據(jù)及已知的數(shù)據(jù)散點(diǎn)圖。

x1=0:0.02:1.10; y1=212.6826*x1./(0.0641+x1); plot(x,y, 'o',x1,y1)

圖1利用nlinfit( )數(shù)據(jù)擬合生成圖

在[p,r]=nlinfit(x,y,f,[200,0.1])中，“r”為殘差。殘差是通過建立的擬合函數(shù)求出的值和實(shí)際值之間的差值。

下面應(yīng)用lsqnonlin函數(shù)進(jìn)行擬合，首先編寫目標(biāo)函數(shù) (curve_fun.m)

function y=curve_fun(p)%非線性最小二乘擬合函數(shù)

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()函數(shù)求解，輸入

[p,resnorm,residual]=lsqnonlin(@curve_fun,[200,0.1])

運(yùn)行結(jié)果為

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

作圖如下：

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];

y1=212.6836*x1./(0.0641+x1); plot(x,y,'o',x1,y1)

圖2 利用lsqnonlin ( )數(shù)據(jù)擬合生成圖

雖然lsqnonlin函數(shù)較nlinfit函數(shù)全面，但在一些基本數(shù)據(jù)擬合問題上特別在不是解病態(tài)奇異的方程組的情況下，應(yīng)用nlinfit函數(shù)與lsqnonlin函數(shù)擬合沒有明顯的差別，幾乎可以省略。Matlab用lsqcurvefit函數(shù)進(jìn)行非線性數(shù)據(jù)擬合，其算法與lsqnonlin函數(shù)的算法相同。

4 結(jié)語

基于最小二乘法闡述了MATLAB中三個(gè)適于數(shù)據(jù)擬合命令的算法原理、格式，并通過應(yīng)用實(shí)例對(duì)這三個(gè)命令的異同進(jìn)行了比較分析。nlinfit()函數(shù)實(shí)際上是非線性回歸函數(shù)，其基本原理是Guass-Newton法，它的擬合結(jié)果比lsqurvefit()函數(shù)的默認(rèn)控制結(jié)果更精確，但是因?yàn)槠洳辉试S給出像lsqurvefit()函數(shù)的精度控制選項(xiàng)，在某方面不能進(jìn)一步精確。用lsqcurvefit()函數(shù)進(jìn)行非線性數(shù)據(jù)擬合，設(shè)置該函數(shù)的目的是提供一個(gè)更方便于進(jìn)行數(shù)據(jù)擬合的方法，其算法與lsqnonlin()函數(shù)的算法相同，都是Guass-Newton法和Levenberg-Marquardt法。但在一些特定的環(huán)境不能用lsqnonlin()函數(shù)，比如由于大型算法不能求解待定系統(tǒng)，中型算法不能求解有邊界約束的問題等，這樣的問題應(yīng)用lsqurvefit()函數(shù)更適用更精確些。

[參考文獻(xiàn)]

[1] 羅成漢,劉小山.曲線擬合法的Matlab實(shí)現(xiàn)[J].現(xiàn)代電子技術(shù),2003,20：16-18.

[2] 薛毅.最優(yōu)化理論與算法[M].北京:北京工業(yè)大學(xué)出版社,2004：201-202.

[3] 包翠蓮,開小明.Matlab語言在多元線性回歸中的應(yīng)用[J]. 安徽教育學(xué)院學(xué)報(bào),2005,23(3):55.

[4] 禇洪生,杜增吉,閻金華.Matlab7.2優(yōu)化設(shè)計(jì)實(shí)例指導(dǎo)教程[M].北京:機(jī)械工業(yè)出版社,2007：117.