/* Example 17.8.1-----test u1=u2 */
/* p=4, n=50, B=100, alpha=95%   */
libname sriv'/u/mengd/sasuser/ssm';
data sriv.tsq;
input tisq;
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;
sp=((50-1)*s1+(50-1)*s2)/f;
T0sq=(f-4+1)*(50*50)/(f*4*(50+50))*(xxbs-xxbv)*inv(sp)*
(xxbs-xxbv)`;
create sriv.t0sq var{t0sq};
append from T0sq;
us=(50/(50-1))**(1/2)*sp**(1/2)*s1**(-1/2)*(xs-xsb)`;
t=us`;
create sriv.uus var{y1 y2 y3 y4};
append from t;
vv=(50/(50-1))**(1/2)*sp**(1/2)*s2**(-1/2)*(xv-xvb)`;
tt=vv`;
create sriv.vvv var{y1 y2 y3 y4};
append from tt;
  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 bootstrap samples from us and vv*/
proc iml;
use rands&n;
read all var{iid} into ids;
use randv&n;
read all var{iid} into idv;
use sriv.uus;
read all into us;
use sriv.vvv;
read all into vv;
create sriv.us var{y1 y2 y3 y4};
do i=1 to 50;
    ustar=us[ids[i,],];
    append from ustar;
end;
create sriv.vv var{y1 y2 y3 y4};
do i=1 to 50;
    vstar=vv[idv[i,],];
    append from vstar;
end;
/* Calculate the pooled covariance */
use sriv.us;
read all into us;
use sriv.vv;
read all into vv;
ubstar=J(50,1,1)*J(1,50,1/50)*us;
vbstar=J(50,1,1)*J(1,50,1/50)*vv;
s1star=(us-ubstar)`*(us-ubstar)/49;
s2star=(vv-vbstar)`*(vv-vbstar)/49;
ub=J(1,50,1/50)*us;
vb=J(1,50,1/50)*vv;
spstar=((50-1)*s1star+(50-1)*s2star)/(50+50-2);
/* Calculate T*^2 */
Tisq=(98-4+1)*50*50/((50+50)*98*4)*(ub-vb)*
inv(spstar)*(ub-vb)`;
edit sriv.tsq;
append from Tisq;
quit;
%end;
proc iml;
use sriv.tsq;
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 into t0sq;
print T0sq T96;
if T0sq>T96 then print "reject";
             else print "accept";
quit;
%mend ran;
%ran(k=100)