/* Example 17.7.1----test u=u0 */
libname sriv'/u/mengd/sasuser/ssm';
/* p=4, n=50, B=100, alpha=0.95 */
/* u0={ 5.0 3.42 1.4 0.24}  */
data sriv.tisq;
input tisq;
cards;
;
run;
%macro ran(k= );
%do n=1 %to &k;
data rand&n;
retain seed &n&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;
proc iml;
use rand&n;
read all var{iid} into id;
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};

s2={3.5 3.0 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 3.7 3.4 3.0 3.0 4.0 4.4
3.9 3.5 3.8 3.8 3.4 3.7 3.6 3.3 3.4 3.0 3.4 3.5 3.4 3.2 3.1 3.4 4.1
4.2 3.1 3.2 3.5 3.6 3.0 3.4 3.5 2.3 3.2 3.5 3.8 3.0 3.8 3.2 3.7 3.3};

s3={1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 1.5 1.6 1.4 1.1 1.2 1.5
1.3 1.4 1.7 1.5 1.7 1.5 1.0 1.7 1.9 1.6 1.6 1.5 1.4 1.6 1.6 1.5 1.5
1.4 1.5 1.2 1.3 1.4 1.3 1.5 1.3 1.3 1.3 1.6 1.9 1.4 1.6 1.4 1.5 1.4};

s4={0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 0.2 0.2 0.1 0.1 0.2 0.4
0.4 0.3 0.3 0.3 0.2 0.4 0.2 0.5 0.2 0.2 0.4 0.2 0.2 0.2 0.2 0.4 0.1
0.2 0.2 0.2 0.2 0.1 0.2 0.2 0.3 0.3 0.2 0.6 0.4 0.3 0.2 0.2 0.2 0.2};

x=s1`||s2`||s3`||s4`;
xb=J(50,1,1)*J(1,50,1/50)*x;
xxb=J(1,50,1/50)*x;
u0={5.0 3.42 1.4 0.24};
s=(x-xb)`*(x-xb)/49;
T0sq=50*(xxb-u0)*inv(s)*(xxb-u0)`;
create sriv.t0sq var{t0sq};
append;
y0=(50/(50-1))**(1/2)*(x-xb);
/* Draw a bootstrap samples with replacement from y0*/
create sriv.y1 var{y1 y2 y3 y4};
 do i=1 to 50;
      ystar=y0[id[i,],];
      append from ystar;
end;
/* Calculate the mean and covariance from the bootstrap sample */
use sriv.y1;
read all into yy;
yib=J(50,1,1)*J(1,50,1/50)*yy;
si=(yy-yib)`*(yy-yib)/49;
Tisq=50*yib[1,]*inv(si)*yib[1,]`;
edit sriv.tisq;
append from Tisq;
quit;
%end;
proc iml;
use sriv.tisq;
read all into tisq;
ord=rank(tisq);
do i=1 to 100;
   rr=ord[i,];
   if rr=96 then
             T96=tisq[i,];
end;
use sriv.t0sq;
read all var{t0sq} into t0sq;
print T0sq T96;
if T0sq>T96 then print "reject";
             else print "accept";
quit;
%mend ran;
%ran(k=100)