/*Example 4.4.3.----Behrens-Fisher Problem */
/*n1=10, n2=13; p=2, alpha=0.05    */
proc iml;
x1={180, 192, 217, 221, 171, 192, 213, 192, 170, 201};
x2={245, 260, 276, 299, 239, 262, 278, 255, 244, 276};
y1={181, 158, 184, 171, 181, 181, 177, 198, 180, 177, 176, 192,
    164};
y2={305, 237, 300, 273, 297, 308, 301, 308, 286, 299, 317, 312,
    265};
x=x1||x2; y=y1||y2;
n1=10; n2=13; p=2; alpha=0.05;
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);
e=(xb[1,]-yb[1,])*inv(s1/n1+s2/n2);
Tsq=e*(xb[1,]-yb[1,])`;
finv=(e*s1*e`)**2/(Tsq**2*(n1**3-n1**2))+
     (e*s2*e`)**2/(Tsq**2*(n2**3-n2**2));
f=round(1/finv,1);
T0sq=(f-p+1)/(f*p)*Tsq;
Ta=finv(1-0.05,p, f-p+1);
print  T0sq  Ta;
if T0sq>Ta then print "T0sq > Ta:  reject";
           else print "T0sq < Ta:  accept";
c1=s1[1,1]/n1/(s1[1,1]/n1+s2[1,1]/n2);
c2=s1[2,2]/n1/(s1[2,2]/n1+s2[2,2]/n2);
f1inv=c1**2/(n1-1)+(1-c1)**2/(n2-1);
f2inv=c2**2/(n1-1)+(1-c2)**2/(n2-1);
f1=round(1/f1inv,1);
f2=round(1/f2inv,1);
tf1a2k=tinv(1-alpha/(2*2),f1);
tf2a2k=tinv(1-alpha/(2*2),f2);
l1=(xb[1,1]-yb[1,1])-tf1a2k*sqrt(s1[1,1]/n1+s2[1,1]/n2);
u1=(xb[1,1]-yb[1,1])+tf1a2k*sqrt(s1[1,1]/n1+s2[1,1]/n2);
l2=(xb[1,2]-yb[1,2])-tf2a2k*sqrt(s1[2,2]/n1+s2[2,2]/n2);
u2=(xb[1,2]-yb[1,2])+tf2a2k*sqrt(s1[2,2]/n1+s2[2,2]/n2);
print "deltx1y1=("l1","u1"),",
      "deltx2y2=("l2","u2").";
quit;
          T0SQ            TA
      56.309901    3.5218933
          T0sq > Ta:  reject
 deltx1y1=( 0.2838606 , 32.593062 ),
 deltx2y2=( -50.33664 , -8.709514 ).