以數(shù)據(jù)為基礎(chǔ)而建立數(shù)學(xué)模型的方法稱為數(shù)據(jù)建模方法， 包括回歸、統(tǒng)計(jì)、機(jī)器學(xué)習(xí)、深度學(xué)習(xí)、灰色預(yù)測(cè)、主成分分析、神經(jīng)網(wǎng)絡(luò)、時(shí)間序列分析等方法， 其中最常用的方法還是回歸方法。 本講主要介紹在數(shù)學(xué)建模中常用幾種回歸方法的 MATLAB 實(shí)現(xiàn)過程。

根據(jù)回歸方法中因變量的個(gè)數(shù)和回歸函數(shù)的類型（線性或非線性）可將回歸方法分為:一元線性、一元非線性、多元回歸。另外還有兩種特殊的回歸方式，一種在回歸過程中可以調(diào)整變量數(shù)的回歸方法，稱為逐步回歸，另一種是以指數(shù)結(jié)構(gòu)函數(shù)作為回歸模型的回歸方法，稱為 Logistic 回歸。本講將逐一介紹這幾個(gè)回歸方法。

1.一元回歸

1.1一元線性回歸

[ 例1 ]近 10 年來，某市社會(huì)商品零售總額與職工工資總額（單位：億元）的數(shù)據(jù)見表3-1，請(qǐng)建立社會(huì)商品零售總額與職工工資總額數(shù)據(jù)的回歸模型。

表1 商品零售總額與職工工資總額

該問題是典型的一元回歸問題，但先要確定是線性還是非線性，然后就可以利用對(duì)應(yīng)的回歸方法建立他們之間的回歸模型了，具體實(shí)現(xiàn)的 MATLAB 代碼如下：

（1）輸入數(shù)據(jù)

clc, clear all, closeall

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*')%作散點(diǎn)圖

xlabel('x（職工工資總額）','fontsize', 12)%橫坐標(biāo)名

ylabel('y（商品零售總額）', 'fontsize',12)%縱坐標(biāo)名

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

運(yùn)行本節(jié)程序，會(huì)得到如圖 1 所示的回歸圖形。在用最小二乘回歸之前，先繪制了數(shù)據(jù)的散點(diǎn)圖，這樣就可以從圖形上判斷這些數(shù)據(jù)是否近似成線性關(guān)系。當(dāng)發(fā)現(xiàn)它們的確近似在一條線上后，再用線性回歸的方法進(jìn)行回歸，這樣也更符合我們分析數(shù)據(jù)的一般思路。

圖1職工工資總額和商品零售總額關(guān)系趨勢(shì)圖

（3）采用 LinearModel.fit 函數(shù)進(jìn)行線性回歸

m2 = LinearModel.fit(x,y)

運(yùn)行結(jié)果如下：

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 函數(shù)進(jìn)行回歸

Y=y';

