/*Example 3.6.4 ---Test Multivariate Normality*/
/* p=4;n=28; alpha=0.05*/
libname sriv '/u/mengd/sasuser/ssm';
proc iml;
north={ 72 60 56 41 32 30 39 42 37 33 32 63 54 47
        91 56 79 81 78 46 39 32 60 35 39 50 43 48};
east={  66 53 57 29 32 35 39 43 40 29 30 45 46 51
        79 68 65 80 55 38 35 30 50 37 36 34 37 54};
south={ 76 66 64 36 35 34 31 31 31 27 34 74 60 52
       100 47 70 68 67 37 34 30 67 48 39 37 39 57};
west={  77 63 58 38 36 26 27 25 25 36 28 63 52 45
        75 50 61 58 60 38 37 32 54 39 31 40 50 43};
x=north`||east`||south`||west`;
n=28; p=4;
xb=J(1,n,1/n)*x;
xxb=J(n,1,1)*J(1,n,1/n)*x;
S=(x-xxb)`*(x-xxb)/(n-1);
call eigen(Du,H,S);
create sriv.du var{du};
append ;
create sriv.try1 var{y1i };
do i= 1 to n;
y1i=H[,1]`*x[i,]`;
append from y1i;
end;
create sriv.try2 var{y2i};
do i= 1 to n;
y2i=H[,2]`*x[i,]`;
append from y2i;
end;
create sriv.try3 var{y3i };
do i= 1 to n;
y3i=H[,3]`*x[i,]`;
append from y3i;
end;
create sriv.try4 var{y4i};
do i= 1 to n;
y4i=H[,4]`*x[i,]`;
append from y4i;
end;
quit;

proc  iml;
sort sriv.try1 by y1i;
sort sriv.try2 by y2i;
sort sriv.try3 by y3i;
sort sriv.try4 by y4i;
use sriv.try1;
read all into y1i;
use sriv.try2;
read all into y2i;
use sriv.try3;
read all into y3i;
use sriv.try4;
read all into y4i;
yij=y1i||y2i||y3i||y4i;
a={-.4328,-.2992,-.2510,-.2151,-.1857,-.1601,-.1372,
           -.1162,-.0965,-.0778,
 -.0598,-.0424,-.0253,-.0084,.0084,.0253,.0424,.0598,.0778,
           .0965,.1162,.1372,.1601,
 .1857,.2151,.2510,.2992,.4328};
ggama=-5.988; ddelta=1.919; eepsilon=0.1907;
use sriv.du;
read all into du;
n=28; p=4; M1=0; M2=10**5; alpha=0.05;
do i=1 to p;
wi=1/(n*du[i,])*(a`*yij[,i])**2;
Gi=ggama+ddelta*log((wi-eepsilon)/(1-wi));
M1=M1+(-2*log(probnorm(Gi)));
M2=min(M2,wi);
PHiGi=probnorm(Gi);
print "PHi(G("i"))="PHiGi;
end;
alpha0=1-exp(1/p*log(1-alpha));
chi= cinv(1-0.05, 2*p);
print M1 chi;
if M1>chi then print "M1>chi: reject";
         else print "M1<chi: accept";
print M2 alpha0;
quit;
              M1           CHI
       34.916149     15.507313
             M1>chi: reject

             M2          ALPHA0
      0.8629791       0.0127415
      PHi(G(2))= 0.0174384   PHi(G(4))= 0.0073349
      Since  PHi(G(4)) < ALPHA0 < PHi(G(2))
      implies     2 < malpha < 4
      Hence   M2 < 2 < malpha:  reject normality