/*Example 5.3.1.---        */
proc iml;
x1={5,4,5,5,6,3};
x2={10,9,9,9,11,9};
x3={16,13,14,14,15,12};
x=x1||x2||x3;
n=6; p=3; f=n-1; alpha=0.05;
xbi=J(1,n,1/n)*x;
xbj=J(1,p,1/p)*x`;
xb=J(1,3,1/3)*xbi`;
v=0;
do i=1 to p;
do j=1 to n;
v=v+(x[j,i]-xbi[1,i]-xbj[1,j]+xb)**2;
end;
end;
v1=n*(xbi-xb)*(xbi-xb)`;
v2=p*(xbj-xb)*(xbj-xb)`;
F0=(n-1)*v1/v;
Fa=finv(1-alpha, p-1, (n-1)*(p-1));
print  F0 Fa ;
if F0>Fa then print "F0 > Fa:  reject";
         else print "F0 < Fa:  accept";
total=v+v1+v2;
r=3.877*sqrt(v/(n*f*(n-1)));
u1=xbi[1,3]-xbi[1,2];
u2=xbi[1,2]-xbi[1,1];
u3=xbi[1,3]-xbi[1,1];
low1=u1-r; upp1=u1+r;
low2=u2-r; upp2=u2+r;
low3=u3-r; upp3=u3+r;
print "mu3-mu2:  ("low1","upp1")",
      "mu2-mu1:  ("low2","upp2")",
      "mu3-mu1:  ("low3","upp3")";
quit;
               F0           FA
        336.14286     4.102821
             F0 > Fa:  reject
         mu3-mu2:  ( 3.8757437 , 5.1242563 )
         mu2-mu1:  (  4.209077 , 5.4575896 )
         mu3-mu1:  (  8.709077 , 9.9575896 )