SASIML交互式矩陣語言模塊原理與實現(xiàn)

1、桂林電子科技大學(xué)數(shù)學(xué)與計算科學(xué)學(xué)院創(chuàng)新基地統(tǒng)計部SAS/IML交互式矩陣語言模塊版權(quán)：桂林電子科技大學(xué)數(shù)學(xué)與計算科學(xué)學(xué)院創(chuàng)新基地統(tǒng)計部2011-10-18桂林電子科技大學(xué)數(shù)學(xué)與計算科學(xué)學(xué)院創(chuàng)新基地統(tǒng)計部SAS/IML模塊與矩陣運算語言我們的學(xué)習(xí)目的：1 .了解SAS/IML模塊的功能和特點；2 .掌握SAS/IML語言建立矩陣變量及矩陣的運算；3 .掌握IML常用函數(shù)的調(diào)用；4 .掌握自定義函數(shù)與過程的定義與調(diào)用；1、SAS/IML模塊1.1 SAS/IML模塊簡介IML是InteractiveMatrixLanguage的縮寫，即交互式矩陣語言，也是一種程序設(shè)計語言，它提供了條件轉(zhuǎn)移、循環(huán)等

2、流程控制語句，還提供了用戶自定義函數(shù)和子程序的功能，它是用戶研究新算法或解決系統(tǒng)中沒有現(xiàn)成方法可用的問題的工具。1.2 IML變量IML的基本數(shù)據(jù)元素是一個矩陣，最小為1X1階矩陣。矩陣元素可為數(shù)據(jù)型，也可為字符型，缺失信用”.”表示。多維矩陣每一行的元素個數(shù)必須相同。(1)數(shù)值型矩陣中元素可以是整數(shù)、小數(shù)，也可用科學(xué)計數(shù)法表示。(2)字符型矩陣中，若元素中含有空格或某些特殊字符，如“？”、"=”、“*”、“：”等，則須將字符串用單引號或雙引號括起來。在同一個字符型矩陣中，每一個元素的字符串長度不必相等，但存儲的長度最大不超過32676.(3) IML變量的取名規(guī)則： .可由下劃線，

3、字母和數(shù)字組成，但不能以數(shù)字開頭。 .不可超過32個字符 .不區(qū)分大小寫1.3 啟動和退出SAS/IML模塊啟動語句為：ProcIML;退出語句為：Quit;1.4 設(shè)置和關(guān)閉IML交互模型設(shè)置交互模式SAS/IML模塊提供了一個交互式環(huán)境，在此環(huán)境下，用戶按照輸入IML語句，然后系統(tǒng)將交互地輸出每一執(zhí)行語句的計算結(jié)果而不需額外使用輸出語句，這在調(diào)試程序時非常有用。使用這一功能前需將交互模型狀態(tài)打開，設(shè)置交互模型的語句為：Resetprint;關(guān)閉交互模式如果不需要一些中間計算結(jié)果，則應(yīng)將IML交互模型關(guān)閉，使得程序輸出比較簡潔。關(guān)閉交互模型的語句為桂林電子科技大學(xué)數(shù)學(xué)與計算科學(xué)學(xué)院創(chuàng)新基地統(tǒng)

4、計部ResetNoPrint;2 .矩陣變量的建立2.1 一般方式建立矩陣格式：矩P$變量名=矩陣初值;(1)定義1x1階矩陣(標(biāo)量)如a1=12;a2=.;輸出矩陣*/a3='Hello'printala2a3;/*(2)定義行向量如b1=12345);b2=abCDDEF);b3='adCD'DEF);/*比較b2與b3的結(jié)果*/(3)定義列向量如c1=1,2,3);(4)定義2X3階矩陣如c2=123,456);詳細(xì)見程序如下:prociml;resetprint;/*定義1x1階矩陣(標(biāo)量)*/a1=12;a2=.;a3='Hello'/

