/*Example 8.6.1---Canonical variates method */
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;
xbb=J(3,3,1)*xb[{1 11 21},]/3;
B=n1*(xb[{1 11 21},]-xbb)`*(xb[{1 11 21},]-xbb);
E=inv(sp)*B;
call geneig(lambda,a,B,sp);
print  lambda a  ;
quit;
       LAMBDA              A
    472.72109      -0.126179 -0.002696 0.0095824
    2.8962089      0.0571848 0.1327153 -0.058806
    -1.91E-14      0.3946508 -0.095081 0.3801019