/*Example 10.4.1---Testing outliers */
proc iml;
x1={ 0,0,0,7,6,0,0,0,1,0,0,1,9,10,0,0,0,0,9,0,0,0,1,4,0,0,0,7,10};
x2={0,0,0,9,10,0,1,1,8,3,1,1,10,10,0,0,0,7,10,1,4,5,5,10,0,2,1,
    10,10};
x3={0,1,3,9,10,3,6,3,10,5,2,3,10,10,0,0,0,8,10,2,4,7,7,10,5, 6,7,
    10,10};
x4={1,3,5,9,10,7,7,7,10,5,2,3,10,10,0,0,3,8,10,4,6,8,7,10,6,7,10,
    10,10};
x5={1,3,5,9,10,7,8,9,10,5,3,3,10,10,0,0,4,8,10,7,8,10,7,10,6,7,10,
    10,10};
d={10,21,30,100,150,40,75,148,256, 27,21,30,90,150,16,43,49,140,
   290,21,
   27,30,30,100, 20,21,48, 80,145};
n=29; p=5; q=2; f=n-q+m-p-1; m=3;
y=arsin(sqrt((x1||x2||x3||x4||x5)/10))*180/3.14158;
yb=y`*J(n,n,1/n);
A=log10(d)||J(n,1,1);
B=J(1,5,1)//{1 2 3 4 5}//{1 4 9 16 25 };
H=A*inv(A`*A)*A`;
W=y`*(I(n)-H)*y;
e=y`-y`*H;
P=B`*inv(B*inv(W)*B`)*B;
Q=0;
do i=1 to n;
   Fi=(f-m+1)/m*e[,i]`*inv(W)*P*inv(W)*e[,i]/
      (1-H[i,i]-e[,i]`*inv(W)*e[,i]);
   if Fi>Q then do; k=i;Q=Fi;  end;
end;
Fa=finv(1-0.05/n,m,f-m+1);
print  k Q Fa ;
if Q>Fa then print "Q > Fa: reject.",
  "Observation "k " will be declared as an outlier";
         else print "Q < Fa: accept"; */
quit;
        k               Q          FA
       22       6.7512092   6.7504089
                 Q > Fa: reject
  Observation 22  may be declared as an outlier
  as it is a borderline case