//  QQ_NORMAL
//  Input: e = unsorted error residual vector to be checked.
//  Output: eq = sample quantile (sorted e)
//          nq = theoretical N(0,1) quantile.
//  Plot (nq,eq) as (x,y) pairs to make a Q-Q plot.
//  
//    ^ Usage: ^
//  [nq,eq]=QQ_NORMAL(e);
//  plot(nq,eq,'o')
//  xlabel('Theoretical Normal N(0,1) Quantile')
//  ylabel('Sorted data values (user units)')

//   ^ Uses Matlab function "inverf", inverse error fn. ^
//   Tests the elements of vector e to examine whether they
//   are distributed according to a normal distribution.  The
//   output of this function may be used to construct a q-q
//   plot, for visual inspection of the distribution. If a normal
//   distribution is appropriate, the resulting locus of (nq,eq)
//   pairs will approximately follow a straight line, whose
//   intercept and slope respectively correspond to the mean
//   and standard deviation of the normal distribution.
//   Outliers, if present, appear as departures from linearity,
//   usually away from the origin, and are easily identified.
//   Keywords: Outliers, Robust Statistics, Data analysis

//   Copyright (c) 1989 by Peter R. Shaw
//   Permission to copy all or part of this work is granted, provided
//   that the copies are not distributed for resale, and that the
//   copyright notice and this notice are retained.
//  
//   If you make any improvements, please let me know.

//   Peter R. Shaw
//   Woods Hole Oceanographic Institution
//   Woods Hole, MA 02543
//   (508) 457-2000 ext. 2473
//   pshaw@aqua.whoi.edu

//   Translated to RLaB 5/10/93, Ian Searle

//   Reference:
//   Lewis, T. and N. I. Fisher, Graphical methods for investigating
//   the fit of a Fisher distribution to spherical data; Geophys. J.
//   R. Astr. Soc., 69:1-13, 1987.

require erf

qq_normal = function ( e )
{
  local (eq, ivec, n, n2, nq, neg, pos);

  n = e.nr; 
  n2 = e.nc;

  if (n2 > 1) { error("(qq_normal): e must be a vector."); }

  n = length (e);
  eq = sort(e).val;
  ivec = ((1:n)'-.5)./n;
  pos = (ivec >= .5);
  neg = !pos;
  nq = zeros (size (e));
  nq[find (neg)] = -sqrt (2)*inverf (2*(.5-ivec[find (neg)]));
  nq[find (pos)] = +sqrt (2)*inverf (2*(ivec[find (pos)]-.5));

  return [nq, eq];
}