X=[ones(size(x,2),1),x'];

[b, bint, r, rint, s] = regress(Y, X)

運(yùn)行結(jié)果如下：

b =

-23.5493

2.7991

在以上回歸程序中，使用了兩個(gè)回歸函數(shù)LinearModel.fit 和 regress。在實(shí)際使用中，只要根據(jù)自己的需要選用一種就可以了。函數(shù) LinearModel.fit 輸出的內(nèi)容為典型的線性回歸的參數(shù)。關(guān)于 regress，其用法多樣，MATLAB 幫助中關(guān)于 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，是顯著性水平, 缺省時(shí)默認(rèn)0.05）。

輸出

bint 是 β0，β1的置信區(qū)間，r 是殘差（列向量），rint是殘差的置信區(qū)間，s包含4個(gè)統(tǒng)計(jì)量：決定系數(shù) R^2（相關(guān)系數(shù)為R），F(xiàn) 值，F(xiàn)(1,n-2) 分布大于 F 值的概率 p，剩余方差 s^2 的值。也可由程序 sum(r^2)/(n-2) 計(jì)算。其意義和用法如下：R^2 的值越接近 1，變量的線性相關(guān)性越強(qiáng)，說明模型有效；如果滿足

則認(rèn)為變量y與x顯著地有線性關(guān)系，其中 F1-α(1，n-2)的值可查F分布表，或直接用 MATLAB 命令 finv(1-α,1, n-2) 計(jì)算得到；如果 p<α 表示線性模型可用。這三個(gè)值可以相互印證。s^2 的值主要用來比較模型是否有改進(jìn)，其值越小說明模型精度越高。

1.2一元非線性回歸

在一些實(shí)際問題中，變量間的關(guān)系并不都是線性的，此時(shí)就應(yīng)該用非線性回歸。用非線性回歸首先要解決的問題是回歸方程中的參數(shù)如何估計(jì)。下面通過一個(gè)實(shí)例來說明如何利用非線性回歸技術(shù)解決實(shí)例的問題。

[ 例2 ]為了解百貨商店銷售額 x 與流通率（這是反映商業(yè)活動(dòng)的一個(gè)質(zhì)量指標(biāo)，指每元商品流轉(zhuǎn)額所分?jǐn)偟牧魍ㄙM(fèi)用）y 之間的關(guān)系，收集了九個(gè)商店的有關(guān)數(shù)據(jù)（見表2）。請(qǐng)建立它們關(guān)系的數(shù)學(xué)模型。

表2銷售額與流通費(fèi)率數(shù)據(jù)

圖2銷售額與流通費(fèi)率之間的關(guān)系圖

為了得到 x 與 y 之間的關(guān)系，先繪制出它們之間的散點(diǎn)圖，如圖 2 所示的“雪花”點(diǎn)圖。由該圖可以判斷它們之間的關(guān)系近似為對(duì)數(shù)關(guān)系或指數(shù)關(guān)系，為此可以利用這兩種函數(shù)形式進(jìn)行非線性擬合，具體實(shí)現(xiàn)步驟及每個(gè)步驟的結(jié)果如下：

（1）輸入數(shù)據(jù)

clc, clear all, closeall

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('流通費(fèi)率y/%', 'fontsize',12)

（2）對(duì)數(shù)形式非線性回歸

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)

運(yùn)行結(jié)果如下：

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）指數(shù)形式非線性回歸

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('原始數(shù)據(jù)','a+b*lnx','a*x^b')

運(yùn)行結(jié)果如下：

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

在該案例中，選擇兩種函數(shù)形式進(jìn)行非線性回歸，從回歸結(jié)果來看，對(duì)數(shù)形式的決定系數(shù)為 0.973 ，而指數(shù)形式的為 0.993 ，優(yōu)于前者，所以可以認(rèn)為指數(shù)形式的函數(shù)形式更符合 y 與 x 之間的關(guān)系，這樣就可以確定他們之間的函數(shù)關(guān)系形式了。

2.多元回歸

[ 例3 ]某科學(xué)基金會(huì)希望估計(jì)從事某研究的學(xué)者的年薪 Y 與他們的研究成果（論文、著作等）的質(zhì)量指標(biāo) X1、從事研究工作的時(shí)間 X2、能成功獲得資助的指標(biāo) X3之間的關(guān)系，為此按一定的實(shí)驗(yàn)設(shè)計(jì)方法調(diào)查了 24 位研究學(xué)者，得到如表3 所示的數(shù)據(jù)（ i 為學(xué)者序號(hào)），試建立 Y 與 X1, X2, X3之間關(guān)系的數(shù)學(xué)模型，并得出有關(guān)結(jié)論和作統(tǒng)計(jì)分析。

表3從事某種研究的學(xué)者的相關(guān)指標(biāo)數(shù)據(jù)

該問題是典型的多元回歸問題，但能否應(yīng)用多元線性回歸，最好先通過數(shù)據(jù)可視化判斷他們之間的變化趨勢(shì)，如果近似滿足線性關(guān)系，則可以執(zhí)行利用多元線性回歸方法對(duì)該問題進(jìn)行回歸。具體步驟如下：