5、*定義行向量*/b1=12345);b2=abCDDEF);b3='adCD'DEF);/*比較b2與b3的結(jié)果*/*定義列向量*/c1=1,2,3);/*定義2X3階矩陣*/c2=123,456);quit;run;輸出結(jié)果如下圖所示:桂林電子科技大學(xué)數(shù)學(xué)與計算科學(xué)學(xué)院創(chuàng)新基地統(tǒng)計部blIrev5colx(mineric)12345b21row2cols(character,size4)ABCDDEFb31row2coI(chmricter,sise4)adCDDEFcl3rars1col(numeric)2.2 通過賦值建立矩陣變量格式：矩P$變量名=IML表達(dá)式;如：X=

6、1,2,3;丫二3,2,1;Z=X+Y;2.3 行向量的特殊產(chǎn)生方式格式1:矩陣變量名=起始值：結(jié)束值；說明：矩陣的第一個元素值為起始值，起始值小于結(jié)束值時，矩陣元素依次增加1,直到最后一個元素的值小于等于結(jié)束值，反之依次減少1,直到最后一個元素的值大于等于結(jié)束值格式2:矩陣變量名=do（起始值，結(jié)束值，步長）；程序如下：prociml;resetprint;i=1:10;i=10:1;i=do(1.1,10,2);i=do(10,0,-2);quit;run;輸出結(jié)果如下圖:桂林電子科技大學(xué)數(shù)學(xué)與計算科學(xué)學(xué)院創(chuàng)新基地統(tǒng)計部i1rox10cals(numeric)1234567B910i1ro

7、v10cala(numeric)10337654321i1row5cols(numeric)1J9.1i1row6cols(nuneric)10Sg4203 .矩陣運算3.1 矩陣的算術(shù)運算除以通常意義下的矩陣相加、減和乘運算外，SAS/IML還定義了矩陣的除、乘方、對應(yīng)元素相乘、對應(yīng)元素乘方等算術(shù)運算。詳見下表（A、B為矩陣，C為向量，k為實數(shù)）運算表達(dá)式含義加法A+BA+k矩陣A中兀素與矩陣B中兀素對應(yīng)相加矩陣A中元素都加上數(shù)k減法A-BA-k矩陣A中兀素與矩陣B中兀素對應(yīng)相減矩陣A中元素都減去數(shù)k乘法A*BA*k笙固A.f?矩陳一B矩陣A中兀素都乘以實數(shù)k除法A/BA/

8、k矩陣A中兀素與矩陣B中兀素對應(yīng)相除矩陣A中元素都除以數(shù)k乘方A*kA*(-1)矩陣A的k次哥，即為k個矩陣A相乘儕次k為大于等于-1的整數(shù)，若非整數(shù)將自動取整求A的逆矩陣對應(yīng)元I-t系兒系相乘A#BA#CA#k矩陣A中兀素與矩陣B中兀素對應(yīng)相乘矩陣A中每一行或每一列中的元素與向量C中元素對應(yīng)相乘（注：向量C可為行向量也可為列向量）矩陣A中元素都乘以數(shù)k對應(yīng)元I-t系兒系乘方A#BA#k矩陣A中兀素為底，矩陣B中兀素為指數(shù)作乘方運算矩陣A中元素為底，數(shù)k為指數(shù)作乘方運算兒系取相反數(shù)-A將矩陣A中兀素均取相反數(shù)桂林電子科技大學(xué)數(shù)學(xué)與計算科學(xué)學(xué)院創(chuàng)新基地統(tǒng)計部矩陣的基本運算程序如下:prociml

9、;resetprint;A=12,34);B=43,21);C=A+B;D=A+2;e=A/B;f=A/2;g=A*2;h=A*(-1);j=A#B;m=A#2;n=A#B;o=A#2;quit;run;|輸出結(jié)果如下圖：C2rows2cols(numeric)55552rows2cols(njitieric)34562rows2cols(numeric)0,25o.eeeGBe?LB42rows2cols(numeric)0.511.522rows2cols(numeric)7II01522桂林電子科技大學(xué)數(shù)學(xué)與計算科學(xué)學(xué)院創(chuàng)新基地統(tǒng)計部3.2 矩陣的比較運算矩陣的比較運算是指比較矩陣相應(yīng)元

