/* Example 7.2.2.---Profile analysis for two groups */
proc iml;
x1={11,12,16,11,6,12,8,19,11,10};
x2={16,11,15,14,12,8,11,8,9,12};
x3={23,18,22,19,16,20,14,9,21,21};
y1={12,6,12,21,9,8,7,14,10,21};
y2={13,10,18,21,10,18,13,20,20,24};
y3={24,13,22,24,22,24,24,24,22,24};
n1=10; n2=10; p=3; f1=n1-1; f2=n2-1;f=f1+f2;  alpha=0.05;
c={1 -1 0,0 1 -1};
x=x1||x2||x3;
xb=x`*J(n1,1,1/n1);
y=y1||y2||y3;
yb=y`*J(n2,1,1/n2);
u=xb-yb; l={1 1 1}`;
sp= ((x`-xb*J(1,n1,1))*(x`-xb*J(1,n1,1))`+
     (y`-yb*J(1,n2,1))*(y`-yb*J(1,n2,1))`)/f;
bsq=n1*n2/(n1+n2);
F0=(f-p+2)*bsq/(f*(p-1))*u`*c`*inv(c*sp*c`)*c*u;
Fa=finv(1-alpha, p-1,f-p+2);
print F0 Fa;
if F0 >Fa then print "F0 > Fa: reject H1";
          else print "F0 < Fa: accept H1";
Tsq=bsq*u`*c`*inv(c*sp*c`)*c*u;
F02=(f-p+2)/f*bsq*(l`*inv(sp)*u)**2*inv(l`*inv(sp)*l)/
    (1+1/f*Tsq);
Fa2=finv(1-alpha,1,f-p+1);
print F02 Fa2;
if F02 >Fa2 then print "F02 > Fa2: reject H2";
          else print "F02 < Fa2: accept H2";
z=(n1*xb+n2*yb)/(n1+n2);
F03=(n1+n2)*(f-p+3)/(p-1)*z`*c`*inv(c*f*sp*c`+bsq*c*u*u`*c`)*c*z;
Fa3=finv(1-alpha,p-1,n1+n2-p+1);
print F03 Fa3;
if F03 >Fa3 then print "F03 > Fa3: reject H3";
          else print "F03 < Fa3: accept H3";
gamah=l`*inv(sp)*u/(l`*inv(sp)*l);
e=1/sqrt(bsq)*tinv(1-alpha/2,f-p+1)*sqrt(1+Tsq/f)/
                sqrt(l`*inv(sp)*l)*sqrt(f/(f-p+1));
low=gamah-e;
upp=gamah+e;
print  low upp;
quit;
                     F0            FA
              2.8317278     3.5915306
                  F0 < Fa: accept H1
                    F02           FA2
              3.4749449     4.4939985
                 F02 < Fa2: accept H2
                    F03           FA3
              28.927701     3.5545571
                 F03 > Fa3: reject H3
                    LOW           UPP
              -6.649227     0.5271539