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 gsl_rstat.h.

Initializing 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 O(1).

 -- Function: void gsl_rstat_free (gsl_rstat_workspace * W)
     This function frees the memory associated with the workspace W.

 -- Function: int gsl_rstat_reset (gsl_rstat_workspace * W)
     This function resets the workspace 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 X to the statistical accumulator, updating calculations of the mean, variance, standard deviation, skewness, kurtosis, and median.

 -- Function: size_t gsl_rstat_n (gsl_rstat_workspace * W)
     This function returns the number of data so far added to the accumulator.

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

 -- Function: double gsl_rstat_min (gsl_rstat_workspace * W)
     This function returns the minimum value added to the accumulator.

 -- Function: double gsl_rstat_max (gsl_rstat_workspace * W)
     This function returns the maximum value added to the accumulator.

 -- Function: double gsl_rstat_mean (gsl_rstat_workspace * W)
     This function returns the mean of all data added to the accumulator, defined as

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

 -- Function: double gsl_rstat_variance (gsl_rstat_workspace * W)
     This function returns the variance of all data added to the accumulator, defined as

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

 -- Function: double gsl_rstat_sd (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_skew (gsl_rstat_workspace * W)
     This function returns the skewness of all data added to the accumulator, defined as

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

 -- Function: double gsl_rstat_kurtosis (gsl_rstat_workspace * W)
     This function returns the kurtosis of all data added to the accumulator, defined as

          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 p/2-, p-, and (1+p)/2-quantiles. Each time a new data point is added, the marker positions and heights are updated.

 -- Function: gsl_rstat_quantile_workspace * gsl_rstat_quantile_alloc (const double P)
     This function allocates a workspace for the dynamic estimation of P-quantiles, where P is between 0 and 1. The median corresponds to p = 0.5. The size of the workspace is O(1).

 -- Function: void gsl_rstat_quantile_free (gsl_rstat_quantile_workspace * W)
     This function frees the memory associated with the workspace W.

 -- Function: int gsl_rstat_quantile_add (const double X, gsl_rstat_quantile_workspace * W)
     This function updates the estimate of the p-quantile with the new data point X.

 -- Function: double gsl_rstat_quantile_get (gsl_rstat_quantile_workspace * W)
     This function returns the current estimate of the p-quantile.

Examples
========

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

          [file rstat.c not included]

The program should produce the following output,

          [file rstat.txt not included]

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

          [file rquantile.c not included]

The program should produce the following output,

          [file rquantile.txt not included]

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

The algorithm used to dynamically estimate 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.