quantiles with standard errors

[Francesco Potorti`]

Disclaimer: I am not a numerical analyst.  The algorithms used here may
be suboptimal AND EVEN WRONG.  In fact, I have checked it with real
and simulated data and it appears to work correctly, but you are warned.

This function is a sort of wrapper around 'quantile' with the limitation
that it only accepts a vector of samples and the difference that it
always returns a row vector.  These limitations are easily removed.
Little error checking is done, but this is easy to add.

The second value returned is the standard error of the computed quantile
values.  For 1000 samples it is about 50% slower than 'quantile', for
1e6 samples it has the same speed, so standard errors are computed for
free.  From standard errors one can easily get confidence interval
estimates for the quantile values.  Remember disclaimer and enjoy:


## Compute quantiles with their standard errors
## Return row vectors QVAL containing the quantile values
## and QERR containing their standard error.
function [qval qerr] = quantilerr (x, qlev)
  if (!isvector(x) || !isvector(qlev))
    error("function dquant");
  endif
  n = length(x);                     # number of samples
  ## Define intervals for computing PDF at qval points
  intvn = sqrt(n);                   # number of samples per interval
  intvw = intvn/n;                   # interval width
  qlev = qlev(:)';                   # make it a row vector
  qq = [qlev-intvw/2; qlev; qlev+intvw/2];
  quantiles = reshape(quantile(x,qq(:)), 3, []);
  qval = quantiles(2,:);        # quantile values
  pdf = intvw./diff(quantiles([1 3],:));
  qleverr = sqrt(qlev.*(1-qlev)/n);  # standard error on q levels
  qerr = qleverr./pdf;               # standard error on q values
endfunction

-- 
Francesco Potortì (ricercatore)        Voice: +39 050 315 3058 (op.2111)
ISTI - Area della ricerca CNR          Fax:   +39 050 315 2040
via G. Moruzzi 1, I-56124 Pisa         Email: address@hidden
(entrance 20, 1st floor, room C71)     Web:   http://fly.isti.cnr.it/