SAS學(xué)習(xí)系列25.非線性回歸

25.非線性回歸現(xiàn)實(shí)世界中嚴(yán)格的線性模型并不多見，它們或多或少都帶有某種程度的近似；在不少情況下，非線性模型可能更加符合實(shí)際。對(duì)變量間非線性相關(guān)問題的曲線擬合，處理的方法主要有：（1）首先確定非線性模型的函數(shù)類型，對(duì)于其中可線性化問題則通過變量變換將其線性化，從而歸結(jié)為前面的多元線性回歸問題來解決；（2）若實(shí)際問題的曲線類型不易確定時(shí)，由于任意曲線皆可由多項(xiàng)式來逼近，故?？捎枚囗?xiàng)式回歸來擬合曲線；（3）若變量間非線性關(guān)系式已知（多數(shù)未知），且難以用變量變換法將其線性化，則進(jìn)行數(shù)值迭代的非線性回歸分析。一）可變換為線性的非線性回歸表1典型的函數(shù)及線性化方法"函數(shù)君稱.國(guó)數(shù)盍達(dá)式？線性化方法.取曲線函歎J1b—=住一一.■yxv=—u=—VX冠函數(shù)J尹二OY?V=111J11-111X指數(shù)函數(shù).y=*v=Invu=xry=aebx^v=Inv起二丄型l，X對(duì)數(shù)函數(shù).y-a-i-31nx■-v-yu=InW型函數(shù)十二+"v=—u=e~x盧在很多場(chǎng)合，可以對(duì)非線性模型進(jìn)行線性化處理，尤其是可變換為線性的非線性回歸，運(yùn)用最小二乘法進(jìn)行推斷，對(duì)線性化后的線性模型，可以應(yīng)用REG過程步進(jìn)行計(jì)算。例1有實(shí)驗(yàn)數(shù)據(jù)如下：X1.11.21.31.41.51.61.71.81.922.12.22.32.4109.940.420.024.511.07.34.92.71.81.40.80.0.0.2Y55932952292322試分別采用指數(shù)回歸（y二aebx）方法進(jìn)行回歸分析。代碼：dataexam25_l;inputxy;cards;109.95TOC\o"1-5"\h\z40.4520.0924.5311.027.391.74.951.82.721.91.8221.492.10.822.20.32.30.22.40.22run;procsgplotdata=exam25_l;scatterx=xy=y;run;proccorrdata=exam25_1;varxy;run;datanew1;setexam25_1;v=log(y);run;procsgplotdata=new1;scatterx=xy=v;title'變量代換后數(shù)據(jù)’;run;procregdata=newl;varxv;modelv=x;printcli;title'殘差圖’；plotresidual.*predicted.;run;datanew2;setexam25_l;yl=14530.28*exp(-4.73895*x);run;procgplotdata=new2;ploty*x=1yl*x=2/overlay;symbolv=doti=nonecv=red;symbol2i=smcolor=blue;title'指數(shù)回歸圖’；run;運(yùn)行結(jié)果：CORK過程2變量：J卜2變量：I￥簡(jiǎn)單統(tǒng)計(jì)量變量1-1均值標(biāo)準(zhǔn)差總和最小值最大值1.H4XK)0.4183324.500001.KJO(Jt)24000016.1392929.4S037225.950000.2000010fl9500DPearson相關(guān)系數(shù)，N-14Proby|r|under110:Rho-01.00000■0.71104【).■0.711041.000000.￡)014讀取的觀測(cè)數(shù)H使用的觀測(cè)數(shù)14方差分柄自由度平方和均方F值Pr>F模型151.0910651.0910S750.31<0001誤差1?D.810630.06755校正合計(jì)13b1.901心均方根誤差0.25991K方09344因?qū)D量均值1.29U?+調(diào)整R方(1皈丨變異系數(shù)20,13438盞數(shù)估計(jì)值自由度參數(shù)估計(jì)值準(zhǔn)差標(biāo)誤t值Pr>|t|1忤1匕「cepL19,583990.30945<.0001X1-4.J3395().17232-27.50<.0001

