第八章插值 擬合與回歸

1、第八章 插值 擬合與回歸在生產(chǎn)和實(shí)驗(yàn)中，由于函數(shù)的表達(dá)式不便于計算或者沒有表達(dá)式而只有在給定點(diǎn)的函數(shù)值（或其導(dǎo)數(shù)值），為此，我們希望建立一個簡單的而且便于計算的近似函數(shù)，來逼近函數(shù)，這就用到插值和擬合方法。8.1 插值插值就是在離散數(shù)據(jù)的基礎(chǔ)上補(bǔ)插連續(xù)函數(shù)，使得這條連續(xù)曲線通過全部給定的離散數(shù)據(jù)點(diǎn)。插值是離散函數(shù)逼近的重要方法，利用它可通過函數(shù)在有限個點(diǎn)處的取值狀況，估算出函數(shù)在其他點(diǎn)處的近似值。MATLAB中的插值函數(shù)主要有以下幾個。8.1.1 interp1函數(shù)MATLAB中用于一維數(shù)據(jù)插值的函數(shù)是interp1，其調(diào)用格式為：yi = interp1(x, y, xi, 'met

2、hod, 'extrap')該命令用于找出由參量x決定的一元函數(shù)y=y(x)在點(diǎn)xi處的值yi。其中x,y為插值節(jié)點(diǎn)的橫坐標(biāo)和縱坐標(biāo)，yi為在被插值點(diǎn)xi處的插值結(jié)果；x,y為向量， method表示采用的插值方法，MATLAB提供的插值方法有幾種：nearest：最近鄰點(diǎn)插值算法；linear：線性插值（默認(rèn)）；spline：三次樣條函數(shù)插值；pchip：分段三次Hermite插值；cubic：與pchip操作相同；v5cubic：在MATLAB 5.0中的三次插值。缺省時表示線性插值。 extrap表示對于超出x 范圍的xi 中的分量將執(zhí)行特殊的外插值法extrap。注意：

3、所有的插值方法都要求x是單調(diào)的。例1 對函數(shù)從0到10，步長為1進(jìn)行采樣得到插值節(jié)點(diǎn)，利用線性插值給出步長為0.25的插值函數(shù)圖形。在命令窗口輸入：>> x = 0:10;y = x.*sin(x); xx = 0:.25:10; yy = interp1(x,y,xx); plot(x,y,'kd',xx,yy)輸出結(jié)果如圖8.1所示。例2 在一天24小時內(nèi)，從零點(diǎn)開始每間隔2小時測得的環(huán)境溫度數(shù)據(jù)分別為12，9，9，1，0，18 ，24，28，27，25，20，18，15，13推測中午13時的溫度。我們采用三次樣條插值，在命令窗口輸入：>> x=0:

4、2:24;a=13;y=12 9 9 10 18 24 28 27 25 20 18 15 13;y1=interp1(x,y,a,'spline')輸出結(jié)果為：y1 = 27.8725若要得到一天24小時的溫度曲線，則輸入：xi=0:1/3600:24;yi=interp1(x,y,xi, 'spline');plot(x,y,'o' ,xi,yi)溫度函數(shù)曲線圖如圖8.2所示。 圖8.1 線性插值函數(shù)曲線圖 圖8.2 溫度曲線圖 8.1.2 interp2函數(shù)MATLAB中用于二維數(shù)據(jù)插值的函數(shù)是interp2，其調(diào)用格式為：zi = int

5、erp2(x, y, z, xi, yi, method)該命令用于找出由參量x,y決定的二元函數(shù)z=z(x,y)在點(diǎn)（xi,yi）處的值zi。其中返回矩陣為zi，其元素為對應(yīng)于參量xi與yi（可以是向量、或同型矩陣）的元素，若xi與yi 中有在x與y范圍之外的點(diǎn)，則相應(yīng)地返回NaN（Not a Number），method和interp1一樣，常用的是cubic （雙三次插值），缺省為linear（雙線性插值算法）。例3 在命令窗口輸入X,Y = meshgrid(-3:.25:3);Z = peaks(X,Y); %具有兩個變量的采樣函數(shù)，可產(chǎn)生一個凹凸有致的曲面，包含了三個局部極大點(diǎn)及三

