/* Example 3.6.1.----Small's graphical method*/
filename sriv1 'ch3p2.ext';
goptions hsize=8cm vsize=8cm
device=pslepsf
gaccess=sasgaedt
gsfname=sriv1
gsfmode=replace;
libname sriv '/u/mengd/sasuser/ssm';
proc iml;
x1={ 59  80  17  11 134 76 50 85 139 175 55 30
    115 152 102  67  74 30  3 25  35  70 67 55
     31  57  54};
x2={ 40.9 42.0 42.7 306.0 146.5 101.0 21.0 22.7
     72.0 73.0  7.7  28.1 305.0  48.8 77.0 23.0
    302.0 48.8 30.1  16.2   8.7 152.5 36.4 151.5
    25.0 303.0 11.9};
x=x1`||x2`;
p=2; n=27;
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')
value=(h=0.3cm t=1 '0.0'
       t=2 '0.05'
       t=3 '0.1'
       t=4 '0.15'
       t=5 '0.20'
       t=6 '0.25')
minor=none
major=(H=1pct N=6 W=1);
axis2 label=(angle=90 h=1 'ci')
value=(a=90 h=0.3cm t=1 '0.0'
       t=2 '0.05'
       t=3 '0.1'
       t=4 '0.15'
       t=5 '0.20'
       t=6 '0.25')
minor=none
major=(H=1pct N=6 W=1);
plot ci*di /frame vaxis=axis2 haxis=axis1  ;
symbol  v=dot i=rl;
run;
quit;
goptions hsize=0cm vsize=0cm goutmode=replace;