** A small IML simulation - Compute the distribution 
of the gap between the largest and next to largest 
values in a sample of size 50 from a N(mu, sigma_sq) 
distribution; 

Proc IML; 
*** set up parameters and dummy matrices ***; 
sigma = 10 ; mu = 100; 
Y = shape(0,1,10); data={0 0 0}; nextline = data; 
print Y data nextline; 
Y = shape(0,1,100); 

*** Create the data vector looping on rep ***; 
Do rep = 1 to 5000; 
Do i = 1 to 100; 
Y[1,i] = mu + sigma*Normal(123); end; 
rank=rank(Y); 

*** Extract highest and second highest ***; 
Do i = 1 to 100; 
  if rank[1,i]=99 then nextline[1,1]=Y[1,i]; 
  if rank[1,i]=100 then nextline[1,2]=Y[1,i]; 
end; nextline[1,3] = nextline[1,2]-nextline[1,1];
 
*** Show rep 1 ****; 
If rep=1 then do; 
   smallY= Y[1,1:5]; smallr= rank[1,1:5]; 
   data = nextline; 
   print smallY smallR data; end; 

*** Stack results *****; 
Else do; data=data//nextline; 

*** Show results after 10 reps ***; 
If rep = 10 then print data; 

end;  end; 

*** Create labels and data ***; 
labs = {X_99 X_100 gap}; 
create results from data [colname=labs]; 
append from data[colname=labs]; 

*** summarize results with PROCs ***; 
proc print data=results(obs=20); run; 
proc univariate normal plot; 
  var gap; 
proc gchart; vbar gap; 
proc gplot; plot gap*X_99; 
  symbol1 v=dot i=none c=red; 
run;