/* Example 4.5.2---- Subvector (Two sample) */
proc iml;
x1={12.3,12.1,12.8,12,12.1,11.8,12.7,12.5};
x2={2.5,2.2,2.9,2.1,2.2,1.9,2.9,2.7};
x3={2.9,2.5,3,2.2,2.4,2,3.3,3};
y1={12,11.8,12.7,12.4,12.1,12,11.7,12.2};
y2={2.3,2,3.1,2.8,2.5,2.2,2,2.5};
y3={2.7,2.4,3.6,3.2,2.8,2.7,2.4,3};
x=x1||x2||x3;
y=y1||y2||y3;
n1=8; n2=8; p=3; f=n1+n2-2; s=1; t=2; alpha=0.05;
xb=J(n1,n1,1/n1)*x;
yb=J(n2,n2,1/n2)*y;
ss=((x-xb)`*(x-xb)+(y-yb)`*(y-yb))/(n1+n2-2);
Tpsq=n1*n2/(n1+n2)*(xb[1,]-yb[1,])*inv(ss)*(xb[1,]-yb[1,])`;
Tssq=n1*n2/(n1+n2)*(xb[1,1]-yb[1,1])*inv(ss[1,1])*(xb[1,1]
     -yb[1,1]);
F0=(f-p+1)/t*(Tpsq-Tssq)/(f+Tssq);
Fa=finv(1-alpha,t,f-p+1);
print F0 Fa;
if F0>Fa then print "F0 > Fa:  reject";
         else print "F0 < Fa:  reject";
casq=t/(f-p+1)*Fa;
s21=ss[{2 3},{2 3}]-ss[{2 3},1]*inv(ss[1,1])*ss[{2 3},1]`;
tt=1/(f-s)*tinv(1-alpha/4, f-s)**2;
if tt<casq then c=tt; else c=casq;
m=xb[1,{2 3}]`-yb[1,{2 3}]`-ss[{2 3},1]*inv(ss[1,1])*(xb[1,1]
   -yb[1,1]);
e=1/sqrt(n1*n2/(n1+n2))*sqrt(f+Tssq)*sqrt(c);
low1=m[1,1]-e*sqrt(s21[1,1]);
upp1=m[1,1]+e*sqrt(s21[1,1]);
low2=m[2,1]-e*sqrt(s21[2,2]);
upp2=m[2,1]+e*sqrt(s21[2,2]);
print low1 upp1, low2 upp2;
quit;
                 F0           FA
          24.792586    3.8852938
              F0 > Fa:  reject
             LOW1         UPP1
       (-0.283964    -0.108212)
             LOW2      UPP2
       (-0.584771    -0.219924)