/*Example 7.3.1.----Profile analysis for J goups */
proc iml;
y11={41.4, 36.6, 37.0, 36.6, 40.9, 37.2, 40.2, 37.5, 39.2, 37.5};
y12={53.7, 45.0, 50.8, 48.2, 48.3, 51.1, 47.5, 41.8, 47.0, 44.7};
y13={65.9, 68.0, 69.3, 70.1, 66.0, 70.9, 66.0, 63.5, 70.5, 64.8};
y14={57.4, 52.3, 54.9, 54.1, 55.9, 53.3, 55.7, 52.4, 53.3, 53.8};

y21={23.7, 20.7, 20.3, 20.1, 19.0, 19.9, 20.3, 22.1, 21.1, 23.1,
             21.8, 19.3, 19.3, 22.9, 22.5};
y22={28.5, 26.1, 28.4, 30.3, 21.2, 27.8, 27.0, 30.4, 29.4, 31.8,
             30.8, 25.8, 26.4, 27.4, 32.6};
y23={48.9, 44.8, 47.5, 52.8, 47.1, 49.0, 49.7, 49.1, 47.5, 50.6,
             48.7, 50.1, 50.5, 48.0, 52.3};
y24={38.3, 34.8, 38.2, 37.0, 35.9, 36.5, 36.5, 37.9, 36.6, 39.1,
             37.9, 35.8, 35.2, 37.8, 38.2};
y31={32.5, 31.8, 30.8, 28.0, 32.2, 28.7, 34.0, 30.5, 29.3, 27.9,
             29.3, 33.4, 31.7, 27.1};
y32={39.8, 40.9, 37.7, 36.3, 40.6, 35.1, 41.0, 34.5, 37.5, 36.2,
             41.3, 36.4, 39.2, 36.2};
y33={55.7, 63.2, 57.8, 56.3, 62.1, 54.6, 61.2, 56.5, 58.6, 60.2,
             60.2, 55.4, 59.9, 54.8};
y34={47.9, 47.1, 46.4, 44.0, 47.3, 45.3, 48.2, 45.9, 46.0, 45.5,
             46.3, 48.2, 47.4, 44.2};
y41={18.7, 19.5, 18.8, 19.7, 13.4, 20.0, 18.9, 17.7, 23.1, 18.1,
             18.4, 16.0};
y42={20.3, 26.4, 29.2, 25.3, 20.4, 27.8, 26.6, 27.8, 26.7, 26.6,
             25.7, 20.5};
y43={44.1, 50.7, 48.0, 48.0, 42.3, 45.4, 50.5, 47.8, 47.6, 45.4,
             50.6, 43.4};
y44={33.8, 34.6, 36.0, 34.9, 30.0, 33.8, 35.4, 35.8, 36.2, 33.7,
             34.5, 32.5};
n1=10; n2=15; n3=14; n4=12; p=4; J=4; n=n1+n2+n3+n4;
f=n-J; alpha=0.05;
y1=y11||y12||y13||y14;
y2=y21||y22||y23||y24;
y3=y31||y32||y33||y34;
y4=y41||y42||y43||y44;
yy=(y1`||y2`||y3`||y4`);
yb=yy*J(n,1,1)/n;
ybj=(y1`*J(n1,1,1)/n1)||(y2`*J(n2,1,1)/n2)||
    (y3`*J(n3,1,1)/n3)||(y4`*J(n4,1,1)/n4);
sst=yy*yy`-n*yb*yb`;
ybs=(y1`*J(n1,1,1)/sqrt(n1))||(y2`*J(n2,1,1)/sqrt(n2))||
    (y3`*J(n3,1,1)/sqrt(n3))||(y4`*J(n4,1,1)/sqrt(n4));
sstr=ybs*ybs`-n*yb*yb`;
sse=sst-sstr;
c={1 -1 0 0 , 1 -2 1 0, 1 -1 1 -1};
lambda1=det(c*sse*c`)/det(c*(sse+sstr)*c`);
lnlam=-(n-J-(p+1-J)/2)*log(lambda1);
chia1=cinv(1-alpha,(p-1)*(J-1));
print lnlam  chia1;
if lnlam> chia1 then  print "lnlam > chia1: reject H1";
                  else  print "lnlam < chia1: accept H1";
lambda2=det(sse)/det(sse+sstr)/lambda1;
ratiolam=(n-J-p+1)/(J-1)*(1-lambda2)/lambda2;
Fa2=finv(1-alpha, J-1, n-j-p+1);
print ratiolam  Fa2;
if ratiolam> Fa2 then  print "ratiolam > Fa2: reject H2";
                else  print "ratiolam < Fa2: accept H2";
F0=n*(n-p+1)/(p-1)*yb`*c`*inv(c*sse*c`+c*sstr*c`)*c*yb;
Fa3=finv(1-alpha, p-1,n-p+1);
print   F0 Fa3;
if F0> Fa3 then  print "F0 > Fa3: reject H3";
                else  print "F0 < Fa3: accept H3";
ybjd=y4`*J(n4,1,1)/n4*J(1,4,1);
z=(ybj-ybjd)[,1:3];
gamah=z`*inv(sse)*J(4,1,1)/(J(1,4,1)*inv(sse)*J(4,1,1));
AA={0.183 0.083 0.083,
    0.083 0.15 0.083,
    0.083 0.083 0.155};
Tasq=(J-1)/(f-p+1)*finv(1-0.05,J-1,f-p+1);
Tlsq=tinv(1-0.05/6,f-p+1)**2;
/*since Tlsq=6.195 > 0.912=Tasq we use Tasq*/
a1={1, 0, 0};  a2={0, 1, 0}; a3={0, 0, 1};
e1=sqrt(Tasq)/sqrt(J(1,4,1)*inv(sse)*J(4,1,1))*
    sqrt(a1`*(AA+z`*c`*inv(c*sse*c`)*c*z)*a1);
e2=sqrt(Tasq)/sqrt(J(1,4,1)*inv(sse)*J(4,1,1))*
    sqrt(a2`*(AA+z`*c`*inv(c*sse*c`)*c*z)*a2);
e3=sqrt(Tasq)/sqrt(J(1,4,1)*inv(sse)*J(4,1,1))*
    sqrt(a3`*(AA+z`*c`*inv(c*sse*c`)*c*z)*a3);
low1=a1`*gamah-e1;
upp1=a1`*gamah+e1;
low2=a2`*gamah-e2;
upp2=a2`*gamah+e2;
low3=a3`*gamah-e3;
upp3=a3`*gamah+e3;
print "(" low1"," upp1")",
      "(" low2"," upp2")",
      "(" low3"," upp3")";
quit;
               LNLAM        CHIA1
            9.973048    16.918978
          lnlam < chia1: accept H1
            RATIOLAM          FA2
           414.20066    2.8164658
         ratiolam > Fa2: reject H2
                  F0          FA3
           4327.7343    2.7980606
              F0 > Fa3: reject H3
                LOW1       UPP1
           (17.61664,  21.29209)
                LOW2       UPP2
           ( 1.04051    4.16226)
                LOW3       UPP3
           (10.18850   13.39743)