/* Example 3.2.1-----detecting an outlier */
/* p=2, n=32, alpha=5% */
proc iml;
x={26 40 32 36 27 39 29 22 34 25 37 36 20
   26 31 22 25 29 24 34 28 21 40 35 24 35
   31 21 36 36 35 24};
y={35.0 57.0 34.3 38.3 37.5 36.3 31.3 35.0
   30.0 30.0 40.5 47.5 30.0 40.0 30.0 26.7
   34.0 32.5 28.3 40.8 29.3 26.0 47.0 40.0
   30.0 45.0 37.5 25.0 51.1 42.5 45.0 25.0};
z=x`||y`;
zb=J(32,1,1)*J(1,32,1/32)*z;
s=(z-zb)`*(z-zb)/31;
D=(z-zb)*inv(s)*(z-zb)`;
Q=max(vecdiag(d));
p=vecdiag(d)`[,<:>];
ca=2/(32-2-1)*FINV(0.05/32,2,29);
F=ca/(1+ca)*31**2/32;
print F Q p;
if Q >F then print
" N0. "p" observation is an outlier";
quit;
                 F            Q         P
         0.0032385    7.2637302         2
        No. 2  observation is an outlier