助力國賽 | 第2彈 擬合和插值

來自專欄 MATLAB學習11 人贊了文章

前言

上次介紹了Excel和記事本與MATLAB的數據交互。數據讀取之後，就要開始進行分析和處理了，本次將會介紹兩種比較常見的處理方式——擬合和插值。在實際中，需要處理由實驗或測量得到的一些離散數據。插值與擬合方法就是要通過已知數據去確定一些未知數據，或者由離散數據得到連續數據。即尋找某一類已知函數的參數或尋求某個近似函數，使近似函數與已知數據有較高的擬合精度。

由有限個已知數據點，構造一個解析表達式，由此計算數據點之間的函數值，稱之為插值。而擬合是在實驗數據較多的時候，考慮到數據可能存在誤差，因而沒有必要強制曲線通過所有的數據點，只要求近似函數儘可能靠近這些數據點就行了。可以很明顯地看出插值與擬合，一個是要求通過所有給定數據的，一個是不要求。但都是根據一組數據推導另一部分數據的近似值，要求不一樣因而所使用的數學方法也不一樣。一般說來由於實驗所獲得的數據都有些許的誤差，故而擬合比插值在數學建模中使用的較廣泛些。擬合用於比較精確的數據分析，而插值一般用於繪圖特別是二維三維的地形圖。因此，本文將主要介紹擬合相關的操作，插值介紹繪圖操作。

• 擬合
• 插值
• 綜合應用

擬合

曲線擬合也稱曲線逼近，與插值函數有區別，只要求擬合的曲線能合理地反映數據的基本趨勢，而不要求曲線一定通過數據點。一般有幾種不同的判別準則，如使偏差的絕對值之後最小，使偏差的最大絕對值最小和使偏差的平方和最小（即最小二乘法）。常用的方法是最後一種，即基於最小二乘法原理的曲線擬合，我們主要介紹這一種。

一元多項式曲線擬合

如果變數y與x的關係式是n此多項式，即

其中：

是隨機誤差，服從正態分布

，

為擬合係數

那麼稱其為多項式擬合/回歸模型。

其實這個模型是十分有用的，因為根據泰勒公式，任何函數表達式都是可以用多項式來逼近的。多項式擬合有指令語句和圖形窗口兩種方法。

多項式擬合指令

polyfit(X,Y,N)：多項式擬合，返回按降冪排列得多項式係數。

其中，X,Y是數據點的值，分別對象自變數與因變數；N是擬合的最高次冪；P是返回的多項式係數；xi是要求的點的橫坐標。

例如：要對一家公司的工資總額與零售總額進行曲線擬合分析,並預測當工資總額為80時，零售總額y的值：。

作為一個回歸模型來進行分析，我們應當首先做出散點圖觀看x，y是否有具有一定的關係。

從圖中看出，可以進行一元線性擬合。

code:

x=[23.80,27.60,31.60,32.40,33.70,34.90,43.20,52.80,63.80,73.40];y=[41.4,51.8,61.70,67.90,68.70,77.50,95.90,137.40,155.0,175.0];% 作散點圖plot(x,y,r*) %作散點圖xlabel(x(職工工資總額)) %橫坐標名ylabel(y(商品零售總額)) %縱坐標名%數據擬合及預測p=polyfit(x,y,1) % 一元線性擬合x0=80;y0=polyval(p,x0) %預測x0=80%繪製曲線擬合圖x2=[23:0.5:90];y2=polyval(p,x2);hold onplot(x2,y2,b-)

效果圖：

圖形窗口的多項式擬合

在圖像窗口中可以用菜單方式對數據進行簡單、快捷、高效的擬合，用以下數據為例，先畫出散點圖：

>> x = [1 2 3 4 5 6 7 8 9];>> y = [9 7 6 3 -1 2 5 7 20];>> plot(x,y,r*);

然後在圖形窗口單擊工具（T） →基本擬合（F），打開對話框，按如下操作，就可以得到相應效果，其中分布進行了線性、二階、三階對數據進行了多項式擬合，並可以查看相應殘差，相互對比得出最好的擬合效果圖，但是一定注意不要出現過擬合現象。

在對比擬合效果之後，可以發現，在該題中三次曲線擬合效果較好。

制定函數擬合

