datasamplingFrom;
  call streaminit(12345);
  do obs =1 to 100000000;
    s1=rand('table',0.50);
    s2=rand('table',0.30);
    s3=rand('table',0.40);
    s4=rand('table',0.30,0.20,0.20);
   output;
  end;
run;
data samplingProps;
   length list $640.;
   do s1 =1 to 2;
     do s2 =1 to 2;
        do s3 =1 to 2;
           do s4 =1 to 4;
              _alloc_=  ifn(s1 =1,3/6,1-3/6)
                       *ifn(s2 =1,2/3,1-2/3)
                       *ifn(s3 =1,1/2,1-1/2)
                       *ifn(s4 in(14),2/8,ifn(s4 =3,3/8,1/8));
              list =catx(',',list,_alloc_);
              p ++_alloc_;
              output;
           end;
        end;
     end;
   end;
   *Maybe 1% tolerance;
   if fuzz(p)>=1.01 | fuzz(p) <=0.99 thenput"ERROR: Sum of props="p"(^=1),check the proportions.";
  callsymputx('list',list);
run;
%put &list;
%put ///sum of props=:%sysfunc(sum(&list));
data SamplingOut;
  call streaminit(12345);
  totSamples =1000;
  array t[2,2,2,4,2] _temporary_;
  if _n_ =1 then do;
    do p1 =1 tonobs1;
     set SamplingFrom point =p1 nobs =nobs1;
     t[s1,s2,s3,s4,1] ++1;
    end;
    do p2 =1 tonobs2;
     set SamplingProps point =p2 nobs =nobs2;
     t[s1,s2,s3,s4,2] =round(_alloc_*totSamples);
     totS ++t[s1,s2,s3,s4,2];
    end;
    puttotS;
    *adjusttotal sampling size when rounding problem occurs;
    *randomlydistribute rounding discrepancy to all strata;
    if totS^=totSamples then do i =1 toabs(totS-totSamples);   t[rand('table',.5),rand('table',.5),rand('table',.5),rand('table',.25,.25,.25),2]+(1*sign(totSamples-tots));
    end;
  end;
  set SamplingFrom;
  if t[s1,s2,s3,s4,2]>0 then ifrand('uniform') <= t[s1,s2,s3,s4,2]/t[s1,s2,s3,s4,1]then do;
   output;   
   t[s1,s2,s3,s4,2] +-1;
  end;