/* Example 8.4.1.---Classifying into two unknown 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};
n1=10; n2=10;  f=n1+n2-2; p=3;
x1=xx1||xx2||xx3;
x2=yy1||yy2||yy3;
x1b=x1`*J(n1,n1,1)/n1;
x2b=x2`*J(n2,n2,1)/n2;
sp=((x1`-x1b)*(x1-x1b`)+(x2`-x2b)*(x2-x2b`))/f;
xb1=x1`*J(n1,1,1)/n1;
xb2=x2`*J(n2,1,1)/n2;
l=(xb1-xb2)`*inv(sp);
c=(xb1-xb2)`*inv(sp)*(xb1+xb2)/2;
deltsq=(f-p-1)/f*(xb1-xb2)`*inv(sp)*(xb1-xb2);
a1=(deltsq+12*(p-1))/(16*sqrt(deltsq))
    /sqrt(2*3.14159)*exp(-deltsq/2);
a2=(deltsq-4*(p-1))/(16*sqrt(deltsq))
    /sqrt(2*3.14159)*exp(-deltsq/2);
a3=sqrt(deltsq)*(p-1)/sqrt(2*3.14159)*exp(-deltsq/2)/4;
e1=1-probnorm(sqrt(deltsq)/2)+a1/n1+a2/n2+a3/f;
e2=1-probnorm(sqrt(deltsq)/2)+a1/n2+a2/n1+a3/f;
print  l ,c , e1 e2 ;
quit;
                   L
       -0.679225   0.4078591   2.7042804
                    C
             54.09784
              E1            E2
       0.0117693     0.0117693