10、素大小關(guān)系，結(jié)果產(chǎn)生一個與相比較矩陣同階的矩陣。若比較關(guān)系成立，則結(jié)果矩陣中的相應(yīng)元素為1,否則為00運算表達(dá)式含義A>BA<B小于A>=BA<=B小于等于A=BAA=B不等于注：上述矩陣A和B中某一個可用一實數(shù)代替，此時將另一矩陣的每一個元素與此數(shù)值進(jìn)行比較。詳見程序如下。prociml;resetprint;A=137,642;B=082,413;C=A>B;D=a>2;quit;輸出結(jié)果如下圖：C2rows3cols(numeric)2rows3cols(numeric)3.3 矩陣的合并與拆分將兩個滿足一定條件的矩陣按照一定規(guī)則合成一個矩陣的過程為矩

11、陣的合并，矩陣合并有水平合并與垂直合并之分。從一個矩陣中分出一個階數(shù)較小的矩陣的過程稱為矩陣的拆分。運算表達(dá)式含義A|B將矩陣A中的每行兀素分別接在矩陣B中相應(yīng)行的右邊，結(jié)果矩陣的行數(shù)與A,B相同A/B將矩陣A中的每列兀素分別接在矩陣B中相應(yīng)列的下邊，結(jié)果矩陣的列數(shù)與A,B相同Ai,j矩陣A的第i行第j列的兀素Ai,矩陣A的第i行所有兀素桂林電子科技大學(xué)數(shù)學(xué)與計算科學(xué)學(xué)院創(chuàng)新基地統(tǒng)計部A,j矩陣A的第j列所有兀素Ac1：c2,d1：d2矩陣A的第c1行到第c2行和第d1列到第d2列的所有元素A+,矩陣A的各列兀素之和所形成的行向量A,+矩陣A的各列兀素之和所形成的列向量詳見程序如下procim

12、l;resetprint;A=111;777;B=000;888;C=A|B;D=A/B;x=111213,212233;E=x1,2;F=x+,;G=x,+;quit;run;輸出結(jié)果如下圖所示：08A2rows3cols(numeric)11177?B2rows3cols(numeric)c2roffs6cols(numeric)111000777e88D4rows3cols(numeric)11177?0D08882ro>53cols(numeric)111213212233桂林電子科技大學(xué)數(shù)學(xué)與計算科學(xué)學(xué)院創(chuàng)新基地統(tǒng)計部3.4 其它常用矩陣運算運算表達(dá)式含義I(n)產(chǎn)一個n階單位

13、陣J(m,n,k)產(chǎn)一個兀素全為k的mxn階矩陣Diag(C)以向量C的各兀素為主對角線兀素（順序不變）的對角矩陣Trace(A)矩陣A的跡（主對角線兀素之和）Inv(A)或A*(-1).一-1矩陣A的逆矩陣，即AT(A)或A'矩陣A的轉(zhuǎn)置（為右斜撇，一般為鍵盤的左上方）Det(A)矩陣A的行列式值Eigval(A)矩陣A的全部特征值Eigvec(A)矩陣A的每個特征向量為列所形成的矩陣（單位化）例詳見程序如下prociml;resetprint;A1=i(3);A2=j(2,3,2);C=123456;A3=diag(C);A=1-1,24;B1=trace(A);B2=inv(A)

14、;B3=det(A);B4=A'B5=eigval(A);B6=eigvec(A);quit;run;桂林電子科技大學(xué)數(shù)學(xué)與計算科學(xué)學(xué)院創(chuàng)新基地統(tǒng)計部1row8col5(numeric)2rows2si不(numeric)611row1col(numeric)B22rows2cols(minerIc)0.C666S670.166667-0.33333：30J6GGGB74 .矩陣函數(shù)4.1 常用數(shù)學(xué)函數(shù)函數(shù)名功能含義Abs(A)將矩陣A中的每個兀素取絕對值Exp(A)以自然數(shù)e為底，將矩陣A中的每個兀素作為指數(shù)，構(gòu)造同階矩陣Int(A)將矩陣A中的每個兀素取整數(shù)Log(A)將矩陣A中的