6、個局部極小點(diǎn)。XI,YI = meshgrid(-3:.125:3);ZZ = interp2(X,Y,Z,XI,YI);surf(X,Y,Z);hold on;surf(XI,YI,ZZ+15) %為作比較，將插值曲面向上平移15單位。hold off輸出結(jié)果如圖8.3圖8.3 二維插值曲面8.1.3 interp3函數(shù)MATLAB中用于三維數(shù)據(jù)插值的函數(shù)是interp3，其調(diào)用格式為：vi = interp3(x, y, z, v, xi, yi, zi, method)該命令用于找出由參量x,y,z決定的三元函數(shù)v=v(x,y,z)在點(diǎn)（xi,yi,zi）處的值vi。參量xi,yi,zi

7、是同型陣列或向量。若向量參量xi,yi,zi是不同長度，不同方向（行或列）的向量，這時輸出參量vi與y1,y2,y3為同型矩陣。其中y1,y2,y3為用命令meshgrid(XI,YI,ZI)生成的同型陣列。若插值點(diǎn)(xi,yi,zi)中有位于點(diǎn)(x,y,z)之外的點(diǎn)，則相應(yīng)地返回特殊變量值NaN。method和interp1一樣。例4 三維插值舉例命令窗口輸入 x,y,z,v = flow(20);xx,yy,zz = meshgrid(.1:.25:10, -3:.25:3, -3:.25:3);vv = interp3(x,y,z,v,xx,yy,zz);slice(xx,yy,zz,v

8、v,6 9.5,1 2,-2 .2); shading interp;colormap cool輸出圖形如圖8.4所示。 圖8.4 三維插值結(jié)果8.1.4 griddata函數(shù) griddata也是一種常用的二維插值方法，其調(diào)用格式為：zi =griddata (x, y, z, xi, yi, method)該命令用于找出由參量x,y決定的二元函數(shù)z=z(x,y)在點(diǎn)（xi,yi）處的值zi。它和interp2的區(qū)別在于，interp2的插值數(shù)據(jù)必須是矩形域，即已知數(shù)據(jù)點(diǎn)（x,y）組成規(guī)則的矩陣，可使用meshgid生成。而griddata函數(shù)的已知數(shù)據(jù)點(diǎn)（x,y）不要求規(guī)則排列，特別是對試

9、驗(yàn)中隨機(jī)沒有規(guī)律采取的數(shù)據(jù)進(jìn)行插值具有很好的效果。method包括：linear（線性插值，為缺省算法）；natural（自然鄰點(diǎn)插值），nearest（最近鄰點(diǎn)插值）； cubic（ 基于三角形的三次插值）和v4（Matlab4中的griddata 算法）。例5 在命令窗口輸入：rand('seed',0)x = rand(100,1)*4-2; y = rand(100,1)*4-2;z = x.*exp(-x.2-y.2);ti = -2:.125:2; XI,YI = meshgrid(ti,ti);ZI = griddata(x,y,z,XI,YI);mesh(XI,

10、YI,ZI), hold;plot3(x,y,z,'o'), hold off;輸出結(jié)果如圖8.5所示。圖8.5 griddata插值結(jié)果例6 有一組散亂數(shù)據(jù)點(diǎn)矩陣如下：A=1.486,3.059,0.1;2.121,4.041,0.1;2.570,3.959,0.1;3.439,4.396,0.1;4.505,3.012,0.1;3.402,1.604,0.1;2.570,2.065,0.1;2.150,1.970,0.1;1.794,3.059,0.2;2.121,3.615,0.2;2.570,3.473,0.2;3.421,4.160,0.2;4.271,3.036,0

