This chapter describes the statistical functions in the library.  The basic statistical functions include routines to compute the mean, variance and standard deviation. More advanced functions allow you to calculate absolute deviations, skewness, and kurtosis as well as the median and arbitrary percentiles.  Statistical tests for comparing different datasets, such as the t-test, are also included.

All the functions are available in versions for datasets in all the standard floating-point and integer types. The versions for double precision floating-point data have the prefix gsl_stats and the versions for integer data have the prefix gsl_stats_int. The following table gives the prefixes for all the numeric types,

     gsl_stats_               double
     gsl_stats_long_double_   long double
     gsl_stats_float_         float
     gsl_stats_long_          long
     gsl_stats_ulong_         unsigned long
     gsl_stats_int_           int
     gsl_stats_uint_          unsigned int
     gsl_stats_short_         short
     gsl_stats_ushort_        unsigned short
     gsl_stats_char_          char
     gsl_stats_uchar_         unsigned char

All the algorithms use a recurrence relation to compute average quantities in a stable way, without large intermediate values that might overflow.

Mean, Standard Deviation and Variance
=====================================

 -- Statistics: double gsl_stats_mean (const double DATA[], size_t N)

     The function gsl_stats_mean returns the arithmetic mean of DATA, an array of length N.

     The sample mean is an unbiased estimator of the population mean and is denoted by $\Hat\mu$ \Hat\mu.

          {\Hat\mu} = {1 \over N} \sum x_i

     where $x_i$ x_i are the elements of the array DATA.

     For a gaussian distribution the variance of $\Hat\mu$ \Hat\mu is $\sigma^2/N$ \sigma^2/N, so the one standard-deviation error on the mean is $\sigma/\sqrtN$ \sigma/\sqrt{N}.

     If the underlying distribution has long tails then the mean may be subject to large fluctuations.

 -- Statistics: double gsl_stats_est_variance (const double DATA[], size_t N)

     The function gsl_stats_est_variance returns the estimated variance of DATA, an array of length N.

     This is an unbiased estimator of $\sigma^2$ \sigma^2, and is defined by

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

     where $x_i$ x_i are the elements of the array DATA.

     Note that the normalisation factor of 1/(N-1) is chosen to make $\Hat\sigma^2$ \Hat\sigma^2 an unbiased estimator. (A factor of 1/N produces variances which, on the average, are too small since correlations between $x_i$ x_i and $\Hat\mu$ \Hat\mu tend to reduce the sample variance).

     For a gaussian distribution the variance of $\Hat\sigma^2$ \Hat\sigma^2 is $2 \sigma^4/N$ 2 \sigma^4/N. The one standard-deviation error on the variance is $\sqrt2 \sigma^2/\sqrtN$ \sqrt{2} \sigma^2/\sqrt{N}.

     If the underlying distribution has long tails then the variance may be subject to large fluctuations, or may not converge.

     Both these functions compute the mean via a call to gsl_stats_mean. If you have already computed the mean then you can pass it directly to gsl_stats_est_variance_with_mean.

 -- Statistics: double gsl_stats_est_variance_with_mean (const double DATA[], size_t N, double MEAN)

     The function gsl_stats_est_variance_with_mean returns the sample variance of DATA relative to the given value of MEAN.  The function is computed with $\Hat\mu$ \Hat\mu replaced by the value of MEAN that you supply,

          {\Hat\sigma}^2 = {1 \over (N-1)} \sum (x_i - mean)^2

 -- Statistics: double gsl_stats_est_sd (const double DATA[], size_t N)
 -- Statistics: double gsl_stats_est_sd_with_mean (const double DATA[], size_t N, double MEAN)

     The standard deviation is defined as the square root of the variance. These functions return the square root of the corresponding variance functions above.

 -- Statistics: double gsl_stats_variance (const double DATA[], size_t N)
 -- Statistics: double gsl_stats_variance_with_mean (const double DATA[], size_t N, double MEAN)

     These functions compute an unbiased estimate of the variance when the population mean of a distribution is known a priori.

     In this case the estimator for the variance requires a factor 1/N and the sample mean $\Hat\mu$ \Hat\mu must be replaced by the known population mean $\mu$ \mu,

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

 -- Statistics: double gsl_stats_sd (const double DATA[], size_t N)
 -- Statistics: double gsl_stats_sd_with_mean (const double DATA[], size_t N, double MEAN)

     The functions gsl_stats_sd and gsl_stats_sd_with_mean calculate the standard deviation, which is simply the square root of the corresponding variance function.

