proc iml;
start myQntl(q, x, p); /* definition 5 from UNIVARIATE doc */
y = colvec(x);
call sort(y,1);
n = nrow(y); /* assume nonmissing data */
q = j(ncol(p),1);
do i = 1 to ncol(p);
j = (n+1)*p; /* position in ordered data */
j1 = floor(j); j2 = ceil(j);/* indices into ordered data */
print j j1 j2;
if j1=j2 then /* return a datum */
q = y[j1];
else /* interpolate between data */
q = (y[j2]+y[j1])/2;
end;
finish;
/* test it */
x = ranuni(j(100,1));
p = {0.1 0.25 0.5 0.75 0.9};
call MyQntl(q, x, p);
print q;