当前位置：文档库 › SAS例题及程序输出6

SAS例题及程序输出6

地质勘探中，在A,B,C 三个地区采集了一些岩石，测量其部分化学成分，其

数据见表3.5。假定这三个地区掩饰的成分遵从()3,(1,2,3)(0.05)i i N i μα∑==（）

。

（1）检验不全01231123:=:,,H H ∑=∑∑∑∑∑；不全等；（2）检验(1)(2)(1)(2)01::H H μμμμ=≠；；

（3）检验(1)(2)(3)()()01::,i j H H i j μμμμμ==≠≠；存在使。

(1)检验假设

01231123:=:,,H H ∑=∑∑∑∑∑；不全等，

在H 0成立时，取近似检验统计量为2()f χ 统计量：

()()*4=121ln d M d ξλ-=--。

由样本值计算三个总体的样本协方差阵：

1(1)(1)(1)(1)

11()()

1111110.243081=0.642649.2855240.014060.020520.00452n S A X X X X n n ααα='==----?? ?- ? ???

∑()()，

1(2)(2)(2)(2)

23()()12211116.30461= 4.756710.672230.05570.23880.006675n S A X X X X

n n ααα='==----?? ?- ? ?-??∑()()， 1(3)(3)(3)(3)

33()()1

3311112.97141=0.63370.342140.00010.002950.001875n S A X X X X

n n ααα='==----?? ? ? ?-??

∑()()。

进一步计算可得

1231

0.0018318,0.0000942,0.0011851,0.0000417,10

S A S S S =

==== 24.52397,0.433333,12,M d f ===

(1)=13.896916d M ξ=-。

对给定显著性水平=0.05α，利用软件SAS9.3进行检验时，首先计算p 值：

p =P {ξ≥13.896916}=0.3073394。因为p 值=0.3073394>0.05，故接收0H ，即认为方差阵之间无显著性差异。

proc iml ; n1=5;n2=4;n3=4; n=n1+n2+n3;k=3;p=3; x1={47.22 5.06 0.1, 47.45 4.35 0.15, 47.52 6.85 0.12, 47.86 4.19 0.17, 47.31 7.57 0.18 };

x2={54.33 6.22 0.12, 56.17 3.31 0.15, 54.4 2.43 0.22, 52.62 5.92 0.12 };

x3={43.12 10.33 0.05, 42.05 9.67 0.08, 42.5 9.62 0.02, 40.77 9.68 0.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;print a1;

a2=x2`*mm2*x2;print a2;

a3=x3`*mm2*x3;print a3;

tt=xx`*mm*xx;print tt;/*总离差阵*/

a=a1+a2+a3;print a;/*组内离差阵*/

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; /*统计量值*/

print da da1 da2 da3 m dd df;

p0=1-probchi(kc,df); /*显著性概率*/

print kc p0;

quit;

(2) 提出假设

(1)(2)(1)(2)01::H H μμμμ=≠，。

取检验统计量为

+1(3,6,9)(2)

n m p F T

p n m n m --=

===+-，

由样本值计算得：

1=(47.472.5.604,0.144)=(54.38,4.47,0.1525)X X ''（）（2）

，

， 120.24308=0.64264

9.285520.014060.020520.004526.3046= 4.7567

10.67220.05570.2388

0.006675A A ??

?- ? ??

?- ? ?-?

，

进一步计算得：

211112(2)()'()()=60.666995D n m X X A A X X -=+--+-（）（2）（）（2），

22134.81554,nm

T D n m

=+ 2

132.098939(2)n m p F T n m p

+--=

=+-。

对给定显著性水平=0.05α，利用软件SAS9.3进行检验时，首先计算p 值：

p =P {F ≥32.098939}=0.0010831。因为p 值=0.0010831<0.05，故否定0H ，即认为A ，B 两地岩石化学成分数据存在显著性差异。在这种情况下，可能犯第一类错误，且犯第一类错误的概率为0.05。

SAS 程序及结果如下：

proc iml;

n=5;m=4; p=3;

x={ 47.22 5.060.1,

47.45 4.350.15,

47.52 6.850.12,

47.86 4.190.17,

47.317.570.18

} ;

ln={[5] 1} ;

x0=(ln*x)`/n; print x0;

mx=i(n)-j(n,n,1)/n;

a1=x`*mx*x; print a1;

y={ 54.33 6.220.12,

56.17 3.310.15,

54.4 2.430.22,

52.62 5.920.12

} ;

lm={[4] 1} ;

y0=(lm*y)`/m; print y0;

my=i(m)-j(m,m,1)/m;

a2=y`*my*y; print a2;

a=a1+a2; xy=x0-y0;

ai=inv(a); print a ai;

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); print d2 t2 f beta;

pp=1-probf(f,p,m+n-p-1);

print pp; quit;

(3) 检验假设

(1)(2)(3)()()01::,i j H H i j μμμμμ==≠≠；存在使；

因似然比统计量~(,,1)p n k k ΛΛ-- ，本题中k-1=2，可以利用Λ统计量与F 统计量的关系，去检验统计量为F 统计量：

3,3,13),F k p n =

===

由样本值计算得：47.947696.5538460.11692=)3(,X '，，

及 (1)

(2)(3)47.47254.3842.115.604 4.479.8250.1440.15250.047,,5X X X ????????????===??????????????????

，

(1)()(1)()

123()()11

(1)(1)

()()11()()9.51908= 4.7656420.299820.069660.215330.01307=()()312.46343132.506284.9823082.5417077 1.5488460.0410769t

n t t t n t A A A A X X X X

T X X X X αααααα===='=++=--??

??-??

'--?=--?∑∑∑∑?

????

???

，

进一步计算得：

1.8318441

=0.0160379114.21942

A T

Λ=

，

22134.81554,nm

T D n m

810.126641

18.390234

30.126641f -=

==。

对给定显著性水平=0.05α，利用软件SAS9.3进行检验时，首先计算p 值：

p =P {F ≥18.390234}=2.3451×10-6

。因为p 值=2.3451×10-6<0.05，故否定0H ，即认为A ，B ，C 三地岩石化学成分数据存在显著性差异。在这种情况下，可能犯第一类错误，且犯第一类错误的

概率为0.05。

proc iml ; n1=5;n2=4;n3=4; n=n1+n2+n3;k=3;p=3; x1={47.22 5.06 0.1, 47.45 4.35 0.15, 47.52 6.85 0.12, 47.86 4.19 0.17, 47.31 7.57 0.18 };

x2={54.33 6.22 0.12, 56.17 3.31 0.15, 54.4 2.43 0.22, 52.62 5.92 0.12 };

x3={43.12 10.33 0.05, 42.05 9.67 0.08, 42.5 9.62 0.02, 40.77 9.68 0.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;print tt;/*总离差阵*/

a=a1+a2+a3; print a;/*组内离差阵*/

da=det(a);/*合并样本协差阵*/

dt=det(tt);

a0=da/dt;

print da dt a0;

b=sqrt(a0); print b;

f=(n-k-p+1)*(1-b)/(b*p);

df1=2*p;df2=2*(n-k-p+1);

p0=1-probf(f,df1,df2); /*显著性概率*/

print f p0;

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);

print fa f1 p1;

quit;