/*Example 8.5.1---Classifying into k normals */
proc iml;
xx1={191,185,200,173,171,160,188,186,174,163};
xx2={131,134,137,127,128,118,134,129,131,115};
xx3={ 53, 50, 52, 50, 49, 47, 54, 51, 52, 47};
yy1={186,211,201,242,184,211,217,223,208,199};
yy2={107,122,144,131,108,118,122,127,125,124};
yy3={ 49, 49, 47, 54, 43, 51, 49, 51, 50, 46};
zz1={158,146,151,122,138,132,131,135,125,130};
zz2={141,119,130,113,121,115,127,123,119,120};
zz3={58,51,51,45,53,49,51,50,51,48};
n1=10; n2=10; n3=10; f=n1+n2+n3-3; p=3;
x=(xx1||xx2||xx3)//(yy1||yy2||yy3)//(zz1||zz2||zz3);
xb=(x`*(I(3)@J(n1,n1,1)/n1))`;
sp=((x-xb)[1:10,]`*(x-xb)[1:10,]
   +(x-xb)[11:20,]`*(x-xb)[11:20,]
   +(x-xb)[21:30,]`*(x-xb)[21:30,])/f;
l=xb[{1 11 21},]*inv(sp);
c=-eigval(diag(xb[{1 11 21},]*inv(sp)*xb[{1 11 21},]`))/2;
print l ,c ;
quit;
                 L
         -0.493395   1.0418004   5.8050957
         0.0917972   0.6885501    4.044093
         -1.134665   1.2471951   7.8778932
                            C
                    -198.6713
                    -169.2787
                    -150.7111