程序說明：(1)調(diào)整后的R2=0.9831，說明擬合程度很好；F檢驗(yàn)的P值=0.0001<a=0.05，拒絕原假設(shè)，故直線回歸的斜率不為0；將線性回歸系數(shù)代入，得到原回歸方程y=14530.28*e4.73895x殘差圖趨勢(shì)，符合殘差隨機(jī)正態(tài)分布的假設(shè)(不帶其它明顯趨勢(shì))。二、多項(xiàng)式回歸一般函數(shù)都可用多項(xiàng)式來逼近，故多項(xiàng)式回歸分析可用來處理相當(dāng)廣泛的非線性問題。對(duì)觀測(cè)數(shù)據(jù)（xt,yt）,t=1,…，N.多項(xiàng)式回歸模型為：yt=+0［兀+仇工;+…++8,,Fl,2,令y=,x=X]—T?—'B=1A-力i?s=左2，1-ryxn"■xjV_電」則模型可寫為：Y=XB+e當(dāng)X列滿秩時(shí)，用最小二乘估計(jì)B=（XX丄1XY可求得其多項(xiàng)式回歸方程。但由于（XX）-1的計(jì)算既復(fù)雜又不穩(wěn)定，故一般采用正交多項(xiàng)式法來進(jìn)行多項(xiàng)式回歸。多項(xiàng)式模型可以直接應(yīng)用GLM（廣義線性模型）求解。例2市種畜場(chǎng)奶牛群1—12月份（X1），產(chǎn)犢母牛平均產(chǎn)奶量（y）的資料如程序數(shù)據(jù)步中，試對(duì)該資料配置一個(gè)合適的回歸方程。代碼：dataexam25_2;inputx1y;x2=x1*x1;datalines;3833.433476.763811.583466.223769.473395.423565.743807.083481.993817.033372.823884.52run;procsgplotdata=exam25_2;scatterx=x1y=y;title'原始數(shù)據(jù)散點(diǎn)圖’；run;procregdata=exam25_2;modely=x1x2;run;運(yùn)行結(jié)果：原始數(shù)據(jù)故點(diǎn)圖MDMDa1REG過稈模型：a1REG過稈模型：MODEL1