像這一類曲線擬合，一般是已知函數表達式的，所要做的工作就是確定其中的相關係數而已，一種比較簡單的想法就是通過函數變化將其轉換為一般的多項式進行擬合，不過MATLAB支持用戶自定義函數擬合，那就通過一個例子來簡單說明一下。

先繪製出散點圖：

code:

x = [0;0.4;1.2;2;2.8;3.6;4.4;5.2;6;7.2;8;9.2;10.4;11.6;12.4;13.6;14.4;15];y = [1;0.85;0.29;-0.27;-0.53;-0.4;-0.12;0.17;0.28;0.15;-0.03;-0.15;-0.071;0.059;0.08;0.032;-0.015;-0.02];plot(x,y,r*);

我們知道阻尼振蕩的函數形式為

，則可以用MATLAB進行擬合。

code:

syms t;x = [0;0.4;1.2;2;2.8;3.6;4.4;5.2;6;7.2;8;9.2;10.4;11.6;12.4;13.6;14.4;15];y = [1;0.85;0.29;-0.27;-0.53;-0.4;-0.12;0.17;0.28;0.15;-0.03;-0.15;-0.071;0.059;0.08;0.032;-0.015;-0.02];plot(x,y,r*);f = fittype(a*cos(k*t)*exp(w*t),independent,t,coefficients,{a,k,w});cfun = fit(x,y,f); %顯示擬合函數xi = 0:.1:20;yi = cfun(xi);plot(x,y,r*,xi,yi,b-);

效果：

程序中，函數fittype是自定義擬合函數；cfun = fit(x,y,f)是工具自定義的擬合函數f來擬合數據的，數據應該為列向量形式。從結果上來看，擬合效果還是可以的，並且給出了各參數的置信區間。至於出現的警告是由於沒有給a，k，w3個參數的初始值，如果擬合效果不理想，可以多運行幾次。

擬合曲線工具箱

MATLAB提供了曲線擬合工具箱，功能十分強大，十分方便。在MATLAB R2016b中只需要依次單擊APP → Curve Fitting即可以。或者在命令行窗口輸入cftool，也可以打開。

下面使用一組數據進行一個簡單的演示：

數據：

x = [1 2 3 4 5 6 7 8 9];y = [9 7 6 3 -1 2 5 7 20];

現在命令行窗口輸入上述數據，然後按如下操作，還是很簡單的。

多元線性回歸

我們上面所討論的都是一元的問題，但如果自變數不再是一個而是兩個或者三個的時候，那又如何處理呢？很顯然多元線性回歸比一元線性回歸要複雜，但強大的MATLAB還是提供了相應的處理函數。

我們通過一個例子來說明其應用。

一個公司的銷售與庫存資金額、廣告投入和員工薪酬總額有關，根據下面數據，試分析出它們之間的相互關係。

處理多元線性回歸問題一般有以下四個步驟：

(1) 分析問題，選擇因變數y與解釋變數X，作出y和xi的散點圖，初步設定多元線性回歸模型的參數個數；

(2)調用多元線性回歸函數，得到部分結果。[b,bint,r,rint,s]=regress(y,X,alpha)

(3)繪製殘差圖rcoplot，分析異常點，改進模型。

(4) 顯著性檢驗，若檢驗通過，則用模型作預測。

其中多元線性回歸建模命令：regeress的一般調用形式為：

b = regress(y,X)

輸入量y表示因變數的觀測向量,X是自變數矩陣，第1列元全部是 「1」，第j列是Xj的觀測向量。向量b是回歸係數估計值

先分別繪製散點圖：

clearclc%錄入數據A=[75.2 30.6 21.1 1090.4;77.6 31.3 21.4 113380.7 33.9 22.9 1242.1;76 29.6 21.4 1003.279.5 32.5 21.5 1283.2;81.8 27.9 21.7 1012.298.3 24.8 21.5 1098.8;67.7 23.6 21 826.374 33.9 22.4 1003.3;151 27.7 24.7 1554.690.8 45.5 23.2 1199;102.3 42.6 24.3 1483.1115.6 40 23.1 1407.1;125 45.8 29.1 1551.3137.8 51.7 24.6 1601.2;175.6 67.2 27.5 2311.7155.2 65 26.5 2126.7;174.3 65.4 26.8 2256.5];%分別繪製散點圖[m,n]=size(A); subplot(3,1,1);plot(A(:,1),A(:,4),+);xlabel(x1(庫存資金額)) ;ylabel(y(銷售額)) subplot(3,1,2);plot(A(:,2),A(:,4),*),;xlabel(x2(廣告投入)) ;ylabel(y(銷售額));subplot(3,1,3);plot(A(:,3),A(:,4),x);xlabel(x3(員工薪酬));ylabel(y(銷售額));