15、每個兀素取對數(shù)Max(A)將矩陣A中的每個兀素取最大值Min(A)將矩陣A中的每個兀素取小值Sqrt(A)將矩陣A中的每個兀素開平方Sum(A1,A2,.)將矩陣A1,A2,.中所有矩陣的所有兀素求和，自變量矩陣最多為15個，如果有缺失將忽略它們。桂林電子科技大學(xué)數(shù)學(xué)與計算科學(xué)學(xué)院創(chuàng)新基地統(tǒng)計部4.2 常用概率統(tǒng)計函數(shù)函數(shù)名功能含義SSQ(A1,A2,.)CuSum(A)求矩陣A1,A2,.中所有矩陣的所有兀素的平方和，自變量矩陣最多為15個。產(chǎn)一個與自變量矩陣A向階的矩陣，其兀素為A中兀素的累積和（按行為從左到右，再從上到下）ProbNorm(x)計算標(biāo)準(zhǔn)正態(tài)分布的分布函數(shù)在x處的值，x為數(shù)

16、值型變量Probit(p)計算標(biāo)準(zhǔn)正態(tài)分布的p分位數(shù),0<p<1ProbChi(x,df<,nc>)計算自由度為df的72分布的分布函數(shù)在x處的值。當(dāng)非中心參數(shù)nc=0或缺省時，即為中心分布的X2分布。Cinv(p,df<,nc>)計算自由度為df的72分布的p分位數(shù)，0wp<1,非中心參數(shù)nc>0.ProbT(x,df<,nc>)計算自由度為df的T分布的分布函數(shù)在x處的值。當(dāng)非中心參數(shù)nc=0或缺省時，即為中心分布的T分布。Tinv(p,df<,nc>)計算自由度為df的T分布的p分位數(shù)，0<p<1,非中心

17、參數(shù)nc>0.ProbF(x,df1,ddf2<,nc>)計算分子、分母自由度為分別為df1和df2的F分布的分布函數(shù)在x處的值。當(dāng)非中心參數(shù)nc=0或缺省時，即為中心分布的T分布。Finv(p,df1,ddf2<,nc>)計算分子、分母自由度為分別為df1和df2的F分布的p分位數(shù)，0<p<1,非中心參數(shù)nc>0.RanNor(seed)產(chǎn)生標(biāo)準(zhǔn)正態(tài)分布隨機數(shù)。種子seed可是一個整數(shù)或一個矩陣，元素值小于231-1,seedw0時將用系統(tǒng)時間作為初始化的種子。若seed為矩陣，則函數(shù)值為與seed同階的矩陣。Uniform(seed)產(chǎn)生（0

18、,1）上均勻分布隨機數(shù)。種子seed可是一個整數(shù)或一個矩陣，元素值小于231-1,seedw0時將用系統(tǒng)時間作為初始化的種子。若seed為矩陣，則函數(shù)值為與seed同階的矩陣。注：其功能等價于數(shù)據(jù)步中的ranuni（）函數(shù)。應(yīng)用舉例如下：prociml;resetprint;x=01;y=0.5;a=ssq(x,y);b1=cusum(123b2=cusum(56,3c1=Probnorm(x);c2=probit(y);d1=probT(x,100);d2=Tinv(y,100);4);4);桂林電子科技大學(xué)數(shù)學(xué)與計算科學(xué)學(xué)院創(chuàng)新基地統(tǒng)計部e1=10+2*rannor(x);quit;run

19、;1row1col(numeric)bl1row1,254cols(numeric)b212rows32colsE10(numeric)cl1row5142cols11IS(numeric)1row0,50J4134471col(numeric)di1row-L16E-172cols(numeric)0.50,84013795 .IML語句與SAS語句相同，類似于C語言，每IML語句以“結(jié)尾，但不區(qū)分大小寫5.1 循環(huán)語句1)計數(shù)型循環(huán)格式：Do循環(huán)變量=初值to終值By步長;循環(huán)體；End;例如：doi=1to5by2；Print'ThevalueofIis:'I；End;2