國(guó)變量：y讀取的觀測(cè)數(shù)I2悚用的觀測(cè)數(shù)I2山差分新]源自由度平方和均方F值Pr>1模型23325941662S71U.890.0009誤差J)886019644.56254校正合計(jì)11421195均方根送差99.21977R方0.7896因?qū)D量均值3540.17167調(diào)整F?方0.7429變異系數(shù)2.72569哄下涮象的殘姜與回歸凰：yAlM過程源白由度平方和r均方卜值hr>卜模型2332594.3259166297.163016.890.0009988601.06299844.5625校正合計(jì)11421195.3888R方變異系數(shù)根MSEy均值0.78S54499.219773640.172源自由度1型SS均方1值Pr>1x111I.仃幀Ii.nosnn.on0.9741x1*x11332583.3210332533.321033.780.0003源自由度111型SS均方r值Pr>rx1I314325.5414314325.541431.S30.0003x1*x11332583.3210332583.32100.0003參數(shù)估計(jì)值標(biāo)準(zhǔn)誤差t值Pr>|t|1nteroept4117.201364102.545490540.15<.0001x1-204,93067836,2683725-!).650.[]0(Hx1*x115.7856992.7158(3775.810.0003程序說明：（1）觀察數(shù)據(jù)的散點(diǎn)圖，更適合二次多項(xiàng)式擬合，也可以測(cè)試幾種不同次數(shù)的多項(xiàng)式擬合選擇其中最優(yōu)的；（2）將回歸系數(shù)代入多項(xiàng)式方程得到：y=4117.20136-204.93668x1+15.78570x12三、不能變換為線性的非線性回歸該類非線性回歸分析就是利用最小二乘準(zhǔn)則來估計(jì)回歸系數(shù)〃，使得殘差平方和最小。一般來用數(shù)值迭代法來進(jìn)行，先選定回歸系數(shù)的初值/30，按照給定的步長(zhǎng)和搜索方向逐步迭代，直到殘差平方和達(dá)到最小。有5種常用的非線性回歸迭代方法：高斯-牛頓法（Gauss-Newton）、最速下降法（梯度法）、牛頓法（Newton）、麥夸特法（Marquardt）、正割法（DUD）。高斯-牛頓法在初值選取適當(dāng)，且可逆時(shí)非常有效，但在其他情形，其求解較為困難，對(duì)此，Marguardt對(duì)其中的正則系數(shù)陣作適當(dāng)修正，得到了改進(jìn)算法。（二）PROCNLIN過程步對(duì)于不能線性化的非線性模型。其估計(jì)不能直接運(yùn)用經(jīng)典的最小二乘法，而需要運(yùn)用其他估計(jì)方法，如加權(quán)最小二乘法、直接搜索法、直接最優(yōu)法與Taylor級(jí)數(shù)展開法進(jìn)行線性逼近此時(shí)可以利用NLIN過程步實(shí)現(xiàn)相應(yīng)的計(jì)算，它是采用最小誤差平方法及迭代推測(cè)法來建立一個(gè)非線性模型，估計(jì)參數(shù)默認(rèn)采用高斯-牛頓迭代法。NLIN過程不保證一定可以算出符合最小誤差平方法之標(biāo)準(zhǔn)的參數(shù)估計(jì)值?；菊Z(yǔ)法：PROCNLINdata=數(shù)據(jù)集</可選項(xiàng)>;PARMS參數(shù)名二數(shù)值；MODEL因變量二表達(dá)式</可選項(xiàng)〉；BOUNDS不等式；〉CONTROL變量；〉der.參數(shù)名二表達(dá)式;〉<OUTPUTout=輸出數(shù)據(jù)集</可選統(tǒng)計(jì)量>；>說明：（1）NLIN的可選項(xiàng)包括：outest=輸出數(shù)據(jù)集輸出每步迭代的結(jié)果；best=n只輸出最好的n組殘差平方和；method=gauss|marquardt|newton|gradient|dud|設(shè)定參數(shù)估計(jì)的迭代方法，默認(rèn)為gauss(沒有der.語(yǔ)句)；PARMS語(yǔ)句指定參數(shù)并賦值，一般包括參數(shù)名、初始值(GridSearch可以幫助選擇合適的初始值)、迭代準(zhǔn)則；例如：parmsb0=0b1=1to10b2=1to10by2b3=1,10,100;bounds語(yǔ)句用于設(shè)定參數(shù)的約束，主要是不等式約束，約束間用逗號(hào)分隔。例如，boundsa<=20,b>30,1<=c<=10;der.語(yǔ)句用于計(jì)算模型關(guān)于各參數(shù)的偏導(dǎo)數(shù)，相應(yīng)格式為：一階偏導(dǎo)數(shù)：der?參數(shù)名=表達(dá)式；二階偏導(dǎo)數(shù)：der?參數(shù)名?參數(shù)名=表達(dá)式；例如，對(duì)于模型modely=b0*(l—exp(-bl*x));二階偏導(dǎo)數(shù)表達(dá)式：der.b0.bl二x*exp(-b1*x);例3根據(jù)對(duì)已有數(shù)據(jù)的XY散點(diǎn)圖的觀察和分析，發(fā)現(xiàn)Y隨X增長(zhǎng)趨勢(shì)是減緩的，并且Y趨向一個(gè)極限值，我們認(rèn)為用負(fù)指數(shù)增長(zhǎng)曲線來擬合模型較為合適。代碼：dataexpd;inputxy;datalines;0200.570300.720400.810500.870600.910700.940800.950900.971000.981100.991201.001300.991400.991501.001601.001700.991801.001901.002000.992101.00=10method=gauss;procnlin=10method=gauss;parmsb0=0to2by0.5b1=0.01to0.09by0.01;modely=b0*(1-exp(-bl*x));der.b0=1-exp(-bl*x);der.bl=b0*x*exp(-bl*x);outputout=expoutp=ygs;run;goptionsreset=globalgunit=pctcback=whiteborderhtitle=6htext=3ftext=swissbcolors=(back);procgplotdata=expout;ploty*xygs*x/haxis=axislvaxis=axis2overlay;symbolli=nonev=pluscv=redh=2.5w=2;symbol2i=joinv=nonel=1h=2.5w=2;axis1order=20to210by10;axis2order=0?5to1.1by0.05;title1'y=b0*(1-exp(-b1*x)';NOTENOTE:Convergenceoriterionnettitle2'procnlinmethod=gauss';run;運(yùn)行結(jié)果：NLIN過程

因變量T網(wǎng)格搜索bObl平方和i.（）（）［）（）().(JO11()L00000.05000.0I6S1.00000.06000.05521.（）（!［）（）0,03000,0666100000.07000.00731.00000.ODDO0.1365i.o(|[)o0,09000,1708100000.02000.41931.50000.01000.9757L00000.01002,1653NLIN過程因豐量y

方法：Gauss~N<mI.on迭代階段迭代bOb1平方和01.uouo0.04000.0U14U10.99010.04190.0005802Q.畫20.04200.OOO&7730.99520.04200.000&7T4t）.腳汐0.0-^(10.OOO&77估計(jì)匯總方法GaussNewton迭代IR2.552E-7PPG(bl)1.()281-8RIJC(b1)匚抽-1「7Object2.56E-10目標(biāo)0.000677讀取的觀測(cè)20愉用的觀測(cè)20缺失的觀測(cè)UNuLe:Ariinteri；epLwuurioLifiedlorIhisniudeI.源自由度平方和均方rva1ue近似Pr>r模型217.67176.9359275733<.0001渓差1U0.0005770.000032未校正合計(jì)2017.6723參數(shù)怙計(jì)值近似標(biāo)準(zhǔn)淒差近似95%査信限bO0.99