/* Example 17.11.1-----test Sigma1=Sigma2 */
/* p=4, n=50, B=100, alpha=95%   */
libname sriv'/u/mengd/sasuser/ssm';
data sriv.tstar;
input tistar;
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};

v1={7.0 6.4 6.9 5.5 6.5 5.7 6.3 4.9 6.6 5.2 5.0 5.9 6.0 6.1 5.6 6.7
5.6 5.8 6.2 5.6 5.9 6.1 6.3 6.1 6.4 6.6 6.8 6.7 6.0 5.7 5.5 5.5 5.8
6.0 5.4 6.0 6.7 6.3 5.6 5.5 5.5 6.1 5.8 5.0 5.6 5.7 5.7 6.2 5.1 5.7};

v2={3.2 3.2 3.1 2.3 2.8 2.8 3.3 2.4 2.9 2.7 2.0 3.0 2.2 2.9 2.9 3.1
3.0 2.7 2.2 2.5 3.2 2.8 2.5 2.8 2.9 3.0 2.8 3.0 2.9 2.6 2.4 2.4 2.7
2.7 3.0 3.4 3.1 2.3 3.0 2.5 2.6 3.0 2.6 2.3 2.7 3.0 2.9 2.9 2.5 2.8};

v3={4.7 4.5 4.9 4.0 4.6 4.5 4.7 3.3 4.6 3.9 3.5 4.2 4.0 4.7 3.6 4.4
4.5 4.1 4.5 3.9 4.8 4.0 4.9 4.7 4.3 4.4 4.8 5.0 4.5 3.5 3.8 3.7 3.9
5.1 4.5 4.5 4.7 4.4 4.1 4.0 4.4 4.6 4.0 3.3 4.2 4.2 4.2 4.3 3.0 4.1};

v4={1.4 1.5 1.5 1.3 1.5 1.3 1.6 1.0 1.3 1.4 1.0 1.5 1.0 1.4 1.3 1.4
1.5 1.0 1.5 1.1 1.8 1.3 1.5 1.2 1.3 1.4 1.4 1.7 1.5 1.0 1.1 1.0 1.2
1.6 1.5 1.6 1.5 1.3 1.3 1.3 1.2 1.4 1.2 1.0 1.3 1.2 1.3 1.3 1.1 1.3};

xs=s1`||s2`||s3`||s4`;
xv=v1`||v2`||v3`||v4`;
xsb=J(50,1,1)*J(1,50,1/50)*xs;
s1=(xs-xsb)`*(xs-xsb)/49;
xvb=J(50,1,1)*J(1,50,1/50)*xv;
s2=(xv-xvb)`*(xv-xvb)/49;
xxbs=J(1,50,1/50)*xs;
xxbv=J(1,50,1/50)*xv;
f=50+50-2; f1=50-1; f2=50-1;
s=(f1*s1+f2*s2)/f;
T0=(det(s1)**(f1/2)*det(s2)**(f2/2))/det(s)**(f/2);
create sriv.t0 var{t0};
append from T0;
us=s**(1/2)*s1**(-1/2)*(xs-xsb)`;
vv=s**(1/2)*s2**(-1/2)*(xv-xvb)`;
ui=us`; vi=vv`;
create sriv.ui var{y1 y2 y3 y4};
append from ui;
create sriv.vi var{y1 y2 y3 y4};
append from vi;
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;
data randv&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;
/* Draw a bootstrap sample with replacement from ui and vi */
proc iml;
use rands&n;
read all var{iid} into ids;
use sriv.ui;
read all into us;
use randv&n;
read all var{iid} into idv;
use sriv.vi;
read all into vs;
sumu=0;
create sriv.us var{y1 y2 y3 y4};
do i=1 to 50;
    ustar=us[ids[i,],];
    append from ustar;
    sumu=sumu+ustar;
end;
ubstar=sumu/50;
sumv=0;
create sriv.vv var{y1 y2 y3 y4};
do i=1 to 50;
    vstar=vs[idv[i,],];
    append from vstar;
    sumv=sumv+vstar;
end;
vbstar=sumv/50;
/* Calculate the sample covariances */
use sriv.us;
read all into us;
use sriv.vv;
read all into vv;
si1star=(us`*us-50*ubstar`*ubstar)/(50-1);
si2star=(vv`*vv-50*vbstar`*vbstar)/(50-1);
f=50+50-2; f1=50-1; f2=50-1;
sistar=(f1*si1star+f2*si2star)/f;
/* Calculate Tistar  and order the values of Tisq */
Tistar=(det(si1star)**(f1/2)*det(si2star)**(f2/2))
/det(sistar)**(f/2);
edit sriv.tstar ;
append from Tistar;
quit;
%end;
proc iml;
use sriv.tstar;
read all into tistar;
ord=rank(tistar);
do i=1 to 100;
   rr=ord[i,];
   if rr=96 then
   T96=tistar[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)