/* Example 4.4.1.----H0: mu1=mu2 */
proc iml;
stridec={131.5,  145,  141,  150};
stridep={ 40.5,   80,   50,   90};
strainc={    9,   12,   30,   36};
strainp={   54, 74.5, 64.5, 60.5};
x=stridec||strainc;
y=stridep||strainp;
p=2; n1=4;n2=4; f=n1+n2-2;
xb=J(n1,1,1)*J(1,n1,1/n1)*x;
yb=J(n2,1,1)*J(1,n2,1/n2)*y;
s1=(x-xb)`*(x-xb)/(n1-1);
s2=(y-yb)`*(y-yb)/(n2-1);
sp=((n1-1)*s1+(n2-1)*s2)/(n1+n2-2);
Tsq=(n1*n2)/(n1+n2)*(xb-yb)[1,]*inv(sp)*(xb-yb)[1,]`;
Tasq=f*p/(f-p+1)*finv(1-0.05,p, f-p+1);
F0=(f-p+1)/(f*p)*Tasq;
pv=1-probf(F0,2,5);
t11=tinv(1-0.05/4, 6)*sqrt((n1+n2)/(n1*n2)*sp[1,1]);
t21=tinv(1-0.05/4, 6)*sqrt((n1+n2)/(n1*n2)*sp[2,2]);
t12=sqrt(Tasq*(n1+n2)/(n1*n2)*sp[1,1]);
t22=sqrt(Tasq*(n1+n2)/(n1*n2)*sp[2,2]);
low11=(xb-yb)[1,1]-t11;
upp11=(xb-yb)[1,1]+t11;
low21=(xb-yb)[1,2]-t21;
upp21=(xb-yb)[1,2]+t21;
low12=(xb-yb)[1,1]-t12;
upp12=(xb-yb)[1,1]+t12;
low22=(xb-yb)[1,2]-t22;
upp22=(xb-yb)[1,2]+t22;
print tsq Tasq pv;
if tsq>Tasq then print "Tsq > Tasq:   reject";
            else print "Tsq < Tasq:   accept";
print  " Tasq at xb1-yb1: (" low12 upp12")",
       " Tasq at xb2-yb2: (" low22 upp22")",
       " Tfa2k at xb1-yb1: (" low11 upp11")",
       " Tfa2k at xb2-yb2: (" low21 upp21")";
quit;
           TSQ        TASQ          PV
     116.70813   13.886724   0.0000637
            Tsq > Tasq:   reject
                        LOW12       UPP12
 Tasq at xb1-yb1:  (30.363375   123.13662)
                        LOW22       UPP22
 Tasq at xb2-yb2: (-71.08426   -12.16574)
                        LOW11       UPP11
 Tfa2k at xb1-yb1: (39.796359   113.70364)
                        LOW21       UPP21
 Tfa2k at xb2-yb2:  (-65.09355   -18.15645)