/* Example 17.3.1-----Confidence intervel for mean */
/* p=1, n=50, B=100, alpha=90%  */
libname sriv'/u/mengd/sasuser/ssm';
data sriv.u;
input ui;
cards;
;
run;
Proc iml;
/* Define: s1=serosa sepal length */
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;
J=J(1,50, 1/49);
ssq=J*((x-xb)#(x-xb));
create sriv.xbssq var{xb ssq};
append;
create sriv.y var{yi};
do i=1 to 50;
   yi=(50/(50-1))**(1/2)*(x[i,]-xb);
   append from yi;
end;
quit;
%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;
/*  Draw a bootstrap samples with replacement from y */
proc iml;
use rands&n;
read all var{iid} into ids;
use sriv.y;
read all into y;
sumy=0;
do i=1 to 50;
      ystar=y[ids[i,],];
      sumy=sumy+ystar;
end;
ybstar=sumy/50;
sums=0;
do i=1 to 50;
     sums=sums+(y[ids[i,],]-ybstar)**2;
end;
sistarsq=sums/(50-1);
ui=sqrt(50)*ybstar/sqrt(sistarsq);
edit sriv.u ;
append from ui;
quit;
%end;
/*  Construct confidence interval */
proc iml;
use sriv.u;
read all into ui;
ord=rank(ui);
do i=1 to 100;
   rr=ord[i,];
   if rr=5 then u5=ui[i,];
           else if rr=96 then u96=ui[i,];
end;
use  sriv.xbssq;
read all into xbssq;
xb=xbssq[1,1];
ssq=xbssq[1,2];
upper=xb+sqrt(ssq)/sqrt(50)*u96;
lower=xb+sqrt(ssq)/sqrt(50)*u5;
print xb ssq, u5 u96,
"90% confidence interval is", lower upper ;
quit;
%mend ran;
%ran(k=100)