（1）作出因變量 Y 與各自變量的樣本散點(diǎn)圖

作散點(diǎn)圖的目的主要是觀察因變量 Y 與各自變量間是否有比較好的線性關(guān)系，以便選擇恰當(dāng)?shù)臄?shù)學(xué)模型形式。圖3 分別為年薪 Y 與成果質(zhì)量指標(biāo) X1、研究工作時(shí)間 X2、獲得資助的指標(biāo) X3之間的散點(diǎn)圖。從圖中可以看出這些點(diǎn)大致分布在一條直線旁邊，因此，有比較好的線性關(guān)系，可以采用線性回歸。繪制圖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與各自變量的樣本散點(diǎn)圖

（2）進(jìn)行多元線性回歸

這里可以直接使用 regress 函數(shù)執(zhí)行多元線性回歸，具體代碼如下：

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

運(yùn)行后即得到結(jié)果如表4所示。

表4對(duì)初步回歸模型的計(jì)算結(jié)果

計(jì)算結(jié)果包括回歸系數(shù) b = (β0,β1,β2,β3) = (18.0157, 1.0817, 0.3212, 1.2835)、回歸系數(shù)的置信區(qū)間，以及統(tǒng)計(jì)變量 stats（它包含四個(gè)檢驗(yàn)統(tǒng)計(jì)量：相關(guān)系數(shù)的平方R^2，假設(shè)檢驗(yàn)統(tǒng)計(jì)量 Ｆ，與 F 對(duì)應(yīng)的概率 p，s^2的值）。因此我們得到初步的回歸方程為：

由結(jié)果對(duì)模型的判斷：

回歸系數(shù)置信區(qū)間不包含零點(diǎn)表示模型較好，殘差在零點(diǎn)附近也表示模型較好，接著就是利用檢驗(yàn)統(tǒng)計(jì)量 Ｒ，Ｆ，p 的值判斷該模型是否可用。

1）相關(guān)系數(shù) Ｒ 的評(píng)價(jià)：本例 Ｒ 的絕對(duì)值為 0.9542 ，表明線性相關(guān)性較強(qiáng)。

2）F 檢驗(yàn)法：當(dāng) F > F1-α(m,n-m-1) ，即認(rèn)為因變量 y 與自變量x1,x2,...,xm之間有顯著的線性相關(guān)關(guān)系；否則認(rèn)為因變量 y 與自變量x1,x2,...,xm之間線性相關(guān)關(guān)系不顯著。本例 Ｆ＝67.919 > F1-0.05( 3,20 ) = 3.10。

3）p 值檢驗(yàn)：若 p < α（α 為預(yù)定顯著水平），則說明因變量 y 與自變量?x1,x2,...,xm之間顯著地有線性相關(guān)關(guān)系。本例輸出結(jié)果，p<0.0001,顯然滿足 p<α=0.05。

以上三種統(tǒng)計(jì)推斷方法推斷的結(jié)果是一致的，說明因變量 y與自變量之間顯著地有線性相關(guān)關(guān)系，所得線性回歸模型可用。s^2當(dāng)然越小越好，這主要在模型改進(jìn)時(shí)作為參考。

3.逐步歸回

[ 例4 ]（Hald,1960）Hald 數(shù)據(jù)是關(guān)于水泥生產(chǎn)的數(shù)據(jù)。某種水泥在凝固時(shí)放出的熱量 Y（單位：卡/克）與水泥中 4 種化學(xué)成品所占的百分比有關(guān)：

在生產(chǎn)中測(cè)得 12 組數(shù)據(jù)，見表5，試建立 Y關(guān)于這些因子的“最優(yōu)”回歸方程。

表5水泥生產(chǎn)的數(shù)據(jù)

對(duì)于例 4 中的問題，可以使用多元線性回歸、多元多項(xiàng)式回歸，但也可以考慮使用逐步回歸。從逐步回歸的原理來看，逐步回歸是以上兩種回歸方法的結(jié)合，可以自動(dòng)使得方程的因子設(shè)置最合理。對(duì)于該問題，逐步回歸的代碼如下：

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];%自變量數(shù)據(jù)

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];%因變量數(shù)據(jù)

