/* Example 17.10.1-----test Sigma=Sigma0 */
/* p=4, n=50, B=100, alpha=95%  Sigma0= */
libname sriv'/u/mengd/sasuser/ssm';
data sriv.tsq;
input tjstar;
cards;
;
run;
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};

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};

xs=s1`||s2`||s3`||s4`;
xsb=J(50,1,1)*J(1,50,1/50)*xs;
s=(xs-xsb)`*(xs-xsb)/49;
sigma0={0.12 0.10 0.01 0.01,
        0.10 0.14 0.01 0.01,
        0.01 0.01 0.03 0.01,
        0.01 0.01 0.01 0.01};
print sigma0;
T0=(50-1)*(-log(det(inv(sigma0)*s))+trace(inv(sigma0)*s-I(4)));
create sriv.t0 var{t0};
append from T0;
ys=(50/(50-1))**(1/2)*s**(-1/2)*(xs-xsb)`;
yi=ys`;
create sriv.yi var{y1 y2 y3 y4};
append from yi;
quit;
%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 bootstrap samples with replacement from yi */
proc iml;
use rands&n;
read all var{iid} into ids;
use sriv.yi;
read all into ys;
sumy=0;
create sriv.us var{y1 y2 y3 y4};
do i=1 to 50;
    ystar=ys[ids[i,],];
    append from ystar;
    sumy=sumy+ystar;
end;
ybstar=sumy/50;
/* Calculate the sample covariances */
use sriv.us;
read all into us;
sistar=(us`*us-50*ybstar`*ybstar)/(50-1);
/*  Calculate Tjstar */
Tjstar=(50-1)*(-log(det(sistar))+trace(sistar-I(4)));
/* Order the values of Tisq */
edit sriv.tsq ;
append from Tjstar;
quit;
%end;
proc iml;
use sriv.tsq;
read all into tjstar;
ord=rank(tjstar);
do i=1 to 100;
   rr=ord[i,];
   if rr=96 then
   T96=tjstar[i,];
end;
use sriv.t0;
read all into t0;
print T0 T96;
if T0>T96 then print "reject";
             else print "accept";
quit;
%mend ran;
%ran(k=100)