Just comment on your friend on the blog, the datastep program may be simpler if using array,

data bites;;
if _n_ le 18 then drink="water";
else drink="beer";
infile datalines delimiter=',';
input bites @@;
datalines;
21, 19, 13, 22, 15, 22, 15, 22, 20, 12, 24, 24, 21, 19, 18, 16, 23, 20
27, 19, 20, 20, 23, 17, 21, 24, 31, 26, 28, 20, 27, 19, 25, 31, 24
28, 24, 29, 21, 21, 18, 27, 20
;

proc sql noprint;
select bites into : data separated by "," from bites;
select count(bites) into : numobs from bites;
select count(bites) into : water from bites where drink="water";
select count (bites) into : beer from bites where drink="beer";
run;

data temp;
array x[&numobs] _TEMPORARY_ (&data);
array y1[&water] _temporary_;
array y2[&beer] _temporary_;
seed=3987674;

do j=1 to dim(y1);
y1[j]=x[j];
end;

do j=dim(y1)+1 to dim(x);
y2[j-dim(y1)]=x[j];
end;

mainmean=mean(of y2[*])-mean(of y1[*]);
call symput('mainmean', mainmean);
do i=1 to 100000;
call ranperm(seed, of x[*]);
do j=1 to dim(y1);
y1[j]=x[j];
end;
do j=dim(y1)+1 to dim(y2);
y2[j-dim(y1)]=x[j];
end;
result=mean(of y2[*])-mean(of y1[*]);
output;
end;
run;

proc sgplot data=temp;
histogram result ;
refline &mainmean/axis=x;
run;