/* Example 3.6.2.----Small's graphical method*/
/*p=4;n=27;*/
filename sriv1 'ch3p4.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);
create sriv.try var{ci};
do i=1 to n;
ci=n/(n-1)**2*(x[i,]-xb)*inv(s)*(x[i,]-xb)`;
append from ci;
end;
a=p/2; b=(n-p-1)/2;
alpha=(a-1)/(2*a); beta=(b-1)/(2*b);
create sriv.try1 var{di};
do i=1 to n;
di=betainv((i-alpha)/(n-alpha-beta+1),a, b);
append from di;
end;
quit;
proc iml;
sort sriv.try by ci;
use sriv.try;
read all into ci;
use sriv.try1;
read all into di;
cidi=ci||di;
create sriv.try2 var{ci di};
append;
quit;
proc gplot data=sriv.try2;
axis1 label=(h=1 'dI')
minor=none
major=(H=1pct N=4 W=1);
axis2 minor=none
major=(H=1pct N=5 W=1);
plot ci*di/frame haxis=axis1 vaxis=axis2;
symbol i=rl v=dot;
run;
quit;
goptions hsize=0cm vsize=0cm goutmode=replace;