/*Example 3.6.3 ---Srivastava's graphical method*/
/* p=4;n=28;*/
filename sriv1 'ch3p51.ext';
goptions hsize=8cm vsize=8cm
device=pslepsf
gaccess=sasgaedt
gsfname=sriv1
gsfmode=replace;
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);
ds=inv(sqrt(diag(s)));
R=ds*S*ds;
call eigen(D,L,R);
z=x*ds*L;
z1=z[,1];z2=z[,2];z3=z[,3];z4=z[,4];
create sriv.try1 var{z1};
append from z1 ;
create sriv.try2 var{z2};
append from z2 ;
create sriv.try3 var{z3};
append from z3 ;
create sriv.try4 var{z4};
append from z4 ;
create sriv.try5 var{xi};
do i= 1 to n;
xi=probit((i-3/8)/(n+1/4));
append from xi;
end;
quit;
proc  iml;
sort sriv.try1 by z1;
use sriv.try1;
read all into z1;
sort sriv.try2 by z2;
use sriv.try2 ;
read all  into  z2;
sort sriv.try3 by z3;
use sriv.try3 ;
read all into z3;
sort sriv.try4 by z4;
use sriv.try4;
read all  into  z4;
use sriv.try5;
read all into xi;
zx=z1||z2||z3||z4||xi;
create sriv.try6 var{z1i z2i z3i z4i xi} ;
append from zx;
quit;
axis1 minor=none;
axis2 minor=none;
symbol i=rl v=dot;
proc gplot data=sriv.try6;
plot z1i*xi/frame haxis=axis1 vaxis=axis2;
run;quit;
filename sriv1 'ch3p52.ext';
proc gplot data=sriv.try6;
plot z2i*xi/frame haxis=axis1 vaxis=axis2;
run;quit;
filename sriv1 'ch3p53.ext';
proc gplot data=sriv.try6;
plot z3i*xi/frame haxis=axis1 vaxis=axis3;
axis3 minor=none major=(n=5 w=1);
run;quit;
filename sriv1 'ch3p54.ext';
proc gplot data=sriv.try6;
plot z4i*xi/frame haxis=axis1 vaxis=axis2;
run;quit;