20、) .直到型循環(huán)桂林電子科技大學(xué)數(shù)學(xué)與計算科學(xué)學(xué)院創(chuàng)新基地統(tǒng)計部格式：DoUntil(表達(dá)式)；循環(huán)體；End;3) .當(dāng)型循環(huán)格式：DoWhile(表達(dá)式)；循環(huán)體；End;應(yīng)用舉例如下：prociml;doi=1to5by2;print'ThevalueofIis:'I;end;quit;prociml;x=1;|dountil(x>100);x=x+1;end;printx;|quit;prociml;x=123;ifmax(z)<2thenp=0;elsep=1;printp;quit;run;iThevalueofIis：1IThevalueofIis：3

21、IIThevaluecfIis；5101桂林電子科技大學(xué)數(shù)學(xué)與計算科學(xué)學(xué)院創(chuàng)新基地統(tǒng)計部5.2 條件判斷語句格式：If表達(dá)式then語句1；Else語句2;5.3 退出語句格式：Quit;執(zhí)行此語句退出IML模塊。5.4 返回語句格式：Return(表達(dá)式)；執(zhí)行此語句后，將計算表達(dá)式的值并返回此值。此語句常用于自定義過程或函數(shù)中。5.5 過程(子程序)調(diào)用語句格式：Run過程名(變量列表)；6.IML自定義函數(shù)和過程6.1 自定義一般格式START過程名或函數(shù)名(形式變量列表)Global(全局變量列表);函數(shù)或過程中語句；Finish過程名；6.2 自定義過程示例/*定義一個名為mymod

22、的過程，其功能是計算a*x+b,其中a和b為形勢變量，x是全局變量，最終結(jié)果通過全局變量y帶出來*/prociml;x=3;startmymod(a,b)global(x,y);y=a*x+b;|finish;runmymod(1,2);printy;|quit;run;6.3 自定義函數(shù)示例startLSE(x,y);b=inv(x'*x)*x'*y;return(b);finish;上述例子定義了一個名為LSE的函數(shù)，其功能是計算y關(guān)于x的線性回歸方程系數(shù)的最小二乘估計，結(jié)果以函數(shù)值的形式返回桂林電子科技大學(xué)數(shù)學(xué)與計算科學(xué)學(xué)院創(chuàng)新基地統(tǒng)計部6.4 線性回歸IML模塊程序線性

23、回歸IML模塊程序prociml;x=111,124,139,1416,1525);y=1,5,9,23,36);b=inv(t(x)*x)*t(x)*y;yhat=x*b;r=y-yhat;sse=ssq(r);dfe=nrow(x)-ncol(x);mse=sse/dfe;startregress;/*定義模塊開始*/xpxi=inv(t(x)*x);/*矩陣X'X的逆*/beta=xpxi*(t(x)*y);/*參數(shù)估計*/yhat=x*beta;/*預(yù)測值*/resid=y-yhat;/*殘差*/sse=ssq(resid);/*SSE*/n=nrow(x);/*觀測值數(shù)*/d

