實(shí)驗(yàn)三-非線性回歸分析

1、實(shí)驗(yàn)三 非線性回歸分析（2學(xué)時(shí)）一、實(shí)驗(yàn)?zāi)康模?）、掌握非線性回歸分析的基本步驟；（2）、熟練用MATLAB軟件實(shí)現(xiàn)可線性化的回歸分析及曲線回歸分析。 二、實(shí)驗(yàn)學(xué)時(shí)：2學(xué)時(shí)三、實(shí)驗(yàn)要求（1）掌握用MATLAB軟件實(shí)現(xiàn)可線性化的回歸分析及曲線回歸分析；（2）掌握非線性回歸分析的基本步驟。 四、實(shí)驗(yàn)原理 1、常見的可線性化模型 (1)雙對(duì)數(shù)模型 (2).半對(duì)數(shù)模型 對(duì)數(shù)-線性模型: 線性-對(duì)數(shù)模型: (3). 倒數(shù)模型 (4) 邏輯(logistic)成長(zhǎng)曲線模型 （5）龔伯斯（Gompertz）成長(zhǎng)曲線模型 (6) 多項(xiàng)式模型 一元二次多項(xiàng)式模型： 一元三次多項(xiàng)式模型： 二元二次多項(xiàng)式模型：2、

2、不可線性化模型（即曲線回歸模型）： 其中f為非線性的，并且模型不可線性化。 一般使用非線性最小二乘估計(jì)方法，并用Newton迭代求解其中的正規(guī)方程組。五、實(shí)驗(yàn)舉例例1、對(duì)GDP（國(guó)內(nèi)生產(chǎn)總值）的擬合。選取GDP指標(biāo)為因變量，單位為億元，擬合GDP關(guān)于時(shí)間t的趨勢(shì)曲線。以1981年為基準(zhǔn)年，取值為t=1,1998年t=18,1991-1998年的數(shù)據(jù)如下： 年份t GDP 年份tGDP14862.41018547.925294.71121617.835934.51226638.1471711334634.458964.41446759.4610202.21558478.1711962.51667

3、884.6814928.31774462.6916909.21879395.7解：第一步：一元線性回歸模型（一）程序如下：建立test3_1_1.m文件如下（具體注釋見m文件）：% 作時(shí)間t 與GDP（y）的回歸分析data,head=xlsread('test3_1.xlsx'); % 導(dǎo)入數(shù)據(jù)t=data(:,1); % 提取年份ty=data(:,2); % 國(guó)內(nèi)生產(chǎn)總值plot(t,y,'k.'); % 畫y與t的散點(diǎn)圖xlabel('年序號(hào)（t）')ylabel('GDP（y）')% 調(diào)用robustfit函數(shù)作穩(wěn)健回歸

4、，返回系數(shù)的估計(jì)值b和相關(guān)統(tǒng)計(jì)量stats b,stats=robustfit(t,y)stats.p % 繪制殘差與權(quán)重的散點(diǎn)圖 %figure;plot(stats.resid,stats.w,'o') % 繪制殘差與權(quán)重的散點(diǎn)圖xlabel('殘差')ylabel('權(quán)重') % 畫robustfit函數(shù)對(duì)應(yīng)的回歸直線 %tdata=ones(size(t,1),1),t; % 在原始數(shù)據(jù)t的左邊加一列1yhat=tdata*b; % 求robustfit 函數(shù)對(duì)應(yīng)的y的估計(jì)值figure;plot(t,y,'ko') %

5、畫原始數(shù)據(jù)散點(diǎn)hold onplot(t,yhat,'r-','linewidth',2) % 畫robustfit函數(shù)對(duì)應(yīng)的回歸直線，紅色虛線xlabel('年序號(hào)（t）')ylabel('GDP（y）')（二）實(shí)驗(yàn)結(jié)果與分析：（1）散點(diǎn)圖圖3.1 年序號(hào)t與GDP(y)的散點(diǎn)圖從散點(diǎn)圖可以看出t與y的線性趨勢(shì)并不明顯，但可以考慮進(jìn)行一元線性回歸。（2）調(diào)用robustfit函數(shù)作穩(wěn)健回歸，返回系數(shù)的估計(jì)值b和相關(guān)統(tǒng)計(jì)量stats的P值 系數(shù)的估計(jì)值 P值1.0e+04 * -1.3410 0.0227 0.4414 0.000

6、0 穩(wěn)健回歸得出的回歸方程為 y=-13410+4414*t 。常數(shù)項(xiàng)和回歸系數(shù)的t檢驗(yàn)的p值只有一個(gè)小于顯著性水平0.0001，可知線性關(guān)系是不顯著的。 （3）下面作出殘差向量和權(quán)重向量的散點(diǎn)圖，從直觀上了解殘差和權(quán)重的關(guān)系。圖3.2 殘差和權(quán)重的散點(diǎn)圖 從圖3.2可以看出殘差絕對(duì)值越大，其權(quán)重就越小，殘差為零時(shí)，相應(yīng)的權(quán)重為1，這就保證了擬合的穩(wěn)健性。此處殘差的數(shù)量級(jí)為104，從而殘差絕對(duì)值太大，可知明顯擬合效果不好。下面作出擬合效果圖，如圖3.3所示。 圖3.3 原始數(shù)據(jù)散點(diǎn)與回歸直線圖 從上圖可知：擬合效果明顯很不好。原始數(shù)據(jù)點(diǎn)偏離回歸直線距離太大。下面進(jìn)行復(fù)合函數(shù)模型的擬合，并進(jìn)行擬

