SAS例题及程序输出6

地质勘探中，在A,B,C三个地区采集了一些岩石，测量其部分化学成分，其数据见表3.5。假定这三个地区掩饰的成分遵从3,(1,2,3)(0.05)iiNi（）。（1）检验不全01231123:=:,,HH；不全等；（2）检验(1)(2)(1)(2)01::HH；；（3）检验(1)(2)(3)()()01::,ijHHij；存在使。表3.5岩石部分化学成分数据SiO2FeOK2OA地区47.225.060.1047.454.350.1547.526.850.1247.864.190.1747.317.570.18B地区54.336.220.1256.173.310.1554.402.430.2252.625.920.12C地区43.1210.330.0542.059.670.0842.509.620.0240.779.680.04解：(1)检验假设01231123:=:,,HH；不全等，在H0成立时，取近似检验统计量为2()f统计量：*4=121lndMd。由样本值计算三个总体的样本协方差阵：1(1)(1)(1)(1)11()()11111110.243081=0.642649.2855240.014060.020520.00452nSAXXXXnn()()，1(2)(2)(2)(2)23()()12211116.30461=4.756710.672230.05570.23880.006675nSAXXXXnn()()，1(3)(3)(3)(3)33()()13311112.97141=0.63370.342140.00010.002950.001875nSAXXXXnn()()。进一步计算可得12310.0018318,0.0000942,0.0011851,0.0000417,10SASSS24.52397,0.433333,12,Mdf(1)=13.896916dM。对给定显著性水平=0.05，利用软件SAS9.3进行检验时，首先计算p值：p=P{ξ≥13.896916}=0.3073394。因为p值=0.30733940.05，故接收0H，即认为方差阵之间无显著性差异。prociml;n1=5;n2=4;n3=4;n=n1+n2+n3;k=3;p=3;x1={47.225.060.1,47.454.350.15,47.526.850.12,47.864.190.17,47.317.570.18};x2={54.336.220.12,56.173.310.15,54.42.430.22,52.625.920.12};x3={43.1210.330.05,42.059.670.08,42.59.620.02,40.779.680.04};xx=x1//x2//x3;/*三组样本纵向拼接*/mm1=i(5)-j(5,5,1)/n1;mm2=i(4)-j(4,4,1)/n2;mm=i(n)-j(n,n,1)/n;a1=x1`*mm1*x1;printa1;a2=x2`*mm2*x2;printa2;a3=x3`*mm2*x3;printa3;tt=xx`*mm*xx;printtt;/*总离差阵*/a=a1+a2+a3;printa;/*组内离差阵*/da=det(a/(n-k));/*合并样本协差阵*/da1=det(a1/(n1-1));/*每个总体的样本协差阵阵*/da2=det(a2/(n2-1));da3=det(a3/(n3-1));m=(n-k)*log(da)-(4*log(da1)+3*log(da2)+3*log(da3));dd=(2*p*p+3*p-1)*(k+1)/(6*(p+1)*(n-k));df=p*(p+1)*(k-1)/2;/*卡方分布自由度*/kc=(1-dd)*m;/*统计量值*/printdada1da2da3mdddf;p0=1-probchi(kc,df);/*显著性概率*/printkcp0;quit;(2)提出假设(1)(2)(1)(2)01::HH，。取检验统计量为2+1(3,6,9)(2)nmpFTpnmnm，由样本值计算得：1=(47.472.5.604,0.144)=(54.38,4.47,0.1525)XX（）（2），，120.24308=0.642649.285520.014060.020520.004526.3046=4.756710.67220.05570.23880.006675AA，，进一步计算得：211112(2)()'()()=60.666995DnmXXAAXX（）（2）（）（2），22134.81554,nmTDnm2132.098939(2)nmpFTnmp。对给定显著性水平=0.05，利用软件SAS9.3进行检验时，首先计算p值：p=P{F≥32.098939}=0.0010831。因为p值=0.00108310.05，故否定0H，即认为A，B两地岩石化学成分数据存在显著性差异。在这种情况下，可能犯第一类错误，且犯第一类错误的概率为0.05。SAS程序及结果如下：prociml;n=5;m=4;p=3;x={47.225.060.1,47.454.350.15,47.526.850.12,47.864.190.17,47.317.570.18};ln={[5]1};x0=(ln*x)`/n;printx0;mx=i(n)-j(n,n,1)/n;a1=x`*mx*x;printa1;y={54.336.220.12,56.173.310.15,54.42.430.22,52.625.920.12};lm={[4]1};y0=(lm*y)`/m;printy0;my=i(m)-j(m,m,1)/m;a2=y`*my*y;printa2;a=a1+a2;xy=x0-y0;ai=inv(a);printaai;dd=xy*ai*xy`;d2=(m+n-2)*dd;t2=n*m*d2/(n+m);f=(n+m-1-p)*t2/((n+m-2)*p);fa=finv(0.95,p,m+n-p-1);beta=probf(f,p,m+n-p-1,t2);printd2t2fbeta;pp=1-probf(f,p,m+n-p-1);printpp;quit;(3)检验假设(1)(2)(3)()()01::,ijHHij；存在使；因似然比统计量~(,,1)pnkk，本题中k-1=2，可以利用统计量与F统计量的关系，去检验统计量为F统计量：()11(3,3,13),nppFkpnp由样本值计算得：47.947696.5538460.11692=)3(,X，，及(1)(2)(3)47.47254.3842.115.6044.479.8250.1440.15250.047,,5XXX，3(1)()(1)()123()()113(1)(1)()()11()()9.51908=4.7656420.299820.069660.215330.01307=()()312.46343132.506284.9823082.54170771.5488460.0410769ttntttntAAAAXXXXTXXXX，进一步计算得：1.8318441=0.0160379114.21942AT，22134.81554,nmTDnm11810.12664118.39023430.126641nkpfp。对给定显著性水平=0.05，利用软件SAS9.3进行检验时，首先计算p值：p=P{F≥18.390234}=2.3451×10-6。因为p值=2.3451×10-60.05，故否定0H，即认为A，B，C三地岩石化学成分数据存在显著性差异。在这种情况下，可能犯第一类错误，且犯第一类错误的概率为0.05。prociml;n1=5;n2=4;n3=4;n=n1+n2+n3;k=3;p=3;x1={47.225.060.1,47.454.350.15,47.526.850.12,47.864.190.17,47.317.570.18};x2={54.336.220.12,56.173.310.15,54.42.430.22,52.625.920.12};x3={43.1210.330.05,42.059.670.08,42.59.620.02,40.779.680.04};xx=x1//x2//x3;/*三组样本纵向拼接*/ln={[5]1};lnn{[4]1};lnnn={[13]1};x10=(ln*x1)`/n1;x20=(lnn*x2)`/n2;x30=(lnn*x3)`/n3;xx0=(lnnn*x1)`/n1;mm1=i(5)-j(5,5,1)/n1;mm2=i(4)-j(4,4,1)/n2;mm=i(n)-j(n,n,1)/n;a1=x1`*mm1*x1;a2=x2`*mm2*x2;a3=x3`*mm2*x3;tt=xx`*mm*xx;printtt;/*总离差阵*/a=a1+a2+a3;printa;/*组内离差阵*/da=det(a);/*合并样本协差阵*/dt=det(tt);a0=da/dt;printdadta0;b=sqrt(a0);printb;f=(n-k-p+1)*(1-b)/(b*p);df1=2*p;df2=2*(n-k-p+1);p0=1-probf(f,df1,df2);/*显著性概率*/printfp0;f1=(tt[1,1]-a[1,1])*(n-k)/((k-1)*a[1,1]);p1=1-probf(f1,k-1,n-k);fa=finv(0.95,k-1,n-k);printfaf1p1;quit;