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

第七節(jié)二次響應(yīng)面分析法當(dāng)試驗(yàn)中考察的指標(biāo)宜于用多元二次回歸方程來擬合因素與指標(biāo)的函數(shù)關(guān)系，就可以分析回歸方程所反應(yīng)的曲面形狀，如果得到的曲面是凸面(像山丘)或凹面(像山谷)這類簡單曲面，那么預(yù)測的最佳指標(biāo)值(極大值或極小值)可以從所估計(jì)的曲面上獲得；如果曲面很復(fù)雜，或者預(yù)測的最佳點(diǎn)遠(yuǎn)離所考察因素的試驗(yàn)范圍，那么可以通過嶺嵴分析來確定重新進(jìn)行試驗(yàn)的方向.這就是應(yīng)用較廣，頗有實(shí)用價(jià)值的響應(yīng)面分析法(ResponseSurfaceAnalysis).響應(yīng)面分析的三個(gè)任務(wù):任務(wù)1.模型擬合和方差分析;任務(wù)2.典型分析(尋找最優(yōu)值點(diǎn));任務(wù)3.嶺嵴分析(若最優(yōu)值點(diǎn)在試驗(yàn)范圍內(nèi)達(dá)到,則不必進(jìn)行嶺嵴分析;否則就要進(jìn)行).例6.7.2求面包酵母高產(chǎn)和高發(fā)酵活力的最優(yōu)工藝參數(shù).指標(biāo):y1---面包酵母產(chǎn)量[指最終發(fā)酵液的菌體濃度(g/L)]y2---發(fā)酵活力[指面團(tuán)從水中浮起的時(shí)間(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下面給出程序運(yùn)行得到的部分結(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)定點(diǎn)不是理想點(diǎn)就要進(jìn)一步作嶺嵴分析,請(qǐng)看示意圖和例子演示(例見datarubber或datapectic)例6.7.3橡膠樹幼苗對(duì)硫35的吸收實(shí)驗(yàn)(2×3析因試驗(yàn))因素:溫度t(℃)---7℃,20℃

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

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

to右端點(diǎn)

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

to右端點(diǎn)

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))是非線性的且無法線性化的模型.此時(shí)SAS采用nlin非線性回歸程序程序形式:procnlindata=…;model因變量

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

parms

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

…;

der.參數(shù)1=偏導(dǎo)表達(dá)式;der.參數(shù)2=偏導(dǎo)表達(dá)式;……………;runs;方法:采用迭代求參數(shù),要給出參數(shù)的一階偏導(dǎo)表達(dá)式der.參數(shù)=…;,要賦參數(shù)初值parms參數(shù)=…;請(qǐng)看例子演示例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進(jìn)行非線性回歸*/zz=exp(b-c*x);/*簡化記號(hào)*/term=exp(-zz);/*簡化記號(hào)*/model

y=a*term;

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

DependentVariableyMethod:Gauss-NewtonSumofIterbcaSquares01.10000.200070.00009161.21-0.1190