試驗設(shè)計與數(shù)據(jù)處理第六章續(xù)(第七,八節(jié))

第七節(jié)二次響應(yīng)面分析法當試驗中考察的指標宜于用多元二次回歸方程來擬合因素與指標的函數(shù)關(guān)系，就可以分析回歸方程所反應(yīng)的曲面形狀，如果得到的曲面是凸面(像山丘)或凹面(像山谷)這類簡單曲面，那么預(yù)測的最佳指標值(極大值或極小值)可以從所估計的曲面上獲得；如果曲面很復(fù)雜，或者預(yù)測的最佳點遠離所考察因素的試驗范圍，那么可以通過嶺嵴分析來確定重新進行試驗的方向.這就是應(yīng)用較廣，頗有實用價值的響應(yīng)面分析法(ResponseSurfaceAnalysis).響應(yīng)面分析的三個任務(wù):任務(wù)1.模型擬合和方差分析;任務(wù)2.典型分析(尋找最優(yōu)值點);任務(wù)3.嶺嵴分析(若最優(yōu)值點在試驗范圍內(nèi)達到,則不必進行嶺嵴分析;否則就要進行).例6.7.2求面包酵母高產(chǎn)和高發(fā)酵活力的最優(yōu)工藝參數(shù).指標:y1---面包酵母產(chǎn)量[指最終發(fā)酵液的菌體濃度(g/L)]y2---發(fā)酵活力[指面團從水中浮起的時間(min)]因素:z1---還原糖濃度,2%~6%;z2---氮源濃度,0.4%~1.2%;z3---磷源濃度,0.2%~0.4%.取三水平,z1:2%,4%,6%;z2:0.4%,0.8%,1.2%;z3:0.2%,0.3%,0.4%.令編碼變量xi=[zi-(zi1+zi2)/2]/(zi2-zi1)/2故可作變換:X1=(z1-4)/2X2=(z2-0.8)/0.4X3=(z3-0.3)/0.1響應(yīng)面分析的SAS程序如下:DataE62;Inputz1-z3x1-x3y1y2@@;Cards;

40.40.20-1-124.513.9

40.40.40-1135.314.1………40.80.300037.58.1

;Proc

RsReg

data=E62;/*響應(yīng)面分析*/Modely1y2=z1-z3;Run;

TheRSREGProcedureCodingCoefficientsfortheIndependentVariablesFactorSubtractedoffDividedbyz14.0000002.000000z20.8000000.400000z30.3000000.100000ResponseSurfaceforVariabley1ResponseMean32.400000RootMSE1.329474R-Square0.9775CoefficientofVariation4.1033下面給出程序運行得到的部分結(jié)果:

TypeISumRegressionDFofSquaresR-SquareFValuePr>FLinear314.5575000.03712.750.1524Quadratic3116.5325000.297221.980.0026Crossproduct3252.1725000.643147.560.0004TotalModel9383.2625000.977524.090.0013SumofResidualDFSquaresMeanSquareTotalError58.8375001.767500

TheRSREGProcedureCanonicalAnalysisofResponseSurfaceBasedonCodedDataCriticalValueFactorCodedUncodedz1-0.0988443.802311z20.2648400.905936z3-0.3158000.268420Predictedvalueatstationarypoint:38.040882EigenvectorsEigenvaluesz1z2z3-0.5185070.184252-0.7614870.621441-1.7733340.894929-0.131446-0.426408-8.0206590.4063900.6347120.657258Stationarypointisamaximum.如果穩(wěn)定點不是理想點就要進一步作嶺嵴分析,請看示意圖和例子演示(例見datarubber或datapectic)例6.7.3橡膠樹幼苗對硫35的吸收實驗(2×3析因試驗)因素:溫度t(℃)---7℃,20℃

取樣時間d---2d,7d,12d響應(yīng)面分析SAS程序如下:datarubber;inputftd@@;cards;數(shù)(略);procsort;bytd;/*對自變量td進行sort由小到大排序*/procrsreg;modelf=td/lackfit;/*選項lackfit要求對有重復(fù)試驗自由度的回歸模型執(zhí)行不適合度檢定(lack-of-fittest)(誤差參照為重復(fù)試驗的均方誤差),預(yù)先應(yīng)先對自變量進行sort由小到大排序*/run;可以作二個因素的響應(yīng)面圖(固定其它因素),E62的響應(yīng)面圖如下:回歸方程作響應(yīng)面,在對數(shù)據(jù)集E62進行響應(yīng)面回歸(Rsreg)后得回歸方程(固定水平x3=0,本例自變量均已經(jīng)編碼處理)y1=37.9-0.5625x1-0.6x2-2.875x12-3.575x1x2-3.45x22SAS程序如下:

DataE62G;Dox1=-1to1by0.1;/*當自變量是原始變量時取類似的二元網(wǎng)格點,即do變量=左端點

to右端點

by步長*/Dox2=-1to1by0.1;/*當自變量是原始變量時取類似的二元網(wǎng)格點,即do變量=左端點

to右端點

by步長*/Y1=37.9-0.5625*x1-0.6*x2-2.875*x1**2-3.575*x1*x2-3.45*x2**2;Output;End;End;Procprint;Run;/*以上為由回歸方程建立作圖數(shù)據(jù)集*/等值線圖(contourplot)37.95/*以下利用程序作三維圖*/procg3ddata=E62G;plotX2*X1=Y1;/*作x1(縱)-x2(橫)-y1(豎)三維圖/run;第八節(jié)非線性回歸非線性回歸是指模型從參數(shù)角度為非線性的且無法線性化的模型,例如y=b0+b1x1+b2x2+b11x12+b12x1x2+b22x22是諸參數(shù)bi及bij的線性模型(從自變量x1,x2角度模型是非線性的),y=aebx是可以線性化的非線性模型(線性化后lny=bx+lna),而y=aexp(-exp(b-cx))是非線性的且無法線性化的模型.此時SAS采用nlin非線性回歸程序程序形式:procnlindata=…;model因變量

=非線性函數(shù)表達式;

parms

參數(shù)1=初值參數(shù)2=初值

…;

der.參數(shù)1=偏導(dǎo)表達式;der.參數(shù)2=偏導(dǎo)表達式;……………;runs;方法:采用迭代求參數(shù),要給出參數(shù)的一階偏導(dǎo)表達式der.參數(shù)=…;,要賦參數(shù)初值parms參數(shù)=…;請看例子演示例6.8.1已知牧草重量y與生長天數(shù)x的關(guān)系是y=aexp(-exp(b-cx))+ε.9次觀察的數(shù)據(jù)如表6.8.1,求a,b,c的值.datahw;inputxy;cards;數(shù)據(jù)(略)

；procnlin

data=hw;/*調(diào)用nlin進行非線性回歸*/zz=exp(b-c*x);/*簡化記號*/term=exp(-zz);/*簡化記號*/model

y=a*term;

/*因變量=非線性函數(shù)表達式,即y=aexp(-exp(b-cx))*/parmsa=70b=1.1c=0.2;/*賦參數(shù)初值*/der.a=term;/*給出參數(shù)的一階偏導(dǎo)表達式*/der.b=-a*zz*term;der.c=a*x*zz*term;run;

DependentVariableyMethod:Gauss-NewtonSumofIterbcaSquares01.10000.200070.00009161.21-0.1190