.. index::
   single: running statistics
   single: online statistics

******************
Running Statistics
******************

This chapter describes routines for computing running statistics,
also known as online statistics, of data. These routines are
suitable for handling large datasets for which it may be
inconvenient or impractical to store in memory all at once.
The data can be processed in a single pass, one point at a time.
Each time a data point is added to the accumulator, internal
parameters are updated in order to compute the current mean, variance,
standard deviation, skewness, and kurtosis. These statistics are
exact, and are updated with numerically stable single-pass algorithms.
The median and arbitrary quantiles are also available, however these
calculations use algorithms which provide approximations, and grow
more accurate as more data is added to the accumulator.

The functions described in this chapter are declared in the header file
:file:`gsl_rstat.h`.

Initializing the Accumulator
============================

.. type:: gsl_rstat_workspace

   This workspace contains parameters used to calculate various statistics
   and are updated after each data point is added to the accumulator.

.. function:: gsl_rstat_workspace * gsl_rstat_alloc (void)

   This function allocates a workspace for computing running statistics.
   The size of the workspace is :math:`O(1)`.

.. function:: void gsl_rstat_free (gsl_rstat_workspace * w)

   This function frees the memory associated with the workspace :data:`w`.

.. function:: int gsl_rstat_reset (gsl_rstat_workspace * w)

   This function resets the workspace :data:`w` to its initial state,
   so it can begin working on a new set of data.

Adding Data to the Accumulator
==============================

.. function:: int gsl_rstat_add (const double x, gsl_rstat_workspace * w)

   This function adds the data point :data:`x` to the statistical
   accumulator, updating calculations of the mean, variance,
   standard deviation, skewness, kurtosis, and median.

.. function:: size_t gsl_rstat_n (const gsl_rstat_workspace * w)

   This function returns the number of data so far added to the accumulator.

Current Statistics
==================

.. function:: double gsl_rstat_min (const gsl_rstat_workspace * w)

   This function returns the minimum value added to the accumulator.

.. function:: double gsl_rstat_max (const gsl_rstat_workspace * w)

   This function returns the maximum value added to the accumulator.

.. function:: double gsl_rstat_mean (const gsl_rstat_workspace * w)

   This function returns the mean of all data added to the accumulator,
   defined as

   .. only:: not texinfo

      .. math:: {\Hat\mu} = {1 \over N} \sum x_i

   .. only:: texinfo

      ::

         \Hat\mu = (1/N) \sum x_i

.. function:: double gsl_rstat_variance (const gsl_rstat_workspace * w)

   This function returns the variance of all data added to the accumulator,
   defined as

   .. only:: not texinfo

      .. math:: {\Hat\sigma}^2 = {1 \over (N-1)} \sum (x_i - {\Hat\mu})^2

   .. only:: texinfo

      ::

         \Hat\sigma^2 = (1/(N-1)) \sum (x_i - \Hat\mu)^2

.. function:: double gsl_rstat_sd (const gsl_rstat_workspace * w)

   This function returns the standard deviation of all data added to the
   accumulator, defined as the square root of the variance given above.

.. function:: double gsl_rstat_sd_mean (const gsl_rstat_workspace * w)

   This function returns the standard deviation of the mean, defined as

   .. only:: not texinfo
   
      .. math:: \Hat\sigma_{\Hat\mu} = {\Hat\sigma \over \sqrt{N}}

   .. only:: texinfo

      ::

         sd_mean = \Hat\sigma / \sqrt{N}

.. function:: double gsl_rstat_rms (const gsl_rstat_workspace * w)

   This function returns the root mean square of all data added to the
   accumulator, defined as

   .. math:: rms = \sqrt{{1 \over N} \sum x_i^2}

.. function:: double gsl_rstat_skew (const gsl_rstat_workspace * w)

   This function returns the skewness of all data added to the accumulator,
   defined as

   .. only:: not texinfo

      .. math::

         skew = {1 \over N} \sum 
          {\left( x_i - {\Hat\mu} \over {\Hat\sigma} \right)}^3

   .. only:: texinfo

      ::

         skew = (1/N) \sum ((x_i - \Hat\mu)/\Hat\sigma)^3

.. function:: double gsl_rstat_kurtosis (const gsl_rstat_workspace * w)

   This function returns the kurtosis of all data added to the accumulator,
   defined as

   .. only:: not texinfo

      .. math::

         kurtosis = \left( {1 \over N} \sum 
          {\left(x_i - {\Hat\mu} \over {\Hat\sigma} \right)}^4 
          \right) 
          - 3

   .. only:: texinfo

      ::

         kurtosis = ((1/N) \sum ((x_i - \Hat\mu)/\Hat\sigma)^4)  - 3

.. function:: double gsl_rstat_median (gsl_rstat_workspace * w)

   This function returns an estimate of the median of the data added to
   the accumulator.

Quantiles
=========

The functions in this section estimate quantiles dynamically without
storing the entire dataset, using the algorithm of Jain and Chlamtec, 1985.
Only five points (markers) are stored which represent the minimum
and maximum of the data, as well as current estimates of the
:math:`p/2`-, :math:`p`-, and :math:`(1+p)/2`-quantiles. Each time
a new data point is added, the marker positions and heights are
updated.

.. type:: gsl_rstat_quantile_workspace

   This workspace contains parameters for estimating quantiles of the
   current dataset

.. function:: gsl_rstat_quantile_workspace * gsl_rstat_quantile_alloc (const double p)

   This function allocates a workspace for the dynamic estimation of
   :data:`p`-quantiles, where :data:`p` is between :math:`0` and :math:`1`.
   The median corresponds to :math:`p = 0.5`. The size of the workspace
   is :math:`O(1)`.

.. function:: void gsl_rstat_quantile_free (gsl_rstat_quantile_workspace * w)

   This function frees the memory associated with the workspace :data:`w`.

.. function:: int gsl_rstat_quantile_reset (gsl_rstat_quantile_workspace * w)

   This function resets the workspace :data:`w` to its initial state,
   so it can begin working on a new set of data.

.. function:: int gsl_rstat_quantile_add (const double x, gsl_rstat_quantile_workspace * w)

   This function updates the estimate of the :math:`p`-quantile with
   the new data point :data:`x`.

.. function:: double gsl_rstat_quantile_get (gsl_rstat_quantile_workspace * w)

   This function returns the current estimate of the :math:`p`-quantile.

Examples
========

Here is a basic example of how to use the statistical functions:

.. include:: examples/rstat.c
   :code:

The program should produce the following output,

.. include:: examples/rstat.txt
   :code:

This next program estimates the lower quartile, median and upper
quartile from 10,000 samples of a random Rayleigh distribution,
using the :math:`P^2` algorithm of Jain and Chlamtec. For
comparison, the exact values are also computed from the sorted
dataset.

.. include:: examples/rquantile.c
   :code:

The program should produce the following output,

.. include:: examples/rquantile.txt
   :code:

References and Further Reading
==============================

The algorithm used to dynamically estimate :math:`p`-quantiles is described
in the paper,

* R. Jain and I. Chlamtac.
  *The P^2 algorithm for dynamic calculation of quantiles and histograms without storing observations*,
  Communications of the ACM, Volume 28 (October), Number 10, 1985,
  p. 1076-1085.