7、合效果比較。第二步：作復(fù)合函數(shù)模型 (1)散點(diǎn)圖：圖3.4 年序號(hào)t與GDP(lny)的散點(diǎn)圖從圖3.4的散點(diǎn)圖表明年序號(hào)t與lny的線性趨勢(shì)比圖3.1的線性趨勢(shì)明顯。（2）調(diào)用robustfit函數(shù)作穩(wěn)健回歸，返回系數(shù)的估計(jì)值b和相關(guān)統(tǒng)計(jì)量stats的P值'系數(shù)的估計(jì)值' '檢驗(yàn)的P值' 8.1896 1.0e-16 * 0.1756 0.00000.9250穩(wěn)健回歸得出的回歸方程為 lny=8.1896+0.1756*t。常數(shù)項(xiàng)和回歸系數(shù)的t檢驗(yàn)的p值分別為0，0.925×10-16，均小于顯著性水平0.0001，可知線性關(guān)系是極顯著的。（2）繪制

8、殘差與權(quán)重的散點(diǎn)圖圖3.5 殘差與權(quán)重的散點(diǎn)圖與圖3.2相比較，從圖3.5可知?dú)埐罱^對(duì)值顯著減小，殘差的數(shù)量級(jí)為10-1，并且權(quán)重集中在1附近的點(diǎn)比之前多。這說(shuō)明擬合的穩(wěn)健性顯著加強(qiáng)。（3）回歸直線與總結(jié)圖3.6 原始數(shù)據(jù)散點(diǎn)與回歸直線圖（lny） 從上圖可知：此次擬合效果是顯著的?？偨Y(jié)：線性回歸方程： y=-13410+4414*t p>0.0001 復(fù)合函數(shù)回歸方程： y=3603.3*(1.192)t即 lny=8.1896+0.1756*t p<0.0001所以由t檢驗(yàn)的p值，再結(jié)合回歸直線圖知：復(fù)合函數(shù)模型比線性回歸模型要好。 例2 、一位藥物學(xué)家是用下面的非線性模型對(duì)藥

9、物反應(yīng)擬合回歸模型 其中，自變量x為藥劑量，用級(jí)別表示； 因變量y為藥物反應(yīng)程度，用百分?jǐn)?shù)表示。三個(gè)參數(shù)c0 ,c1 ,c2都是非負(fù)的, c0 的上限是100%，三個(gè)參數(shù)的初始值取為c0 =100,c1=5 ,c2=4.8.測(cè)得9個(gè)數(shù)據(jù)如下表：x123456789y(%)0.52.33.42454.782.194.896.296.4解：（一）程序如下：x=1:9;y=0.5 2.3 3.4 24 54.7 82.1 94.8 96.2 96.4;plot(x,y,'k.') xlabel('藥劑量（x）') ylabel('藥物反應(yīng)程度（y）')

10、 options=statset;options.Robust='on' % 調(diào)用nlinfit 函數(shù)作非線性回歸beta,r,J,COVB,mse=nlinfit(x,y,fun1,100,5,4.8,options);beta % 查看未知參數(shù)估計(jì)值 mse % 查看均方殘差平方和（誤差方差2的估計(jì)）mse% 求參數(shù)估計(jì)值的95%置信區(qū)間ci1=nlparci(beta,r,'covar',COVB,'alpha',0.05)% 繪制robustfit 函數(shù)對(duì)應(yīng)的回歸直線plot(x,yhat,'linewidth',3)xl

11、abel('藥劑量（x）')ylabel('藥物反應(yīng)程度（y）')legend('原始數(shù)據(jù)散點(diǎn)','非線性回歸曲線','location','northwest')（二）實(shí)驗(yàn)結(jié)果與分析： （1）散點(diǎn)圖圖3.7 藥劑量（x）與藥物反應(yīng)程度(y)的散點(diǎn)圖（2）回歸模型與誤差方差2的估計(jì)參數(shù)C0 C1 C2 誤差方差2的估計(jì) 99.5054 6.7701 4.7987 3.8895 由上述可得非線性回歸模型為 y=99.51-99.51/(1+(x4.80)6.77。誤差方差2的估計(jì)mse=3.8895，

12、這個(gè)數(shù)還是非常小的，擬合效果還是很好的。 （3）非線性回歸曲線圖3.8 原始數(shù)據(jù)散點(diǎn)與非線性回歸曲線從上圖可以看出擬合效果還是很不錯(cuò)的。但非線性回歸，不能像線性回歸那樣對(duì)回歸方程作顯著性檢驗(yàn)。所以，非線性回歸模型為 y=99.51-99.51/(1+(x4.80)6.77對(duì)該藥物反應(yīng)擬合回歸模型是合理的。六、實(shí)驗(yàn)內(nèi)容 Logistic 回歸函數(shù)常用于擬合某種消費(fèi)品的擁有率，下表是北京市每百戶家庭平均擁有的照相機(jī)數(shù)，試針對(duì)以下兩種擬合Logistic回歸函數(shù)：(1) 已知u=100，用線性化方法擬合； (2) u未知，用非線性最小化方法擬合。年份ty年份ty197817.519881159.6197929.819891262.21980311.419901366.51981413.319911472.71982517.219921577.21983620.619931682.41984729.119941785.41985834.619951886.81986947.419961987.219871055.5 七、思考練習(xí)某種商品的流通率y(%