Absolute deviation
==================

 -- Statistics: double gsl_stats_absdev (const double DATA[], size_t N)

     The function gsl_stats_absdev computes the absolute deviation from the mean of DATA, an array of length N.

     The absolute deviation from the mean is defined as

          absdev  = {1 \over N} \sum |x_i - {\Hat\mu}|

     where $x_i$ x_i are the elements of the array DATA.

     The absolute deviation from the mean provides a more robust measure of the width of a distribution than the variance.

     Both these functions compute the mean of DATA via a call to gsl_stats_mean or gsl_stats_int_mean.

 -- Statistics: double gsl_stats_absdev_with_mean (const double DATA[], size_t N, double MEAN)

     The function gsl_stats_absdev_with_mean computes the absolute deviation of the array DATA relative to the given value of MEAN,

          absdev  = {1 \over N} \sum |x_i - mean|

     These functions are useful if you have already computed the mean of DATA (and want to avoid recomputing it), or wish to calculate the absolute deviation relative to another value (such as zero, or the median).

Higher moments (skewness and kurtosis)
======================================

 -- Statistics: double gsl_stats_skew (const double DATA[], size_t N)
     The functions gsl_stats_skew and gsl_stats_int_skew compute the skewness of DATA, an array of length N.

     The skewness is defined as

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

     where $x_i$ x_i are the elements of the array DATA.

     The skewness measures the asymmetry of the tails of a distribution.

     Both these functions compute the mean and standard deviation of DATA via calls to gsl_stats_mean and gsl_stats_est_sd or gsl_stats_int_mean and gsl_stats_int_est_sd.

 -- Statistics: double gsl_stats_skew_with_mean_and_sd (const double DATA[], size_t N, double MEAN, double SD)

     The function gsl_stats_skew_with_mean_and_sd computes the skewness of the array DATA using the given values of the mean MEAN and standard deviation SD.

          skew = {1 \over N}
               \sum {( x_i - mean \over sd )}^3

     These functions are useful if you have already computed the mean and standard deviation of DATA and want to avoid recomputing them.

 -- Statistics: double gsl_stats_kurtosis (const double DATA[], size_t N)

     The function gsl_stats_kurtosis computes the kurtosis of DATA, an array of length N.

     The kurtosis is defined as

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

     The kurtosis measures how sharply peaked a distribution is, relative to its width. The kurtosis is normalised to zero for a gaussian distribution by the -3 in the definition.

 -- Statistics: double gsl_stats_kurtosis_with_mean_and_sd (const double DATA[], size_t N, double MEAN, double SD)

     The function gsl_stats_kurtosis_with_mean_and_sd computes the kurtosis of the array DATA using the given values of the mean MEAN and standard deviation SD.

          kurtosis = {1 \over N}
            ( \sum {(x_i - mean \over sd )}^4 ) - 3

     These functions are useful if you have already computed the mean and standard deviation of DATA and want to avoid recomputing them.

Maximum and Minimum values
==========================

 -- Statistics: double gsl_stats_max (const double DATA[], size_t N)

     The function gsl_stats_max returns the maximum value in DATA, an array of length N. The maximum value is defined as the value of the element $x_i$ x_i which satisfies $x_i \ge x_j$ x_i >= x_j for all j.

     If you want instead to find the element with the largest absolute magnitude you will need to apply fabs or abs to your data before calling this function.

 -- Statistics: double gsl_stats_min (const double DATA[], size_t N)

     The function gsl_stats_min returns the minimum value in DATA, an array of length N.  The minimum value is defined as the value of the element $x_i$ x_i which satisfies $x_i \le x_j$ x_i <= x_j for all j.

     If you want instead to find the element with the smallest absolute magnitude you will need to apply fabs or abs to your data before calling this function.

 -- Statistics: size_t gsl_stats_max_index (const double DATA[], size_t N)

     The function gsl_stats_max_index returns the index of the maximum value in DATA, an array of length N. The maximum value is defined as the value of the element $x_i$ x_i which satisfies $x_i \ge x_j$ x_i >= x_j for all j. When there are several equal maximum elements then the first one is chosen.

 -- Statistics: size_t gsl_stats_min_index (const double DATA[], size_t N)

     The function gsl_stats_min_index returns the index of the minimum value in DATA, a double-precision array of length N. The minimum value is defined as the value of the element $x_i$ x_i which satisfies $x_i \ge x_j$ x_i >= x_j for all j. When there are several equal minimum elements then the first one is chosen.

Median and Percentiles
======================