11、.2;3.411,1.876,0.2;2.561,2.562,0.2;2.179,2.420,0.2;2.757,3.024,0.3;3.439,3.970,0.3;4.084,3.036,0.3;3.402,2.077,0.3;2.879,3.036,0.4;3.421,3.793,0.4;3.953,3.036,0.4;3.402,2.219,0.4;3.000,3.047,0.5;3.430,3.639,0.5;3.822,3.012,0.5;3.411,2.385,0.5;3.103,3.012,0.6;3.430,3.462,0.6;3.710,3.036,0.6;3.402,2.5

12、62,0.6;3.224,3.047,0.7;3.411,3.260,0.7;3.542,3.024,0.7;3.393,2.763,0.7;x=A(:,1);y=A(:,2);z=A(:,3);X,Y,Z=griddata(x,y,z,linspace(min(x),max(x),20)',linspace(min(y),max(y),20),'v4');mesh(X,Y,Z); hold onplot3(x,y,z,'o'); hold off;輸出結(jié)果如圖8.6所示。圖8.6 散亂數(shù)據(jù)插值圖8.1.5 spline函數(shù)該函數(shù)是利用三次樣條對數(shù)據(jù)進(jìn)行

13、插值，其調(diào)用格式為：yy = spline(x,y,xx)該命令用三次樣條插值計算出由向量x與y 確定的一元函數(shù)y=f(x)在點(diǎn)xx 處的值。若參量y 是一矩陣，則以y 的每一列和x 配對，再分別計算由它們確定的函數(shù)在點(diǎn)xx 處的值。則yy 是一個階數(shù)為length(xx)*size(y,2)的矩陣。例7 對離散地分布在y=exp(x)sin(x)函數(shù)曲線上的數(shù)據(jù)點(diǎn)進(jìn)行樣條插值。在命令窗口輸入：x = 0 2 4 5 8 12 12.8 17.2 19.9 20;y = exp(x).*sin(x);xx = 0:.25:20;yy = spline(x,y,xx);plot(x,y,'

14、;o',xx,yy) 輸出結(jié)果如圖8.7。圖8.7 三次樣條插值8.2 擬合擬合就是用連續(xù)曲線近似地刻畫或比擬平面上離散點(diǎn)組所表示的坐標(biāo)之間的函數(shù)關(guān)系的一種數(shù)據(jù)處理方法。如果已知某函數(shù)的若干離散函數(shù)值f1,f2,fn，通過調(diào)整該函數(shù)中若干待定系數(shù)f(1, 2,n)，使得該函數(shù)與已知點(diǎn)集的差別（最小二乘意義）最小。如果待定函數(shù)是線性，就叫線性擬合或者線性回歸(主要在統(tǒng)計中)，否則叫做非線性擬合或者非線性回歸。MATLAB中提供了線性最小二乘擬合和非線性最小二乘擬合函數(shù)。8.2.1 polyfit函數(shù)ployfit函數(shù)是多項(xiàng)式擬合函數(shù)，其調(diào)用格式為：p = polyfit(x,y,n)其中

15、x,y為長度相同的向量，n為擬合多項(xiàng)式的次數(shù)，返回值p是擬合多項(xiàng)式的系數(shù)向量，冪次由高到低，擬合多項(xiàng)式在x處的值可以通過y=polyval(p,x)來計算。例8，2004年全國大學(xué)生數(shù)學(xué)建模競賽C題，飲酒駕車問題，時間和血液中酒精濃度的函數(shù)關(guān)系，可以利用polyfit進(jìn)行擬合。可輸入：>> time=0.25, 0.5 0.75, 1 1.5 2 2.5 3 3.5 4 4.5 5 6 7 8 9 10 11 12 13 14 15 16'vol=20 68 75 82 82 77 68 68 58 51 50 41 38 35 28 25 18 15 12 10 7 7

16、4'plot(time,vol,'o')p=polyfit(time,vol,7) % 7次多項(xiàng)式擬合f=polyval(p,time); %求多項(xiàng)式的值hold on plot(time,f)hold off輸出7多項(xiàng)次的系數(shù)向量為：p = 0.0002 -0.0113 0.2817 -3.6307 25.5345 -94.6710 153.8757 -1.2592擬合曲線如圖8.8所示。圖8.8 血液酒精濃度擬合曲線MATLAB提供了兩個非線性最小二乘擬合函數(shù)：lsqcurvefit和lsqnonlin。兩個命令都要先建立M-文件fun.m，在其中定義函數(shù)f(x)，

17、但兩者定義f(x)的方式是不同的。8.2.2 lsqcurvefit函數(shù)該函數(shù)用來進(jìn)行非線性擬合，其調(diào)用格式為：x = lsqcurvefit ('fun',x0,xdata,ydata,options);其中，fun為事先建立的擬合函數(shù)F(x,xdata)，其中自變量x表示擬合函數(shù)中的待定參數(shù)，xdata為已知擬合節(jié)點(diǎn)的x坐標(biāo)，x0為待定參數(shù)x的迭代初始值，xdata,ydata為已知數(shù)據(jù)點(diǎn)的x和y坐標(biāo)，options是一些控制參數(shù)。lsqcurvefit函數(shù)用來求含參數(shù)x（向量）的向量值函數(shù)F(x,xdata)=f(x,xdata1), f(x,xdata2), f(x,x

18、datan)中的參數(shù)x（向量），使得最小。例9 根據(jù)表1中的數(shù)據(jù)，利用lsqcurvefit函數(shù)擬合。表1 已知數(shù)據(jù)點(diǎn)X1002003004005006007008009001000y×1034.544.995.355.655.906.106.266.396.506.59首先建立擬合函數(shù)的M文件myfit1.m，其內(nèi)容如下：function f = myfit1(x,xdata) f = x(1)+x(2)*exp(-0.02*x(3)*xdata); 其中x(1),x(2),x(3)分別表示擬合曲線中的參數(shù)a,b,k。然后在命令窗口輸入：>> xdata = 100:1

19、00:1000;ydata = 1e-03*4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59;x0=0.2,0.05,0.05;x = lsqcurvefit ('myfit1',x0,xdata,ydata)f= myfit1(x,xdata) 輸出結(jié)果為：x = 0.0069 -0.0029 0.0809f = Columns 1 through 9 0.0044 0.0048 0.0051 0.0054 0.0056 0.0058 0.0060 0.0061 0.0063 Column 10 0.0064 %擬合函數(shù)在數(shù)據(jù)

20、點(diǎn)的函數(shù)值擬合曲線如圖8.9。>> plot(xdata,ydata,'o',xdata,f)圖8.9 擬合曲線8.2.3 lsqnonlin函數(shù)該函數(shù)用來進(jìn)行非線性擬合，其調(diào)用格式為：x = lsqnonlin ('fun',x0, options);其中fun為事先建立的擬合函數(shù)F(x)，其中自變量x表示擬合函數(shù)中的待定參數(shù)，x0為待定參數(shù)x的迭代初始值，options是一些控制參數(shù)。由于lsqnonlin中定義的擬合函數(shù)的自變量是x，所以已知參數(shù)xdata,ydata應(yīng)寫在該函數(shù)中。lsqnonlin函數(shù)用來求含參量x（向量）的向量值函數(shù)f(x

21、)=f1(x), f2(x),，fn(x)中的參量x，使得最小。其中fi(x)=f(x,xdatai,ydaytai)=F(x, xdatai)- ydatai。例10 根據(jù)表1中的數(shù)據(jù)，利用lsqnonlin函數(shù)擬合。首先建立擬合函數(shù)的M文件myfit2.m，其內(nèi)容如下：function f=myfit2(x)xdata=100:100:1000;ydata=1e-03*4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59;f=x(1)+x(2)*exp(-0.02*x(3)*xdata)- ydata;然后在命令窗口輸入：>> x0

22、=0.2,0.05,0.05;x = lsqnonlin('myfit2',x0)f = myfit2(x)輸出結(jié)果為：x = 0.0069 -0.0029 0.0809f = 1.0e-03 * Columns 1 through 9 -0.0971 -0.1739 -0.2164 -0.2463 -0.2666 -0.2713 -0.2651 -0.2538 -0.2436 Column 10 -0.2313該結(jié)果與lsqcurvefit函數(shù)擬合的結(jié)果相同。8.3 回歸分析MATLAB中提供了一些線性和非線性回歸分析函數(shù)。8.3.1 regress函數(shù)一般地，稱由確定的模型

23、為一元線性回歸模型，記為固定的未知參數(shù)、稱為回歸系數(shù)，自變量x也稱為回歸變量。一般地，稱為高斯馬爾柯夫線性模型(k元線性回歸模型)。，MATLAB提供了多元線性回歸函數(shù)regress，采用的是最小二乘估計，其調(diào)用格式有：b = regress (y,x)返回值為線性模型y = x*b的回歸系數(shù)向量。其中x為n×(k+1) 矩陣，行對應(yīng)于觀測值，列對應(yīng)于預(yù)測變量，y為n維向量，為因變量，一元線性回歸可取k=1。b,bint,r,rint,stats = regress(y,x,alpha)其中bint是回歸系數(shù)的區(qū)間估計，r是殘差，rint是置信區(qū)間，stats是用于檢驗(yàn)回歸模型的統(tǒng)計

24、量，有四個數(shù)值：相關(guān)系數(shù)r2，F(xiàn)值，與F對應(yīng)的概率P和誤差方差估計，alpha是顯著性水平(缺省的時候?yàn)?.05)。相關(guān)系數(shù)r2越接近1，說明回歸方程越顯著；與F對應(yīng)的概率P<alpha時候拒絕H0，回歸模型成立。例11輸入>> y=8.8818 8.9487 9.0541 9.1545 9.2693 9.4289 9.6160 9.8150 9.9825 10.1558 10.3193'k= 7.8381 7.9167 8.0048 8.1026 8.2556 8.5822 8.8287 9.0756 9.2175 9.4148 9.6198'l= 8.38

25、71 8.3872 8.3935 8.3971 8.4025 8.4048 8.4079 8.4141 8.4261 8.4377 8.4444't= 9.9551 9.9057 10.0972 9.9537 9.9370 9.9449 9.9636 10.1291 10.1573 10.2944 10.2093'x= ones(size(k) k l t; b,bint,r,rint,stats=regress(y,x)輸出結(jié)果為b = -55.4988 0.5644 7.1254 0.0222bint = -83.6008 -27.3968 0.4693 0.6594 3.

26、5740 10.6769 -0.1914 0.2359r = -0.0262 -0.0032 0.0033 0.0261 0.0164 -0.0249 0.0006 0.0124 0.0137 -0.0101 -0.0081rint = -0.0620 0.0097 -0.0479 0.0415 -0.0306 0.0372 -0.0111 0.0632 -0.0218 0.0546 -0.0621 0.0124 -0.0385 0.0397 -0.0190 0.0438 -0.0304 0.0577 -0.0460 0.0259 -0.0420 0.0258stats = 1.0e+003

27、* 0.0010 2.1022 0.0000 0.0000該結(jié)果說明y=-55.4988+0.5644k+7.1254l+0.0222t，stats中的數(shù)據(jù)說明r2=1，F(xiàn)=2102.2，p=0，由于p<0.05可知回歸模型成立。可以利用rcoplot函數(shù)畫出殘差及其置信區(qū)間，紅色的表示超出期望值的數(shù)據(jù)，圓圈代表殘差的值，豎線代表置信區(qū)間的范圍。輸入>> rcoplot(r,rint)輸出殘差圖見8.10。圖8.10 殘差圖 8.3.2 rstool函數(shù)該函數(shù)是多元二項(xiàng)式回歸函數(shù)，其調(diào)用格式為rstool(x,y, 'model',alpha)其中x為n

28、15;m為矩陣，y為n維列向量，'model'為以下4種模型：'linear'（線性，缺?。?；'interaction'（交叉）：；'quadratic'（完全二次）：；'purequadratic'（純二次）：。alpha為顯著性水平，默認(rèn)值為0.05。例13 設(shè)某商品的需求量與消費(fèi)者的平均收入、商品價格的統(tǒng)計數(shù)據(jù)如表2，建立多元二項(xiàng)式純二次回歸模型，并預(yù)測平均收入為1000、價格為6時的商品需求量。表2 需求量、平均收入和價格統(tǒng)計表需求量10075807050659010011060收入10006001200

29、500300400130011001300300價格5766875439可以直接使用多元二項(xiàng)式回歸，在命令窗口輸入：>> x1=1000 600 1200 500 300 400 1300 1100 1300 300;x2=5 7 6 6 8 7 5 4 3 9;y=100 75 80 70 50 65 90 100 110 60'x=x1' x2'rstool(x,y,'purequadratic')其輸出如圖8.11。圖8.11 多元二項(xiàng)式回歸在圖8.11中x1上面的方框中輸入1000，x2上面的方框中輸入6，在圖形框左側(cè)的“Predic

30、ted Y1”下方的數(shù)據(jù)變?yōu)?8.47981，即預(yù)測出平均收入為1000，價格為6時的商品需求量為88.4791。單擊圖形框左邊Export，則出現(xiàn)圖8.12對話框，可以將回歸參數(shù)beta、剩余標(biāo)準(zhǔn)差rmse和殘差residuals傳送到MATALB的工作區(qū)中。圖8.12輸出對話框在命令窗口輸入：>> betarmseresiduals輸出結(jié)果為：beta = 110.5313 0.1464 -26.5709 -0.0001 1.8475rmse = 4.5362residuals = 5.2724 -0.7162 -4.5158 -1.9390 -3.3315 3.4566 3.

31、4843 -3.4452 -0.0976 1.8320故回歸模型為：，剩余標(biāo)準(zhǔn)差為4.5362, 說明此回歸模型的顯著性較好。還可將該模型轉(zhuǎn)化為多元線性回歸模型，利用regress函數(shù)進(jìn)求解?？奢斎耄?gt;> x1=1000 600 1200 500 300 400 1300 1100 1300 300;x2=5 7 6 6 8 7 5 4 3 9;X=ones(10,1) x1' x2' (x1.2)' (x2.2)'b,bint,r,rint,stats=regress(y,X);b,stats輸出結(jié)果為： b = 110.5313 0.1464 -

32、26.5709 -0.0001 1.8475stats = 0.9702 40.6656 0.0005 20.5771輸出的回歸參數(shù)和rstool函數(shù)的結(jié)果相同。8.3.3 nlinfit函數(shù)和nlintool函數(shù)nlinfit函數(shù)用來確定非線性回歸系數(shù)，調(diào)用格式為：beta, r, J = nlinfit (x, y, 'modelfun', beta0)其中輸入數(shù)據(jù)x,y分別為n×p維矩陣和n維列向量，對于一元非線性回歸，取p=1,即可；'modelfun'為事先定義的非線性回歸函數(shù)的M文件，是回歸系數(shù)beta和x的函數(shù)；beta0是回歸系數(shù)的初值

33、，輸出參數(shù)beta是估計出的回歸系數(shù)，r為殘差，J為Jacobain矩陣。例14 根據(jù)表1中的數(shù)據(jù)，利用nlinfit函數(shù)進(jìn)行非線性回歸，回歸函數(shù)為。首先建立回歸函數(shù)的M文件myfit3.m，內(nèi)容如下：function f = myfit3(beta,xdata) f = beta(1)+beta(2)*exp(-0.02*beta(3)*xdata); 在命令窗口輸入：>> xdata = 100:100:1000;ydata = 1e-03*4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59;beta0 = 0.1 0.1 0.2;beta,r,