/* Example 17.6.1-----test sigma=sigma0 */
libname sriv'/u/mengd/sasuser/ssm';
/* p=1, n=50, B=100, alpha=90%  sigma0=0.115*/
/* Chi-square test */
data sriv.t;
input tistarsq;
cards;
;
run;
%macro ran(k= );
%do n=1 %to &k;
data rands&n;
retain seed &n&n;
do i=1 to 50;
call rantbl(seed,
 1/50, 1/50, 1/50, 1/50, 1/50, 1/50, 1/50, 1/50, 1/50, 1/50,
 1/50, 1/50, 1/50, 1/50, 1/50, 1/50, 1/50, 1/50, 1/50, 1/50,
 1/50, 1/50, 1/50, 1/50, 1/50, 1/50, 1/50, 1/50, 1/50, 1/50,
 1/50, 1/50, 1/50, 1/50, 1/50, 1/50, 1/50, 1/50, 1/50, 1/50,
 1/50, 1/50, 1/50, 1/50, 1/50, 1/50, 1/50, 1/50, 1/50, 1/50,iid);
output;
 end;
run;
/* Draw a bootstrap samples with replacement from yi */
proc iml;
s1={5.1 4.9 4.7 4.6 5.0 5.4 4.6 5.0 4.4 4.9 5.4 4.8 4.8 4.3 5.8 
    5.7 5.4 5.1 5.7 5.1 5.4 5.1 4.6 5.1 4.8 5.0 5.0 5.2 5.2 4.7 
    4.8 5.4 5.2 5.5 4.9 5.0 5.5 4.9 4.4 5.1 5.0 4.5 4.4 5.0 5.1
    4.8 5.1 4.6 5.3 5.0};
x=s1`;
J=J(1,50,1/50);
xb=J*x;
J=J(1,50, 1/49);
ssq=J*((x-xb)#(x-xb));
sigma0=0.115;
T0=ssq*(50-1)/sigma0;
create sriv.t0 var{ssq sigma0 T0};
append;
/* Define the residuals vectors */
xbs=j(50,1,xb);
yi=(50/(50-1))**(1/2)*(x-xbs)/ssq**(1/2);
use rands&n;
read all var{iid} into ids;
sumy=0;
create sriv.ys var{y1};
do i=1 to 50;
    ystar=yi[ids[i,],];
    append from ystar;
    sumy=sumy+ystar;
end;
ybstar=sumy/50;
use sriv.ys;
read all into y1;
Tistarsq=y1`*y1-50*ybstar**2;
edit sriv.t ;
append from Tistarsq;
quit;
%end;
proc iml;
use sriv.t;
read all into tistarsq;
ord=rank(tistarsq);
do i=1 to 100;
   rr=ord[i,];
   if rr=5 then
   T5=tistarsq[i,];
   else if rr=96 then
   T96=tistarsq[i,];
end;
/* test */
use sriv.t0;
read all into t00;
ssq=t00[1,1];
sigma0=t00[1,2];
T0=t00[1,3];
print sigma0 ssq;
print T0 T5 T96;
if T5<T0 & T0<T96 then print "accept";
             else print "reject";
quit;
%mend ran;
%ran(k=100)