24、fe=nrow(x)-ncol(x);/*誤差自由度DF*/mse=sse/dfe;/*MSE*/cssy=ssq(y-sum(y)/n);/*校正平方和*/rsquare=(cssy-sse)/cssy;/*RSQUARE*/print,"RegressionResults",ssedfemsersquare;stdb=sqrt(vecdiag(xpxi)*mse);/*參數(shù)估計的標(biāo)準(zhǔn)差*/t=beta/stdb;/*參數(shù)的t檢驗*/prob=1-probf(t#t,1,dfe);/*p-值*/print,"ParameterEstimates",be

25、tastdbtprob;print,yyhatresid;finishregress;/*模塊結(jié)束*/resetnoprint;runregress;/*執(zhí)行REGRESS模塊*/resetprint;/*打開自動打印*/covb=xpxi*mse;/*參數(shù)估計協(xié)方差陣*/s=1/sqrt(vecdiag(covb);|corrb=diag(s)*covb*diag(s);/*參數(shù)估計的相關(guān)系數(shù)*/x1=1,2,3,4,5);/*矩陣X的第二列*/x=orpol(x1,2);/*產(chǎn)生正交多項式*/桂林電子科技大學(xué)數(shù)學(xué)與計算科學(xué)學(xué)院創(chuàng)新基地統(tǒng)計部resetnoprint;/*關(guān)閉自動打印*/ru

26、nregress;/*運行REGRESS*/resetprint;covb=xpxi*mse;|s=1/sqrt(vecdiag(covb);corrb=diag(s)*covb*diag(s);resetfuzz;corrb=diag(s)*covb*diag(s);x=orpol(x1,2);注釋；產(chǎn)生矩陣a的正交多項式6.5綜合例子某長途運輸公司對同類卡車統(tǒng)計了行行駛里程y與行駛天數(shù)x的數(shù)據(jù)如表:x3.51.04.0y8252151070希望通過此數(shù)據(jù)建立行駛里程SAS/IML程序如下:2.01.05504803.04.51.590013503253.06705.01215y關(guān)于行駛天數(shù)

27、x的線性回歸方程。prociml;startLSE(x,y);/*定義函數(shù)LSE用于求統(tǒng)計回歸方程*/|b=inv(x'*x)*x'*y;/*求統(tǒng)計回歸方程的最小二乘估計*/return(b);/*以函數(shù)值的形式返回求得系數(shù)的最小二乘估計*/finishstartmain;/*定義主程序*/x=13.5,11.0,14.0,12.0,11.0,13.0,14.5,11.5,13.0,15.0;y=825,215,1070,5500,480,900,1350,325,670,1251;beat=LSE(x,y);printbeat;finish;run;/*運行主程序*/quit

28、;輸出結(jié)果如下：beat1178.815428.064777桂林電子科技大學(xué)數(shù)學(xué)與計算科學(xué)學(xué)院創(chuàng)新基地統(tǒng)計部7.打開與使用SAS數(shù)據(jù)集7.1 打開一個SAS數(shù)據(jù)集在訪問一個SAS數(shù)據(jù)集前，你必須發(fā)布一個命令來打開它。如果只是從已存在的數(shù)據(jù)集中簡單地讀取數(shù)據(jù)，可使用USE語句來打開數(shù)據(jù)集，一般格式為：USE數(shù)據(jù)集名VAR變量名WHERE（達(dá)式）;例：prociml;usesashelp.class;listall;/*顯示數(shù)據(jù)集的觀測值*/quit;7.2 由SAS數(shù)據(jù)集創(chuàng)建矩陣可以使用REA曲句，從SASa據(jù)集傳遞數(shù)據(jù)到矩陣?yán)铮⒁庠摂?shù)據(jù)集必須先打開，READg句的一般格式為READ范圍VAR

29、變量名WHERE俵達(dá)式）INTO矩陣名;prociml;usesashelp.class;|listall;/*顯示數(shù)據(jù)及觀測值*/readallvarnamesexage;/*把數(shù)據(jù)集class中的變量namesexage讀取出來并形成一個矩陣*/printnamesexage;quit;prociml;usesashelp.class;|listall;readallvarnamesexageintotemp;|printtemp;quit;run;注：第二段程序執(zhí)行后會在Log窗口中顯示警告信息：WARNIN&ariableAGEexcludedbecauseoftypemism

30、atchwithA.,提示age變量沒有導(dǎo)入矩陣temp中，因為它的數(shù)據(jù)類型與其它兩個變量不一致。思考：如果只要求將女性的觀測讀取出來，則需加入where選項。readallvarnamesexagewhere（sex='F'）;7.3 由矩陣創(chuàng)建SAS數(shù)據(jù)集SAS/IML提供了從矩陣創(chuàng)建一個新SA繳據(jù)集的能力，主要使用create和append語句來實現(xiàn)，有兩種方式分別來創(chuàng)建只有一個變量和多個變量數(shù)據(jù)集。格式為：方式一：CREAT敢據(jù)集名FROM巨陣名COLNAM網(wǎng)名ROWNAMfi;AppendFrom矩陣名；桂林電子科技大學(xué)數(shù)學(xué)與計算科學(xué)學(xué)院創(chuàng)新基地統(tǒng)計部方式二：CREATE數(shù)據(jù)集名VAR變量名;Append;例：/*創(chuàng)建只有一個變量的SAS數(shù)據(jù)集*/prociml;usesashelp.class;|listall;readallvarheighweight;/*從數(shù)據(jù)集中讀取變量，形成矩陣*/printheightweight;htwt=height/weight;