The median and percentile functions described in this section operate on sorted data. For convenience we use quantiles, measured on a scale of 0 to 1, instead of percentiles (which use a scale of 0 to 100).

 -- Statistics: void gsl_stats_sort_data (double DATA[], size_t N)

     The function gsl_stats_sort_data sorts the elements of DATA, an array of length N, into ascending numerical order in-place (i.e. after sorting, the minimum value of data is in DATA[0] and the maximum value is in DATA[N-1]). The system qsort function is used to perform the sorting.

     The function gsl_stats_int_sort_data sorts an integer array.

 -- Statistics: double gsl_stats_median_from_sorted_data (const double SORTED_DATA[], size_t N)

     The function gsl_stats_median_from_sorted_data returns the median value of SORTED_DATA, an array of length N with its elements in ascending numerical order. There are no checks to see whether the data are sorted, so the function gsl_stats_sort_data should always be used first.

     When the array has an odd number of elements the median is the value of element (n-1)/2. When the array has an even number of elements the median is the mean of the two nearest middle values, elements (n-1)/2 and n/2.

     The function gsl_stats_int_median_from_sorted_data returns the median of an integer array. Since the median may have to be found by interpolation the function gsl_stats_int_median_from_sorted_data always returns a floating-point number.

 -- Statistics: double gsl_stats_quantile_from_sorted_data (const double SORTED_DATA[], size_t N, double F)
 -- Statistics: double gsl_stats_int_quantile_from_sorted_data (const int SORTED_DATA[], size_t N, double F)

     The function gsl_stats_median_from_sorted_data returns a quantile value of SORTED_DATA, a double-precision array of length N with its elements in ascending numerical order.  The quantile is determined by the F, a fraction between 0 and 1. For example, to compute the value of the 75th percentile F should have the value 0.75.

     There are no checks to see whether the data are sorted, so the function gsl_stats_sort_data should always be used first.

     The quantile is found by interpolation, using the formula

          quantile = (1 - \delta) data[i] + \delta data[i+1]

     where i is floor($(n - 1)f$ (n - 1)f) and $\delta$ \delta is $(n-1)f - i$. (n-1)f - i.

     Thus the minimum value of the array (DATA[0]) is given by F equal to zero, the maximum value (DATA[N-1]) is given by F equal to one and the median value is given by F equal to 0.5.

     The function gsl_stats_int_quantile_from_sorted_data returns a quantile of an integer array. Since the quantile is found by interpolation the function gsl_stats_int_quantile_from_sorted_data always returns a floating-point number.

Statistical tests
=================

FIXME, do more work on the statistical tests

Example statistical programs
============================

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

     #include <stdio.h>
     #include <gsl_statistics.h>

     int main()
     {
       double data[5] = {17.2, 18.1, 16.5, 18.3, 12.6} ;
       double mean, variance, largest, smallest;

       mean     = gsl_stats_mean(data, 5);
       variance = gsl_stats_variance(data, 5);
       largest  = gsl_stats_max(data, 5);
       smallest = gsl_stats_min(data, 5);

       printf("The dataset is %g, %g, %g, %g, %g\n",
             data[0], data[1], data[2], data[3], data[4]);

       printf("The sample mean is %g\n", mean) ;
       printf("The estimated variance is %g\n", variance) ;
       printf("The largest value is %g\n", largest) ;
       printf("The smallest value is %g\n", smallest) ;
     }

The program should produce the following output,

     The dataset is 17.2, 18.1, 16.5, 18.3, 12.6
     The sample mean is 16.54
     The estimated variance is 4.2984
     The largest value is 18.3
     The smallest value is 12.6

Here is an example using sorted data,

     #include <stdio.h>
     #include <gsl_statistics.h>

     int main()
     {
       double data[5] = {17.2, 18.1, 16.5, 18.3, 12.6} ;
       double median, upperq, lowerq;

       printf("The original dataset is %g, %g, %g, %g, %g\n",
             data[0], data[1], data[2], data[3], data[4]);

       gsl_stats_sort_data(data, 5) ;

       printf("The sorted dataset is %g, %g, %g, %g, %g\n",
             data[0], data[1], data[2], data[3], data[4]);

       median   = gsl_stats_median_from_sorted_data(data, 5);
       upperq   = gsl_stats_quantile_from_sorted_data(data, 5, 0.75);
       lowerq   = gsl_stats_quantile_from_sorted_data(data, 5, 0.25);

       printf("The median is %g\n", median) ;
       printf("The upper quartile is %g\n", upperq) ;
       printf("The lower quartile is %g\n", lowerq) ;
     }

This program should produce the following output,

     The original dataset is 17.2, 18.1, 16.5, 18.3, 12.6
     The sorted dataset is 12.6, 16.5, 17.2, 18.1, 18.3
     The median is 17.2
     The upper quartile is 18.1
     The lower quartile is 16.5

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

The standard reference for almost any topic in statistics is the multi-volume Advanced Theory of Statistics by Kendall and Stuart.

  Maurice Kendall, Alan Stuart, and J. Keith Ord. The Advanced Theory of Statistics (multiple volumes) reprinted as Kendall's Advanced Theory of Statistics. Wiley, ISBN 047023380X.

Many statistical concepts can be more easily understood by a Bayesian approach. The recent book by Gelman, Carlin, Stern and Rubin gives a comprehensive coverage of the subject, while Jeffreys' Theory of Probability gives an older but more pedagogical introduction.

  Andrew Gelman, John B. Carlin, Hal S. Stern, Donald B. Rubin. Bayesian Data Analysis. Chapman & Hall, ISBN 0412039915.

  Sir Harold Jeffreys. Theory of Probability. Clarendon Press, ISBN 0198531931

For physicists the Particle Data Group provides useful reviews of Probability and Statistics in the "Mathematical Tools" section of its Annual Review of Particle Physics.

  Review of Particle Properties R.M. Barnett et al., Physical Review D54, 1 (1996)

The Review of Particle Physics is available online in postscript and pdf format at http://pdg.lbl.gov/.