數學建模專欄 | 第三篇：MATLAB數據建模方法(上) —常用方法
作 者 簡 介
卓金武，MathWorks中國高級工程師，教育業務經理，在數據分析、數據挖掘、機器學習、數學建模、量化投資和優化等科學計算方面有多年工作經驗，現主要負責MATLAB校園版業務。曾2次獲全國大學生數學建模競賽一等獎，1次獲全國研究生數學建模競賽一等獎。
以數據為基礎而建立數學模型的方法稱為數據建模方法， 包括回歸、統計、機器學習、深度學習、灰色預測、主成分分析、神經網絡、時間序列分析等方法， 其中最常用的方法還是回歸方法。 本講主要介紹在數學建模中常用幾種回歸方法的 MATLAB 實現過程。
根據回歸方法中因變量的個數和回歸函數的類型（線性或非線性）可將回歸方法分為:一元線性、一元非線性、多元回歸。另外還有兩種特殊的回歸方式，一種在回歸過程中可以調整變量數的回歸方法，稱為逐步回歸，另一種是以指數結構函數作為回歸模型的回歸方法，稱為 Logistic 回歸。本講將逐一介紹這幾個回歸方法。
1. 一元回歸
1.1 一元線性回歸
[ 例1 ] 近 10 年來，某市社會商品零售總額與職工工資總額（單位：億元）的數據見表3-1，請建立社會商品零售總額與職工工資總額數據的回歸模型。
表1 商品零售總額與職工工資總額
該問題是典型的一元回歸問題，但先要確定是線性還是非線性，然后就可以利用對應的回歸方法建立他們之間的回歸模型了，具體實現的 MATLAB 代碼如下：
（1）輸入數據
clc, clear all, close all
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];
（2）采用最小二乘回歸
Figure
plot(x,y,'r*') %作散點圖
xlabel('x（職工工資總額）','fontsize', 12) %橫坐標名
ylabel('y（商品零售總額）', 'fontsize',12) %縱坐標名
set(gca,'linewidth',2);
% 采用最小二乘擬合
Lxx=sum((x-mean(x)).^2);
Lxy=sum((x-mean(x)).*(y-mean(y)));
b1=Lxy/Lxx;
b0=mean(y)-b1*mean(x);
y1=b1*x+b0;
hold on
plot(x, y1,'linewidth',2);
運行本節程序，會得到如圖 1 所示的回歸圖形。在用最小二乘回歸之前，先繪制了數據的散點圖，這樣就可以從圖形上判斷這些數據是否近似成線性關系。當發現它們的確近似在一條線上后，再用線性回歸的方法進行回歸，這樣也更符合我們分析數據的一般思路。
圖1 職工工資總額和商品零售總額關系趨勢圖
（3）采用 LinearModel.fit 函數進行線性回歸
m2 = LinearModel.fit(x,y)
運行結果如下：
m2 =
Linear regression model:
y ~ 1 + x1
Estimated Coefficients:Estimate SE tStat pValue
(Intercept) -23.549 5.1028 -4.615 0.0017215
x1 2.7991 0.11456 24.435 8.4014e-09
R-squared: 0.987, Adjusted R-Squared 0.985
F-statistic vs. constant model: 597, p-value = 8.4e-09
（4）采用 regress 函數進行回歸
Y=y';
X=[ones(size(x,2),1),x'];
[b, bint, r, rint, s] = regress(Y, X)
運行結果如下：
b =
-23.5493
2.7991
在以上回歸程序中，使用了兩個回歸函數 LinearModel.fit 和 regress。在實際使用中，只要根據自己的需要選用一種就可以了。函數 LinearModel.fit 輸出的內容為典型的線性回歸的參數。關於 regress，其用法多樣，MATLAB 幫助中關於 regress 的用法，有以下幾種：
b = regress(y,X)
[b,bint] = regress(y,X)
[b,bint,r] = regress(y,X)
[b,bint,r,rint] = regress(y,X)
[b,bint,r,rint,stats] = regress(y,X)
[...] = regress(y,X,alpha)
輸入 y（因變量，列向量），X（1與自變量組成的矩陣）和（alpha，是顯著性水平, 缺省時默認0.05）。
輸出
bint 是 β0，β1 的置信區間，r 是殘差（列向量），rint是殘差的置信區間，s包含4個統計量：決定系數 R^2（相關系數為R），F 值，F(1,n-2) 分布大於 F 值的概率 p，剩余方差 s^2 的值。也可由程序 sum(r^2)/(n-2) 計算。其意義和用法如下：R^2 的值越接近 1，變量的線性相關性越強，說明模型有效；如果滿足
則認為變量y與x顯著地有線性關系，其中 F1-α(1，n-2) 的值可查F分布表，或直接用 MATLAB 命令 finv(1-α,1, n-2) 計算得到；如果 p<α 表示線性模型可用。這三個值可以相互印證。s^2 的值主要用來比較模型是否有改進，其值越小說明模型精度越高。
1.2 一元非線性回歸
在一些實際問題中，變量間的關系並不都是線性的，此時就應該用非線性回歸。用非線性回歸首先要解決的問題是回歸方程中的參數如何估計。下面通過一個實例來說明如何利用非線性回歸技術解決實例的問題。
[ 例2 ] 為了解百貨商店銷售額 x 與流通率（這是反映商業活動的一個質量指標，指每元商品流轉額所分攤的流通費用）y 之間的關系，收集了九個商店的有關數據（見表2）。請建立它們關系的數學模型。
表2 銷售額與流通費率數據
圖2 銷售額與流通費率之間的關系圖
為了得到 x 與 y 之間的關系，先繪制出它們之間的散點圖，如圖 2 所示的“雪花”點圖。由該圖可以判斷它們之間的關系近似為對數關系或指數關系，為此可以利用這兩種函數形式進行非線性擬合，具體實現步驟及每個步驟的結果如下：
（1）輸入數據
clc, clear all, close all
x=[1.5, 4.5, 7.5,10.5,13.5,16.5,19.5,22.5,25.5];
y=[7.0,4.8,3.6,3.1,2.7,2.5,2.4,2.3,2.2];
plot(x,y,'*','linewidth',2);
set(gca,'linewidth',2);
xlabel('銷售額x/萬元','fontsize', 12)
ylabel('流通費率y/%', 'fontsize',12)
（2）對數形式非線性回歸
m1 = @(b,x) b(1) + b(2)*log(x);
nonlinfit1 = fitnlm(x,y,m1,[0.01;0.01])
b=nonlinfit1.Coefficients.Estimate;
Y1=b(1,1)+b(2,1)*log(x);
hold on
plot(x,Y1,'--k','linewidth',2)
運行結果如下：
nonlinfit1 =
Nonlinear regression model:
y ~ b1 + b2*log(x)
Estimated Coefficients:
Estimate SE tStat pValue
b1 7.3979 0.26667 27.742 2.0303e-08
b2 -1.713 0.10724 -15.974 9.1465e-07
R-Squared: 0.973, Adjusted R-Squared 0.969
F-statistic vs. constant model: 255, p-value = 9.15e-07
（3）指數形式非線性回歸
m2 = 'y ~ b1*x^b2';
nonlinfit2 = fitnlm(x,y,m2,[1;1])
b1=nonlinfit2.Coefficients.Estimate(1,1);
b2=nonlinfit2.Coefficients.Estimate(2,1);
Y2=b1*x.^b2;
hold on
plot(x,Y2,'r','linewidth',2)
legend('原始數據','a+b*lnx','a*x^b')
運行結果如下：
nonlinfit2 =
Nonlinear regression model:
y ~ b1*x^b2
Estimated Coefficients:
Estimate SE tStat pValue
b1 8.4112 0.19176 43.862 8.3606e-10
b2 -0.41893 0.012382 -33.834 5.1061e-09
R-Squared: 0.993, Adjusted R-Squared 0.992
F-statistic vs. zero model: 3.05e+03, p-value = 5.1e-11
在該案例中，選擇兩種函數形式進行非線性回歸，從回歸結果來看，對數形式的決定系數為 0.973 ，而指數形式的為 0.993 ，優於前者，所以可以認為指數形式的函數形式更符合 y 與 x 之間的關系，這樣就可以確定他們之間的函數關系形式了。
2. 多元回歸
[ 例3 ] 某科學基金會希望估計從事某研究的學者的年薪 Y 與他們的研究成果（論文、著作等）的質量指標 X1、從事研究工作的時間 X2、能成功獲得資助的指標 X3 之間的關系，為此按一定的實驗設計方法調查了 24 位研究學者，得到如表3 所示的數據（ i 為學者序號），試建立 Y 與 X1 , X2 , X3 之間關系的數學模型，並得出有關結論和作統計分析。
表3 從事某種研究的學者的相關指標數據
該問題是典型的多元回歸問題，但能否應用多元線性回歸，最好先通過數據可視化判斷他們之間的變化趨勢，如果近似滿足線性關系，則可以執行利用多元線性回歸方法對該問題進行回歸。具體步驟如下：
（1）作出因變量 Y 與各自變量的樣本散點圖
作散點圖的目的主要是觀察因變量 Y 與各自變量間是否有比較好的線性關系，以便選擇恰當的數學模型形式。圖3 分別為年薪 Y 與成果質量指標 X1、研究工作時間 X2、獲得資助的指標 X3 之間的散點圖。從圖中可以看出這些點大致分布在一條直線旁邊，因此，有比較好的線性關系，可以采用線性回歸。繪制圖3的代碼如下：
subplot(1,3,1),plot(x1,Y,'g*'),
subplot(1,3,2),plot(x2,Y,'k+'),
subplot(1,3,3),plot(x3,Y,'ro'),
圖3 因變量Y與各自變量的樣本散點圖
（2）進行多元線性回歸
這里可以直接使用 regress 函數執行多元線性回歸，具體代碼如下：
x1=[3.5 5.3 5.1 5.8 4.2 6.0 6.8 5.5 3.1 7.2 4.5 4.9 8.0 6.5 6.5 3.7 6.2 7.0 4.0 4.5 5.9 5.6 4.8 3.9];
x2=[9 20 18 33 31 13 25 30 5 47 25 11 23 35 39 21 7 40 35 23 33 27 34 15];
x3=[6.1 6.4 7.4 6.7 7.5 5.9 6.0 4.0 5.8 8.3 5.0 6.4 7.6 7.0 5.0 4.0 5.5 7.0 6.0 3.5 4.9 4.3 8.0 5.0];
Y=[33.2 40.3 38.7 46.8 41.4 37.5 39.0 40.7 30.1 52.9 38.2 31.8 43.3 44.1 42.5 33.6 34.2 48.0 38.0 35.9 40.4 36.8 45.2 35.1];
n=24; m=3;
X=[ones(n,1),x1',x2',x3'];
[b,bint,r,rint,s]=regress(Y',X,0.05);
運行后即得到結果如表4所示。
表4 對初步回歸模型的計算結果
計算結果包括回歸系數 b = (β0,β1,β2,β3) = (18.0157, 1.0817, 0.3212, 1.2835)、回歸系數的置信區間，以及統計變量 stats（它包含四個檢驗統計量：相關系數的平方R^2，假設檢驗統計量 Ｆ，與 F 對應的概率 p，s^2 的值）。因此我們得到初步的回歸方程為：
由結果對模型的判斷：
回歸系數置信區間不包含零點表示模型較好，殘差在零點附近也表示模型較好，接着就是利用檢驗統計量 Ｒ，Ｆ，p 的值判斷該模型是否可用。
1）相關系數 Ｒ 的評價：本例 Ｒ 的絕對值為 0.9542 ，表明線性相關性較強。
2）F 檢驗法：當 F > F1-α(m,n-m-1) ，即認為因變量 y 與自變量 x1,x2,...,xm 之間有顯著的線性相關關系；否則認為因變量 y 與自變量 x1,x2,...,xm 之間線性相關關系不顯著。本例 Ｆ＝67.919 > F1-0.05( 3,20 ) = 3.10。
3）p 值檢驗：若 p < α（α 為預定顯著水平），則說明因變量 y 與自變量 x1,x2,...,xm之間顯著地有線性相關關系。本例輸出結果，p<0.0001,顯然滿足 p<α=0.05。
以上三種統計推斷方法推斷的結果是一致的，說明因變量 y 與自變量之間顯著地有線性相關關系，所得線性回歸模型可用。s^2 當然越小越好，這主要在模型改進時作為參考。
3. 逐步歸回
[ 例4 ] （Hald,1960）Hald 數據是關於水泥生產的數據。某種水泥在凝固時放出的熱量 Y（單位：卡/克）與水泥中 4 種化學成品所占的百分比有關：
在生產中測得 12 組數據，見表5，試建立 Y 關於這些因子的“最優”回歸方程。
表5 水泥生產的數據
對於例 4 中的問題，可以使用多元線性回歸、多元多項式回歸，但也可以考慮使用逐步回歸。從逐步回歸的原理來看，逐步回歸是以上兩種回歸方法的結合，可以自動使得方程的因子設置最合理。對於該問題，逐步回歸的代碼如下：
X=[7,26,6,60;1,29,15,52;11,56,8,20;11,31,8,47;7,52,6,33;11,55,9,22;3,71,17,6;1,31,22,44;2,54,18,22;21,47,4,26;1,40,23,34;11,66,9,12]; %自變量數據
Y=[78.5,74.3,104.3,87.6,95.9,109.2,102.7,72.5,93.1,115.9,83.8,113.3]; %因變量數據
Stepwise（X,Y,[1,2,3,4],0.05,0.10）% in=[1,2,3,4]表示X1、X2、X3、X4均保留在模型中
程序執行后得到下列逐步回歸的窗口，如圖 4 所示。
圖4 逐步回歸操作界面
在圖 4 中，用藍色行顯示變量 X1、X2、X3、X4 均保留在模型中，窗口的右側按鈕上方提示：將變量X3剔除回歸方程（Move X3 out），單擊 Next Step 按鈕，即進行下一步運算，將第 3 列數據對應的變量 X3 剔除回歸方程。單擊 Next Step 按鈕后，剔除的變量 X3 所對應的行用紅色表示，同時又得到提示：將變量 X4 剔除回歸方程（Move X4 out），單擊 Next Step 按鈕，這樣一直重復操作，直到 “Next Step” 按鈕變灰，表明逐步回歸結束，此時得到的模型即為逐步回歸最終的結果。
4. Logistic 回歸
[ 例5 ] 企業到金融商業機構貸款，金融商業機構需要對企業進行評估。評估結果為 0 , 1 兩種形式，0 表示企業兩年后破產，將拒絕貸款，而 1 表示企業 2 年后具備還款能力，可以貸款。在表 6 中，已知前 20 家企業的三項評價指標值和評估結果，試建立模型對其他 5 家企業（企業 21-25）進行評估。
表6 企業還款能力評價表
對於該問題，很明顯可以用 Logistic 模型來回歸，具體求解程序如下：
% logistic回歸Matlab實現程序
%% 數據准備
clc, clear, close all
X0=xlsread('logistic_ex1.xlsx', 'A2:C21'); % 回歸模型的輸入
Y0=xlsread('logistic_ex1.xlsx', 'D2:D21'); % 回歸模型的輸出
X1=xlsread('logistic_ex1.xlsx', 'A2:C26'); % 預測數據輸入
%% logistics函數
GM = fitglm(X0,Y0,'Distribution','binomial');
Y1 = predict(GM,X1);
%% 模型的評估
N0 =1:size(Y0,1); N1= 1:size(Y1,1);
plot(N0', Y0, '-kd');
hold on; scatter(N1', Y1, 'b')
xlabel('數據點編號'); ylabel('輸出值');
得到的回歸結果與原始數據的比較如圖5所示。
圖5 回歸結果與原始數據的比較圖
5. 小結
本講主要介紹數學建模中常用的幾種回歸方法。在使用回歸方法的時候，首先可以判斷自變量的個數，如果超過 2 個，則需要用到多元回歸的方法，否則考慮用一元回歸。然后判斷是線性還是非線性，這對於一元回歸是比較容易的，而對於多元，往往是將其他變量保持不變，將多元轉化為一元再去判斷是線性還是非線性。如果變量很多，而且復雜，則可以首先考慮多元線性回歸，檢驗回歸效果，也可以用逐步回歸。總之，用回歸方法比較靈活，根據具體情景還是比較容易找到合適的方法的。
往期 | 數模專欄
Scan QR Code via WeChat
to follow Official Account