/* PERCTILE calculates percentile values
for a given series of n real values. It does
the interpolation between points as you would
expect [e.g., in  the median when n is even,
reporting (x[n/2] + x[n/2+1])/2] rather than just
taking the closest point.
Example:
  n=100;
  pc = perctile(rndu(n,5), 0.25|0.5|0.75);
  Calculates the quartiles for each of the
  5 columns of 100 random numbers.

The code is provided for public non-commercial
use, without performance guarantees.
Posted April 1996 by Dr. David Baird.
____________________________________________________________________
Dr David Baird, Biometrician        EMail:  BairdD@AgResearch.CRI.NZ
Mail: AgResearch, PO Box 60, Gerald St, Lincoln, NEW ZEALAND
Phone: +64 3 325 6900   Fax: +64 3 325 2946
*/
/* PERCTILE: Estimates the percentiles of the columns of a matrix

Usage:  pc = perctile(y,p);

Input:
      y - R x C matrix of data
      p - 1 x N or N x 1 matrix of proportions (0 =< x =< 1)

Output:
     pc - a N x C matrix with each column giving the
          percentiles corresponding  to p for each column of y

Example:
     pc = perctile(rndu(100,5), 0.25|0.5|0.75);
     Calculates the quartiles of the 5 columns of random numbers. */

proc perctile(y,p) ;
     local i,n,x,k,r,np,ny,pc,z ;
     n = rows(y) ;
     x = 1 + p.*(n-1) ;
     k = floor(x) ;
     r = x - k ;
     np = rows(p) ;
     ny = cols(y) ;
     pc = zeros(np,ny) ;
     i = 1 ;
     do until(i > ny) ;
        z = sortc(y[.,i],1) ;
        pc[.,i] = (1-r).*z[k] + r.*z[minc((k'+1)|(n.*ones(1,np)))] ;
        i = i + 1 ;
     endo ;
     retp(pc) ;
endp ;