/* Example 17.4.1-----Confidence intervel for mean */
/* p=1, n=50, B=100, alpha=5%  */
libname sriv'/u/mengd/sasuser/ssm';
data sriv.zstar;
input zistar;
cards;
;
run;
%macro ran(k= );
%do n=1 %to &k;
data rands&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;
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;
use rands&n;
read all var{iid} into ids;
sum=0;
do i=1 to 50;
    xstar=x[ids[i,],];
    sum=sum+xstar;
end;
xbstar=sum/50;
zistar=sqrt(50)*(xbstar-xb);
edit sriv.zstar ;
append from zistar;
quit;
%end;

/*  Construct confidence interval */
proc iml;
use sriv.zstar;
read all into zistar;
ord=rank(zistar);
do i=1 to 100;
   rr=ord[i,];
   if rr=96 then
   u96=zistar[i,];
end;
use  sriv.xb;
read all into xb;
upper=xb+1/sqrt(50)*u96;
print u96 xb,
"upper limit of 95% confidence interval is",
       upper;
quit;
%mend ran;
%ran(k=100)