/* Example 8.3.1---Classifying two normals with sigma1=sigma2*/
proc iml;
sigma={2.3008 .2516 .4742, .2516 .6075 .0358, .4742 .0358 .5951};
xb={2.9298 1.6670 .7281,    3.0303 1.2424 .5455};
p=3;   n1=114; n2=33;
l=(xb[1,]-xb[2,])*inv(sigma);
c=(xb[1,]-xb[2,])*inv(sigma)*(xb[1,]+xb[2,])`/2;
deltsq=(xb[1,]-xb[2,])*inv(sigma)*(xb[1,]-xb[2,])`;
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);
e1=1-probnorm(sqrt(deltsq)/2)+a1/n1+a2/n2;
e2=1-probnorm(sqrt(deltsq)/2)+a1/n2+a2/n1;
print  l, c ,  e1 e2 ;
 quit;
                          L
                  -0.216447  0.7630322  0.4334109
                                     C
                             0.7409547
                          E1            E2
                   0.3717719     0.3930834