Index: /trunk/doc/pslib/psLibADD.tex
===================================================================
--- /trunk/doc/pslib/psLibADD.tex	(revision 302)
+++ /trunk/doc/pslib/psLibADD.tex	(revision 303)
@@ -12,6 +12,9 @@
 \docnumber{PSDC-430-006}
 
+\newcommand\citealt{}
+
 \begin{document}
 \maketitle
+
 
 % -- Revision History --
@@ -252,3 +255,182 @@
 \end{center}
 
+\subsection{Sorting}
+
+A variety of sorting algorithms exist, with a wide range in speed for
+different set sizes.  Three of the best sorting algorithms are
+``heapsort'', ``quicksort'', and ``mergesort'' (see, e.g.,
+\citealt{knuth}, \citealt{sedgewick}, \citealt{press}).  These three
+sorting algorithms all run in a time of $O(n \log n)$ in the best
+case, but have different worse-case run times.  Implementations of all
+three sorting algorithms are described in references above and in
+various places online (e.g.,
+http://www.iti.fh-flensburg.de/lang/algorithmen/sortieren/algoen.htm).
+The linux \code{qsort} function orignially mades use of the quicksort
+algorithm, but now uses \tbd{what algorithm?}.  The GSL function
+\code{gsl_heapsort} uses the heapsort algorithm.  Both of these
+implemenations require a pointer to the comparison function, which may
+result in a slower code.  
+
+For PSLib, the sorting functions shall be implemented via
+\code{qsort}.  The function \code{psSort} shall return the sorted
+result of the input array without over-writing the input array.  The
+function \code{psSortIndex} shall return an integer index to the
+sequence of the input array without overwriting the input array.
+Given the following line of code, \code{out = psSortIndex (NULL,
+&in);}, the elements of the array \code{out} are in the sequence
+\code{in.arr[out->arr[0]]} to \code{in.arr[out->arr[in.n - 1]]}. 
+
+\subsection{Smoothing: Boxcar and Gaussian}
+\label{smooth}
+
+Smoothing may occasionally be perfomed on data.  We present the
+algorithms for two typical versions: boxcar and Gaussian smoothing.
+In both smoothing techniques, given a series of data values $f_i$, the
+smoothed values $g_i$ are determined by calculating a linear
+combination based on the input data point and its nearest $2N$
+neighbors in the form:
+
+\[ g_i = \sum_{n=-N}^N c_n f_i \]
+%
+where the values of $c_n$ determine the filter type.  For boxcar
+smoothing, the values $c_n$ are constant, and must be equal to $1/(2N
++ 1)$ to maintain the zeroth moment of the data (care must be taken at
+the ends of the data range to reduce the value of $c_n$ as fewer input
+data points may be used).  For Gaussian smoothing, the crucial
+parameter is $\sigma$, the standard deviation.  The
+value of $N$ should be chosen to be large enough, $N = 5\sigma$, and
+the values of $c_n$ are then just the Gaussian curve:
+
+\[ c_n = \frac{e^{\frac{-n^2}{2\sigma^2}}}{\sqrt{2\pi\sigma^2}} \]
+
+\subsection{Statistics}
+
+The general statistics function \code{psStats} performs a variety of
+statistical calculations on a collection of floating point numbers.
+These statistics include both sample values, in which the formal
+statistic is calculated for the complete sample, and robust values, in
+which the statistic is estimated on the basis of an assumption of an
+underlying population.  While more reliable in general, the robust
+estimators may be somewhat slower than the sample statisic, so their
+use may vary depending on the context.  In addition, certain
+sigma-clipped values are defined as an intermediate choice between the
+sample and robust estimators.  
+
+\subsubsection{Sample Statistics}
+
+We define the following statistical terms, assuming there is a set of
+data elements $x_i$.
+
+\paragraph{Mean}
+
+The mean is defined as: \[ \bar{x} = \frac{1}{N} \sum_{i = 1}^{N} x_i \]
+
+\paragraph{Median}
+
+The median is defined as the value for which 50\% of the data values
+are larger and 50\% are smaller.  For a sample, the values are sorted
+and the median is given by the value of the middle element, if the
+number of values is odd, and by the average of the two middle values
+if the number of values is even.  This median should be avoided for
+samples which are large (e.g., $N > 10000$ elements) as the basic
+robust median is quicker and more accurate.  
+
+\paragraph{Upper and Lower Quartiles}
+
+The upper and lower quartiles ($U_{\frac{1}{4}}$ and
+$L_{\frac{1}{4}}$) are first and last quarter equivalents to the
+median.  The upper quartile is the value for which 25\% of the data
+are larger and 75\% are smaller.  The lower quartile is the value for
+which 75\% of the data are larger and 25\% are smaller.  For a sample,
+the values are sorted and the quartiles are given by the values of
+elements $N/4$ and $3N/4$.  For the purpose of a sample upper and
+lower quartile, it is sufficient to provide the closest integer entry
+to these values.  The sample quartiles should be avoided for samples
+which are large (e.g., $N > 10000$ elements) as the robust quartiles
+are quicker and more accurate.
+
+\paragraph{Standard Deviation}
+
+The standard deviation of the sample is given by:
+
+\[ \sigma = \sqrt{\sum_{i = 1}^N \frac{(x_i - \bar{x})^2}{N - 1}} \]
+
+To minimize the numerical rounding error, this should be calculated
+numerically as:
+
+\[ \sigma = \sqrt{\frac{1}{N - 1} [ \sum_{i = 1}^{N} (x_i - \bar{x})^2 - \frac{1}{N} (\sum_{i = 1}^{N} (x_i - \bar{x}))^2 ]} \]
+
+\subsubsection{Clipped Statistics}
+
+The clipped statistics are used to determine the mean and $\sigma$ of
+a distribution in the presence of outliers. The clipped statistics are
+quicker than the complete robust statistical estimators and more
+reliable than the sample statistics.  The clipped statistics algorithm
+produces the clipped mean and clipped standard deviation.  The
+algorithm has 2 parameters: $N$, the number of iterations and $k$, the
+multiplying factor of $\sigma$ used to exclude outliers.  Typical
+values for both $N$ and $k$ are 3.  The algorithm is as follows:
+
+\begin{enumerate}
+\item Compute the sample median. 
+\item Compute the sample standard deviation.
+\item Use the sample median as the first estimator of the mean, $\bar{x}$.
+\item Use the sample standard deviation as the first estimator of the
+  true standard deviation, $\sigma$.
+\item Repeat the following N times:
+  \begin{enumerate}
+  \item Exclude all values $x_{i}$ for while $|x_{i} - \bar{x}| > k \sigma$.
+  \item Compute the mean and standard deviation of the sub-sample.
+  \item Use the new mean for $\bar{x}$.
+  \item Use the new standard deviations for $\sigma$.
+\end{enumerate}
+\item The last calulcated value of $\bar{x}$ is the clipped mean.
+\item The last calulcated value of $\sigma$ is the clipped standard deviation.
+\end{enumerate}
+
+\subsubsection{Robust Statistics}
+
+The robust version of the statistics provide estimators of basic
+statistical concepts which are reliable even for data samples with
+significant contamination.  A typical case is the situation in which
+the data of interest represent a primary population of interest with a
+single-valued mean and standard deviation and a secondary population
+of data with a substantially different distribution.  For example, an
+image of an uncrowded night-time field may consist of a sparse
+collection of stars and an overall background level.  The majority of
+pixels have the background value, with some variance due to noise
+sources, but many are significantly higher (contributed by stars) or
+significantly lower (dead pixels).  If we want to measure the mean of
+the background, a robust mean is necessary as the outliers will
+strongly bias the sample statistics.
+
+The robust statistics are calculated by constructing a histogram of
+the values and performing the measurements on the histogram.  The
+choice of a bin size requires some care.  If the data are integer
+valued, the natural bin size is an integer.  Otherwise, the bin should
+be a fraction of an estimate of the standard deviation.  Use the
+sample upper and lower quartiles to determine an estimate of the
+standard deviation: $\sigma_e = (U_{\frac{1}{4}} - L_{\frac{1}{4}}) /
+1.34$.  The bin size shall be set at $\sigma_e / 10$.  The remaining
+steps of the algorithm are as follows:
+
+\begin{itemize}
+\item Construct the histogram with the specified bin size.
+\item Smooth the histogram by a Gaussian with $\sigma_s = \sigma_e /
+  4$.  
+\item Find the bin with the peak value in the range $L_{\frac{1}{4}}$
+  to $U_{\frac{1}{4}}$; this is the robust mode $\mbox{mode}_r$.  
+\item Determine $dL = (U_{\frac{1}{4}} - L_{\frac{1}{4}}) / 8$.
+\item Fit a Gaussian to the bins in the range $\mbox{mode}_r - dL$ to
+  $\mbox{mode}_r + dL$.
+\item The resulting fit parameters are the robust mean,
+$\mbox{mean}_r$ and the robust standard deviation, $\sigma_r$.
+\end{itemize}
+
+To determine the robust median, construct the cumulative histogram
+from the histogram above. Select the bin which contains the 50th
+percentile value and its two neighbors.  Fit a quadratic to these
+three points.  The robust median value is the coordinate of the
+quadratic which returns the 50\% value.
+
 \end{document}