Stepwise（X,Y,[1,2,3,4],0.05,0.10）% in=[1,2,3,4]表示X1、X2、X3、X4均保留在模型中

程序執(zhí)行后得到下列逐步回歸的窗口，如圖 4 所示。

圖4逐步回歸操作界面

在圖 4 中，用藍(lán)色行顯示變量 X1、X2、X3、X4 均保留在模型中，窗口的右側(cè)按鈕上方提示：將變量X3剔除回歸方程（Move X3 out），單擊 Next Step 按鈕，即進(jìn)行下一步運(yùn)算，將第 3 列數(shù)據(jù)對(duì)應(yīng)的變量 X3 剔除回歸方程。單擊 Next Step 按鈕后，剔除的變量 X3 所對(duì)應(yīng)的行用紅色表示，同時(shí)又得到提示：將變量 X4 剔除回歸方程（Move X4 out），單擊 Next Step 按鈕，這樣一直重復(fù)操作，直到 “Next Step” 按鈕變灰，表明逐步回歸結(jié)束，此時(shí)得到的模型即為逐步回歸最終的結(jié)果。

4.Logistic 回歸

[ 例5 ]企業(yè)到金融商業(yè)機(jī)構(gòu)貸款，金融商業(yè)機(jī)構(gòu)需要對(duì)企業(yè)進(jìn)行評(píng)估。評(píng)估結(jié)果為 0 , 1 兩種形式，0 表示企業(yè)兩年后破產(chǎn)，將拒絕貸款，而 1 表示企業(yè) 2 年后具備還款能力，可以貸款。在表 6 中，已知前 20 家企業(yè)的三項(xiàng)評(píng)價(jià)指標(biāo)值和評(píng)估結(jié)果，試建立模型對(duì)其他 5 家企業(yè)（企業(yè) 21-25）進(jìn)行評(píng)估。

表6企業(yè)還款能力評(píng)價(jià)表

對(duì)于該問題，很明顯可以用 Logistic 模型來回歸，具體求解程序如下：

% logistic回歸Matlab實(shí)現(xiàn)程序

%% 數(shù)據(jù)準(zhǔn)備

clc, clear, closeall

X0=xlsread('logistic_ex1.xlsx', 'A2:C21');% 回歸模型的輸入

Y0=xlsread('logistic_ex1.xlsx', 'D2:D21');% 回歸模型的輸出

X1=xlsread('logistic_ex1.xlsx', 'A2:C26');% 預(yù)測(cè)數(shù)據(jù)輸入

%% logistics函數(shù)

GM = fitglm(X0,Y0,'Distribution','binomial');

Y1 = predict(GM,X1);

%% 模型的評(píng)估

N0 =1:size(Y0,1); N1= 1:size(Y1,1);

plot(N0', Y0, '-kd');

hold on; scatter(N1', Y1, 'b')

xlabel('數(shù)據(jù)點(diǎn)編號(hào)'); ylabel('輸出值');

得到的回歸結(jié)果與原始數(shù)據(jù)的比較如圖5所示。

圖5 回歸結(jié)果與原始數(shù)據(jù)的比較圖

5.小結(jié)

本講主要介紹數(shù)學(xué)建模中常用的幾種回歸方法。在使用回歸方法的時(shí)候，首先可以判斷自變量的個(gè)數(shù)，如果超過 2 個(gè)，則需要用到多元回歸的方法，否則考慮用一元回歸。然后判斷是線性還是非線性，這對(duì)于一元回歸是比較容易的，而對(duì)于多元，往往是將其他變量保持不變，將多元轉(zhuǎn)化為一元再去判斷是線性還是非線性。如果變量很多，而且復(fù)雜，則可以首先考慮多元線性回歸，檢驗(yàn)回歸效果，也可以用逐步回歸?？傊?，用回歸方法比較靈活，根據(jù)具體情景還是比較容易找到合適的方法的。