從圖中看出可以進行多元線性回歸分析。

然後調用regeress函數

x=[ones(m,1), A(:,1), A(:,2), A(:,3)];y=A(:,4);b=regress(y,x);

可以得到b的值如下，也就是對應的各變數的係數和常數項。

也就是說回歸模型為：

但是我們還需要進一步分析，來得出更精確的模型。

將上述調用regeress函數的代碼改為：

>> [b,bint,r,rint,stats]=regress(y,x);b,bint,stats, % 輸出結果b = 162.0632 7.2739 13.9575 -4.3996bint = -580.3603 904.4867 4.3734 10.1743 7.1649 20.7501 -46.7796 37.9805stats = 1.0e+04 * 0.0001 0.0105 0.0000 1.0078

其中bint的各行分別為各參數的95%的置信區間。

繪製殘差圖，踢出異常點，然後重新分析。

在原有基礎上輸入下述代碼：

>> x=[ones(m,1), A(:,1), A(:,2), A(:,3)];y=A(:,4)[b,bint,r,rint,stats]=regress(y,x);b,bint,statsrcoplot(r,rint)y = 1.0e+03 * 1.0904 1.1330 1.2421 1.0032 1.2832 1.0122 1.0988 0.8263 1.0033 1.5546 1.1990 1.4831 1.4071 1.5513 1.6012 2.3117 2.1267 2.2565b = 162.0632 7.2739 13.9575 -4.3996bint = -580.3603 904.4867 4.3734 10.1743 7.1649 20.7501 -46.7796 37.9805stats = 1.0e+04 * 0.0001 0.0105 0.0000 1.0078

繪製出殘差圖：

發現第5個數據點異常，踢出掉，重複上述操作，得出合理模型。

擬合優度檢驗

我們通過MATLAB擬合出來了函數，那麼如何判斷其效果呢？

插值

插值是一種古老的數學方法，它是要求過所有數據的，在MATLAB中分為一維插值和二維插值。

一維插值

當被插值的函數為一元函數時，為一維插值。MATLAB使用interp1函數來實現一維插值。其調用格式如下：

Vq = interp1(X,V,Xq,METHOD);

其中，X為自變數的取值範圍；V為函數值，；Xq為插值點向量或數組；METHOD是字元串變數，用來設定插值方法。

MATLAB提供了以下幾種插值方法：

• METHOD = 『nearest』：臨近點插值。插值點函數值估計為該插值點最近的數據點函數值。
• METHOD = 『linear』：線性插值。根據相鄰數據點的線性函數估計落在該區域內的插值數據點的函數值。
• METHOD = 『spline』：三次樣條插值。這種方法在相鄰數據點間建立三次多項式函數，根據多項式函數確定數據點的函數值。
• METHOD = 『pchip』或『cubic』：立方插值。通過分段立方Hermite插值方法計算結果。

在上述的插值方法中，三次樣條插值最為常用，雖運行時間最長，但平滑性最好。

二維插值

當被插值的函數為二元函數的時，為二維插值。MATLAB用interp2函數來實現。其調用格式如下：

Vq = interp2(X,Y,V,Xq,Yq,METHOD);

其中，X、Y、V是具有相同大小的矩陣；V(i,j)是數據點[X(i,j),Y(i,j)]上的函數值；Xq，Yq為待插值數據網路；METHOD是字元串變數，用來設定插值方法。

MATLAB提供了以下幾種插值方法：

• METHOD = 『nearest』：臨近點插值。將插值點周圍4個數據點中離該插值點最近的數據點作為該插值點函數值的估計值。
• METHOD = 『linear』：雙線性插值。是MATLAB的interp2函數的默認使用的插值方法。該方法將插值點周圍4個數據點的函數值的線性組合作為插值點的函數值的估計值
• METHOD = 『spline』：三次樣條插值。該方法的計算效率高，得到的曲面光滑，是經常使用的插值方法。
• METHOD = 『cubic』：雙立方插值。該方法利用插值點周圍16個數據點，相對於前兩個插值方法，需要消耗更多的內存和計算時間，因而計算效率不高，但是所得到的曲面更加光滑。

