/* Appendix A */

   *** The macro SRS selects a random sample of size n     ;
   *** from dataset pop without replacement:               ;
   ***    pop: population data set                         ;
   ***    sample: resulting simple random sample           ;
   ***    n: sample size                                   ;
   *** *************************************************** ;
%macro SRS(pop,sample,n);
   %local i j;
   data &sample(drop=i j count);
      count=0;
      array sel_obs(&n) _temporary_;
      do i=1 to &n;
         redo: select=ceil(ranuni(0)*total);
         do j=1 to count;
            if sel_obs(j)=select then goto redo;
         end;
         set &pop point=select nobs=total;
         count=count+1;
         sel_obs(count)=select;
         output;
      end;
      stop;
   run;
%mend;


/* Appendix B */

   * Cluster and simple random population (POP_C)          ;
data pop_c;
   do obn=1 to 500;
      x=normal(12345);
      if mod(obn,50)<9 then
         cluster='ABC0'!!input(left(put(mod(obn,50)+1,2.)),$1.);
      else cluster='ABC'!!input(left(put(mod(obn,50)+1,2.)),$2.);
      output;
   end;
run;

proc means data=pop_c noprint;
   var x;
   output out=means mean=xmean std=stdx;

proc sort data=pop_c;
   by cluster;

data pop_c(keep=cluster x);
   if _n_=1 then set means;
   set pop_c;
   retain xmean stdx;
   x=(x-xmean)/stdx;
   x=100 + sqrt(50)*x;
run;


   * Stratification population (POP_S)                     ;
data pop_s;
   do obn=1 to 500;
      x=normal(12345);
      if obn<=100 then group='groupA';
      else if obn<=250 then group='groupB';
      else                  group='groupC';
      output;
   end;
run;

      * Summary statistics for each stratum.               ;

      * Group A                                            ;
   %let xmeana=100;
   %let xvara=25;

      * Group B                                            ;
   %let xmeanb=90;
   %let xvarb=20;

      * Group C                                            ;
   %let xmeanc=106;
   %let xvarc=30;

proc means data=pop_s nway noprint;
   class group;
   var x;
   output out=means mean=xmean std=stdx;

data pop_s(keep=group x);
   merge pop_s means;
   by group;
   x=(x-xmean)/stdx;

   if group='groupA' then x=&xmeana+sqrt(&xvara)*x;
   else if group='groupB' then x=&xmeanb+sqrt(&xvarb)*x;
   else if group='groupC' then x=&xmeanc+sqrt(&xvarc)*x;

run;