給出如下一個例子：

clear all[X,Y] = meshgrid(-3:0.25:3);Z = peaks(X,Y);[XI,YI] = meshgrid(-3:.125:3);ZI = interp2(X,Y,Z,XI,YI);mesh(X,Y,Z),hold,mesh(XI,YI,ZI+15)hold offaxis([-3 3 -3 3 -5 20])

綜合應用

人口模型

下表是某地區1971-2000年的人口數據，是給出該地區人口增長的數據模型。

通過散點圖，選取Logistic模型來進行曲線擬合，其基本形式為：

MATLAB求解如下：

clearclc% 讀入人口數據（1971－2000年）Y=[33815 33981 34004 34165 34212 34327 34344 34458 34498 34476 34483 34488 34513 34497 34511 34520 34507 34509 34521 34513 34515 34517 34519 34519 34521 34521 34523 34525 34525 34527]% 讀入時間變數數據（t＝年份－1970）T=[1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30]% 線性化處理for t = 1:30, x(t)=exp(-t); y(t)=1/Y(t);end% 計算，並輸出回歸係數Bc=zeros(30,1)+1;X=[c,x];B=inv(X*X)*X*yfor i=1:30,% 計算回歸擬合值 z(i)=B(1,1)+B(2,1)*x(i);% 計算離差 s(i)=y(i)-sum(y)/30;% 計算誤差 w(i)=z(i)-y(i);end% 計算離差平方和SS=s*s;% 回歸誤差平方和QQ=w*w;% 計算回歸平方和UU=S-Q;% 計算，並輸出F檢驗值F=28*U/Q% 計算非線性回歸模型的擬合值for j=1:30, Y(j)=1/(B(1,1)+B(2,1)*exp(-j));end% 輸出非線性回歸模型的擬合曲線（Logisic曲線）plot(T,Y)

繪製地形區

利用插值函數來繪製地形圖。

插值代碼如下：

clearclc[x,y] = meshgrid(1:10);h = [0,0.02,-0.12,0,-2.09,0,-0.58,-0.08,0,0;... 0.02,0,0,-2.38,0,-4.96,0,0,0,-0.1;... 0,0.10,1.00,0,-3.04,0,-0.53,0,0.10,0;... 0,0,0,3.52,0,0,0,0,0,0;... -0.43,-1.98,0,0,0,0.77,0,2.17,0,0;... 0,0,-2.29,0,0.69,0,2.59,0,0.30,0;... -0.09,-0.31,0,0,0,4.27,0,0,0,-0.01;... 0,0,0,5.13,7.40,0,1.89,0,0.04,0;... 0.1,0,0.58,0,0,1.75,0,-0.11,0,0;... 0,-0.01,0,0,0.3,0,0,0,0,0.01];[xi,yi] = meshgrid(1:0.1:10);hi = interp2(x,y,h,xi,yi,spline);surf(hi);xlabel(x),ylabel(y),zlabel(h);

再看一個插值地形圖的例子：

x=0:.5:5;y=0:.5:6;z=[89 90 87 85 92 91 96 93 90 87 82 92 96 98 99 95 91 89 86 84 82 84 96 98 95 92 90 88 85 84 83 81 85 80 81 82 89 95 96 93 92 89 86 86 82 85 87 98 99 96 97 88 85 82 83 82 85 89 94 95 93 92 91 86 84 88 88 92 93 94 95 89 87 86 83 81 92 92 96 97 98 96 93 95 84 82 81 84 85 85 81 82 80 80 81 85 90 93 95 84 86 81 98 99 98 97 96 95 84 87 80 81 85 82 83 84 87 90 95 86 88 80 82 81 84 85 86 83 82 81 80 82 87 88 89 98 99 97 96 98 94 92 87];mesh(x,y,z) %繪原始數據圖xi=linspace(0,5,50); %加密橫坐標數據到50個yi=linspace(0,6,80); %加密縱坐標數據到60個[xii,yii]=meshgrid(xi,yi); %生成網格數據zii=interp2(x,y,z,xii,yii,cubic); %插值mesh(xii,yii,zii) %加密後的地貌圖hold on % 保持圖形[xx,yy]=meshgrid(x,y); %生成網格數據plot3(xx,yy,z+0.1,ob) %原始數據用『O』繪出

編輯不易，歡迎推廣