This chapter describes functions for creating histograms.  Histograms provide a convenient way of summarizing the distribution of a set of data. A histogram consists of a set of bins which count the number of events falling into a given range of a continuous variable x. In GSL the bins of a histogram contain floating-point numbers, so they can be used to record both integer and non-integer distributions. The bins can use arbitrary sets of ranges (uniformly spaced bins are the default). Both one and two-dimensional histograms are supported.

Once a histogram has been created it can also be converted into a probability distribution function. The library provides efficient routines for selecting random samples from probability distributions. This can be useful for generating simulations based real data.

The histogram struct
====================

A histogram is defined by the following struct,

 -- Data Type: gsl_histogram

     size_t n
          This is the number of histogram bins

     double * range
          The ranges of the bins are stored in an array of N+1 elements pointed to by RANGE.

     double * bin
          The counts for each bin are stored in an array of N elements pointed to by BIN. The bins are floating-point numbers, so you can increment them by non-integer values if necessary.

The range for BIN[i] is given by RANGE[i] to RANGE[i+1]. For n bins there are n+1 entries in the array RANGE. Each bin is inclusive at the lower end and exclusive at the upper end. Mathematically this means that the bins are defined by the following inequality,

     bin[i] corresponds to range[i] <= x < range[i+1]

Here is a diagram of the correspondence between ranges and bins on the number-line for x,

        r[0]      r[1]      r[2]      r[3]      r[4]      r[5]
       ---|---------|---------|---------|---------|---------|---  x
          [ bin[0] )[ bin[1] )[ bin[2] )[ bin[3] )[ bin[5] )

In this picture the values of the RANGE array are denoted by r. On the left-hand side of each bin the square bracket "[" denotes an inclusive lower bound ($r \le x$ r <= x), and the round parentheses ")" on the right-hand side denote an exclusive upper bound (x < r). Thus any samples which fall on the upper end of the histogram are excluded. If you want to include this value for the last bin you will need to add an extra bin to your histogram.

The gsl_histogram struct and its associated functions are defined in the header file gsl_histogram.h.

Histogram allocation
====================

The functions for allocating memory to a histogram follow the style of malloc and free. In addition they also perform their own error checking. If there is insufficient memory available to allocate a histogram then the functions call the GSL error handler (with an error number of GSL_ENOMEM) in addition to returning a null pointer. Thus if you use the library error handler to abort your program then it isn't necessary to check every histogram alloc.

 -- Function: gsl_histogram * gsl_histogram_calloc (size_t N)
     This function allocates memory for a histogram with N bins, and returns a pointer to its newly initialized gsl_histogram struct. The bins are uniformly spaced with a total range of $0 \le  x < n$ 0 <=  x < n, as shown in the table below.

          bin[0] corresponds to 0 \le  x < 1
          bin[1] corresponds to 1 \le  x < 2
          ......
          bin[n-1] corresponds to n-1 \le  x < n

     The bins are initialized to zero so the histogram is ready for use.

     If insufficient memory is available a null pointer is returned and the error handler is invoked with an error code of GSL_ENOMEM.

 -- Function: gsl_histogram * gsl_histogram_calloc_uniform (size_t N, double XMIN, double XMAX)
     This function allocates memory for a histogram with N uniformly spaced bins from XMIN to XMAX, and returns a pointer to the newly initialized gsl_histogram struct. The bins are shown in the table below,

          bin[0] corresponds to xmin \le  x < xmin + d
          bin[1] corresponds to xmin + d \le  x < xmin + 2 d
          ......
          bin[n-1] corresponds to xmin + (n-1)d \le  x < xmax

     where d is the bin spacing, (xmax-xmin)/n. Each bin is initialized to zero.

     If insufficient memory is available a null pointer is returned and the error handler is invoked with an error code of GSL_ENOMEM.

To create a histogram with non-uniform bins you will need to call gsl_histogram_calloc to prepare a new histogram struct and then modify the range array to use your desired bin limits. The ranges can be arbitrary, subject to the restriction that they are monotonically increasing.

For example, the following code fragment shows how to create a histogram with logarithmic bins from 1---10, 10---100 and 100---1000.

     gsl_histogram * h = gsl_histogram_calloc (3) ;

     h->range[0] = 1.0 ;     /* bin[0] covers the range 1 <= x < 10 */
     h->range[1] = 10.0 ;    /* bin[1] covers the range 10 <= x < 100 */
     h->range[2] = 100.0 ;   /* bin[2] covers the range 100 <= x < 1000 */
     h->range[3] = 1000.0 ;

Note that the size of the RANGE array is automatically defined as double range[4] by gsl_histogram_calloc, and is one element bigger than the array of bins double bin[3]. Thus the range array safely includes extra space for the final upper value, RANGE[3].

 -- Function: void gsl_histogram_free (gsl_histogram * h)
     This function frees the histogram H and all of the memory associated with it.

Updating and accessing histogram elements
=========================================

There are two ways to access histogram bins, either by specifying an x coordinate or by using the bin-index directly. The functions for accessing the histogram through x coordinates use a binary search to identify the bin which covers the appropriate range.

 -- Function: int gsl_histogram_increment (gsl_histogram * h, double x)
     This function updates the histogram H by adding one (1.0) to the bin whose range contains the coordinate X.

     If X lies in the valid range of the histogram then the function returns zero to indicate success. If X is less than the lower limit of the histogram then the function returns GSL_EDOM, and none of bins are modified. Similarly, if the value of X is greater than or equal to the upper limit of the histogram then the function returns GSL_EDOM, and none of the bins are modified. The error handler is not called, however, since it is often necessary to compute histogram for a small range of a larger dataset, ignoring the values outside the range of interest.

 -- Function: int gsl_histogram_accumulate (gsl_histogram * H, double X, double WEIGHT)
     This function is similar to gsl_histogram_increment but increases the value of the appropriate bin in the histogram H by the floating-point number WEIGHT.

 -- Function: double gsl_histogram_get (const gsl_histogram * H, size_t I)
     This function returns the contents of the Ith bin of the histogram H. If I lies outside the valid range of indices for the histogram then the error handler is called with an error code of GSL_EDOM and the function returns 0.

 -- Function: int gsl_histogram_get_range (const gsl_histogram * H, size_t I, double * LOWER, double * UPPER)
     This function finds the upper and lower range limits of the Ith bin of the histogram H. If the index I is valid then the corresponding range limits are stored in LOWER and UPPER. The lower limit is inclusive (i.e. events with this coordinate are included in the bin) and the upper limit is exclusive (i.e. events with the coordinate of the upper limit are excluded and fall in the neighboring higher bin, if it exists). The function returns 0 to indicate success. If I lies outside the valid range of indices for the histogram then the error handler is called and the function returns an error code of GSL_EDOM.

 -- Function: double gsl_histogram_max (const gsl_histogram * H)
 -- Function: double gsl_histogram_min (const gsl_histogram * H)
 -- Function: size_t gsl_histogram_bins (const gsl_histogram * H)
     These functions return the maximum upper and minimum lower range limits and the number of bins of the histogram H. They provide a way of determining these values without accessing the gsl_histogram struct directly.

 -- Function: void gsl_histogram_reset (gsl_histogram * H)
     This function resets all the bins in the histogram H to zero.

Searching histogram ranges
==========================

The following functions are used by the access and update routines to locate the bin which corresponds to a given x coordinate.

 -- Function: int gsl_histogram_find_impl (size_t N, const double RANGE[], double X, size_t * I)
     This function finds and sets the index I to the offset in the array RANGE of size N which bounds the value of X, such that range[i] \le  x < range[i+1]. The binary search function bsearch from the system C-library is used to locate the appropriate range. If a suitable value of I is found then the function returns 0 to indicate success. If X is less than the lower limit the function returns -1, and if X is greater than or equal to the upper limit it returns +1. The error handler is not called.

 -- Function: int gsl_histogram_find (const gsl_histogram * H, double X, size_t * I)
     This function uses gsl_histogram_find_impl to set the index I to the bin number which covers the coordinate X in the histogram H. If X is found then the function sets the index I and returns zero to indicate success. If X lies outside the valid range of the histogram then the function returns GSL_EDOM and the error handler is invoked.

Reading and writing histograms
==============================

The library provides functions for reading and writing histograms to a file as binary data or formatted text.

 -- Function: int gsl_histogram_fwrite (FILE * STREAM, const gsl_histogram * H)
     This function writes the ranges and bins of the histogram H to the stream STREAM in binary format. The return value is 0 for success and GSL_EFAILED if there was a problem writing to the file. Since the data is written in the native binary format it may not be portable between different architectures.

 -- Function: int gsl_histogram_fread (FILE * STREAM, gsl_histogram * H)
     This function reads into the histogram H from the open stream STREAM in binary format. The histogram H must be preallocated with the correct size since the function uses the number of bins in H to determine how many bytes to read. The return value is 0 for success and GSL_EFAILED if there was a problem reading from the file. The data is assumed to have been written in the native binary format on the same architecture.

 -- Function: int gsl_histogram_fprintf (FILE * STREAM, const gsl_histogram * H, const char * RANGE_FORMAT, const char * BIN_FORMAT)
     This function writes the ranges and bins of the histogram H line-by-line to the stream STREAM using the format specifiers RANGE_FORMAT and BIN_FORMAT. These should be one of the %g, %e or %f formats for floating point numbers. The function returns 0 for success and GSL_EFAILED if there was a problem writing to the file. The histogram output is formatted in three columns, and the columns are separated by spaces, like this,

          range[0] range[1] bin[0]
          range[1] range[2] bin[1]
          range[2] range[3] bin[2]
          ....
          range[n-1] range[n] bin[n-1]

     The values of the ranges are formatted using RANGE_FORMAT and the value of the bins are formatted using BIN_FORMAT. Each line contains the lower and upper limit of the range of the bins and the value of the bin itself. Since the upper limit of one bin is the lower limit of the next there is duplication of these values between lines but this allows the histogram to be manipulated with line-oriented tools.

 -- Function: int gsl_histogram_fscanf (FILE * STREAM, gsl_histogram * H)
     This function reads formatted data from the stream STREAM into the histogram H. The data is assumed to be in the three-column format used by gsl_histogram_fprintf. The histogram H must be preallocated with the correct length since the function uses the size of H to determine how many numbers to read. The function returns 0 for success and GSL_EFAILED if there was a problem reading from the file.

Resampling from histograms
==========================

A histogram made by counting events can be regarded as a measurement of a probability distribution. Allowing for statistical error, the height of each bin represents the probability of an event where the value of x falls in the range of that bin. The probability distribution function has the one-dimensional form p(x)dx where,

     p(x) = n_i/ (N w_i)

In this equation n_i is the number of events in the bin which contains x, w_i is the width of the bin and N is the total number of events. The distribution of events within each bin is assumed to be uniform.

The histogram probability distribution struct
=============================================

The probability distribution function for a histogram consists of a set of bins which measure the probability of an event falling into a given range of a continuous variable x. A probability distribution function is defined by the following struct, which actually stores the cumulative probability distribution function. This is the natural quantity for generating samples via the inverse transform method, because there is a one-to-one mapping between the cumulative probability distribution and the range [0,1]. It can be shown that by taking a uniform random number in this range and finding its corresponding coordinate in the cumulative probability distribution we obtain samples with the desired probability distribution.

 -- Data Type: gsl_histogram_pdf

     size_t n
          This is the number of bins used to approximate the probability distribution function.

     double * range
          The ranges of the bins are stored in an array of N+1 elements pointed to by RANGE.

     double * sum
          The cumulative probability for the bins is stored in an array of N elements pointed to by SUM.

The following functions allow you to create a gsl_histogram_pdf struct which represents this probability distribution and generate random samples from it.

 -- Function: gsl_histogram_pdf * gsl_histogram_pdf_alloc (const gsl_histogram * H)
     This function allocates memory for a probability distribution calculated from the histogram H and returns a pointer to a newly initialized gsl_histogram_pdf struct. If any of the bins of H are negative then a null pointer is returned and the error handler is invoked with an error code of GSL_EDOM because a probability distribution cannot contain negative values.

     If insufficient memory is available a null pointer is returned and the error handler is invoked with an error code of GSL_ENOMEM.

 -- Function: void gsl_histogram_pdf_free (gsl_histogram_pdf * P)
     This function frees the probability distribution function P and all of the memory associated with it.

 -- Function: double gsl_histogram_pdf_sample (const gsl_histogram_pdf * P, double R)
     This function uses R, a uniform random number between zero and one, to compute a single random sample from the probability distribution P. The algorithm used to compute the sample s is given by the following formula,

          s = range[i] + delta * (range[i+1] - range[i])

     where i is the index which satisfies $sum[i] \le  r < sum[i+1]$ sum[i] <=  r < sum[i+1] and $\delta$ delta is $(r - sum[i])/(sum[i+1] - sum[i])$ (r - sum[i])/(sum[i+1] - sum[i]).

Example programs for histograms
===============================

The following program shows how to make a simple histogram of a column of numerical data supplied on stdin. The program takes three arguments, specifying the upper and lower bounds of the histogram and the number of bins. It then reads numbers from stdin, one line at a time, and adds them to the histogram. When there is no more data to read it prints out the accumulated histogram using gsl_histogram_fprintf.

     #include <stdio.h>
     #include <stdlib.h>
     #include <gsl_histogram.h>

     int
     main (int argc, char **argv)
     {
       double a, b ;
       size_t n;

       if (argc != 4)
         {
           printf ("Usage: gsl-histogram xmin xmax n\n"
                   "Computes a histogram of the data on stdin"
                   "using n bins from xmin to xmax\n");
           exit (0);
         }

       a = atof (argv[1]);
       b = atof (argv[2]);
       n = atoi (argv[3]);

       {
         gsl_histogram * h = gsl_histogram_calloc_uniform (n, a, b) ;
         int status ;

         do {
           double x ;
           status = fscanf(stdin, "%lg", &x) ;

           gsl_histogram_increment (h, x) ;

         } while (status == 1) ;

         gsl_histogram_fprintf (stdout, h, "%g", "%g") ;
       }

       exit (0) ;
     }

Here is an example of the program in use. We generate 10000 random samples from a lorentz distribution with a width of 30 and histogram them over the range -100 to 100, using 200 bins.

     gsl-lorentz 30 10000 | gsl-histogram -100 100 200 > histogram.dat

A plot of the resulting histogram shows the familiar shape of the lorentz distribution and the fluctuations caused by the finite sample size.

     gnuplot> plot 'histogram.dat' using 1:3 with step

Two dimensional histograms
==========================

A two dimensional histogram consists of a set of bins which count the number of events falling in a given area of the (x,y) plane. The simplest way to use a two dimensional histogram is to record two-dimensional position information, n(x,y). Another possibility is to form a joint distribution by recording related variables. For example a detector might record both the position of an event (x) and the amount of energy it deposited E. These could be histogrammed as the joint distribution n(x,E).

The 2D histogram struct
=======================

Two dimensional histograms are defined by the following struct,

 -- Data Type: gsl_histogram2d

     size_t nx, ny
          This is the number of histogram bins in the x and y directions.

     double * xrange
          The ranges of the bins in the x-direction are stored in an array of NX + 1 elements pointed to by XRANGE.

     double * yrange
          The ranges of the bins in the y-direction are stored in an array of NY + 1 pointed to by YRANGE.

     double * bin
          The counts for each bin are stored in an array pointed to by BIN. The bins are floating-point numbers, so you can increment them by non-integer values if necessary. The array BIN stores the two dimensional array of bins in a single block of memory according to the mapping BIN(I,J) = bin[i * ny + j].

The range for BIN(I,J) is given by XRANGE[I] to XRANGE[I+1] in the x-direction and YRANGE[J] to YRANGE[J+1] in the y-direction. Each bin is inclusive at the lower end and exclusive at the upper end. Mathematically this means that the bins are defined by the following inequality,

     bin(i,j) corresponds to xrange[i] <= x < xrange[i+1]
                         and yrange[j] <= y < yrange[j+1]

Note that any samples which fall on the upper sides of the histogram are excluded. If you want to include these values for the side bins you will need to add an extra row or column to your histogram.

The gsl_histogram2d struct and its associated functions are defined in the header file gsl_histogram2d.h.

2D Histogram allocation
=======================

The functions for allocating memory to a 2D histogram follow the style of malloc and free. In addition they also perform their own error checking. If there is insufficient memory available to allocate a histogram then the functions call the GSL error handler (with an error number of GSL_ENOMEM) in addition to returning a null pointer. Thus if you use the library error handler to abort your program then it isn't necessary to check every 2D histogram alloc.

 -- Function: gsl_histogram2d * gsl_histogram2d_calloc (size_t NX, size_t NY)
     This function allocates memory for a two-dimensional histogram with NX bins in the x direction and NY bins in the y direction. The function returns a pointer to a newly initialized gsl_histogram2d struct. The bins are uniformly spaced with a total range of $0 \le  x < nx$ 0 <= x < nx in the x-direction and $0 \le  y < ny$ 0 <=  y < ny in the y-direction, as shown in the table below.

     The bins are initialized to zero so the histogram is ready for use.

     If insufficient memory is available a null pointer is returned and the error handler is invoked with an error code of GSL_ENOMEM.

 -- Function: gsl_histogram2d * gsl_histogram2d_calloc_uniform (size_t NX, size_t NY, double XMIN, double XMAX, double YMIN, double YMAX)

 -- Function: void gsl_histogram2d_free (gsl_histogram2d * H)
     This function frees the 2D histogram H and all of the memory associated with it.

Updating and accessing 2D histogram elements
============================================

You can access the bins of a two-dimensional histogram either by specifying a pair of (x,y) coordinates or by using the bin indices (i,j) directly. The functions for accessing the histogram through (x,y) coordinates use binary searches in the x and y directions to identify the bin which covers the appropriate range.

 -- Function: int gsl_histogram2d_increment (gsl_histogram2d * H, double X, double Y)
     This function updates the histogram H by adding one (1.0) to the bin whose x and y ranges contain the coordinates (X,Y).

     If the point (x,y) lies inside the valid ranges of the histogram then the function returns zero to indicate success. If (x,y) lies outside the limits of the histogram then the function returns GSL_EDOM, and none of bins are modified. The error handler is not called, since it is often necessary to compute histogram for a small range of a larger dataset, ignoring any coordinates outside the range of interest.

 -- Function: int gsl_histogram2d_accumulate (gsl_histogram2d * H, double X, double Y, double WEIGHT)
     This function is similar to gsl_histogram2d_increment but increases the value of the appropriate bin in the histogram H by the floating-point number WEIGHT.

 -- Function: double gsl_histogram2d_get (const gsl_histogram2d * H, size_t I, size_t J)
     This function returns the contents of the (I,J)th bin of the histogram H. If (I,J) lies outside the valid range of indices for the histogram then the error handler is called with an error code of GSL_EDOM and the function returns 0.

 -- Function: int gsl_histogram2d_get_xrange (const gsl_histogram2d * H, size_t I, double * XLOWER, double * XUPPER)
 -- Function: int gsl_histogram2d_get_yrange (const gsl_histogram2d * H, size_t J, double * YLOWER, double * YUPPER)
     These functions find the upper and lower range limits of the Ith and Jth bins in the x and y directions of the histogram H. The range limits are stored in XLOWER and XUPPER or YLOWER and YUPPER. The lower limits are inclusive (i.e. events with these coordinates are included in the bin) and the upper limits are exclusive (i.e. events with the value of the upper limit are not included and fall in the neighboring higher bin, if it exists). The functions return 0 to indicate success. If I or J lies outside the valid range of indices for the histogram then the error handler is called with an error code of GSL_EDOM.

 -- Function: double gsl_histogram2d_xmax (const gsl_histogram2d * H)
 -- Function: double gsl_histogram2d_xmin (const gsl_histogram2d * H)
 -- Function: size_t gsl_histogram2d_nx (const gsl_histogram2d * H)
 -- Function: double gsl_histogram2d_ymax (const gsl_histogram2d * H)
 -- Function: double gsl_histogram2d_ymin (const gsl_histogram2d * H)
 -- Function: size_t gsl_histogram2d_ny (const gsl_histogram2d * H)
     These functions return the maximum upper and minimum lower range limits and the number of bins for the x and y directions of the histogram H. They provide a way of determining these values without accessing the gsl_histogram2d struct directly.

 -- Function: void gsl_histogram2d_reset (gsl_histogram2d * H)
     This function resets all the bins of the histogram H to zero.

Searching 2D histogram ranges
=============================

The following functions are used by the access and update routines to locate the bin which corresponds to a given (x\,y) coordinate.

 -- Function: int gsl_histogram2d_find_impl (const gsl_histogram2d * H, double X, double Y, size_t * I, size_t * J)
     This function finds and sets the indices I and J to the offsets in the arrays XRANGE and YRANGE which bound the value of (x,y), such that $xrange[i] \le x < xrange[i+1]$ xrange[i] <= x < xrange[i+1] and $yrange[j] \le y < yrange[j+1]$ yrange[j] <= y < yrange[j+1]. The binary search function bsearch from the system C-library is used to locate the appropriate ranges. If suitable values of I and J are found then the function returns 0 to indicate success. If X is less than the lower limit the function returns -1, and if X is greater than or equal to the upper limit it returns +1. If X is valid then the same checks are applied to Y, with return values of -1 or +1 indicating and error. The error handler is not called.

 -- Function: int gsl_histogram2d_find (const gsl_histogram2d * H, double X, double Y, size_t * I, size_t * J)
     This function uses gsl_histogram2d_find_impl to set the indices I and J to the bin which covers the coordinates (x,y) in the histogram H. If (x,y) is found then the function sets the indices (I,J) and returns zero to indicate success. If (x,y) lies outside the valid range of the histogram then the function returns GSL_EDOM and the error handler is invoked.

Reading and writing 2D histograms
=================================

The library provides functions for reading and writing two dimensional histograms to a file as binary data or formatted text.

 -- Function: int gsl_histogram2d_fwrite (FILE * STREAM, const gsl_histogram2d * H)
     This function writes the ranges and bins of the histogram H to the stream STREAM in binary format. The return value is 0 for success and GSL_EFAILED if there was a problem writing to the file. Since the data is written in the native binary format it may not be portable between different architectures.

 -- Function: int gsl_histogram2d_fread (FILE * STREAM, gsl_histogram2d * H)
     This function reads into the histogram H from the stream STREAM in binary format. The histogram H must be preallocated with the correct size since the function uses the number of x and y bins in H to determine how many bytes to read. The return value is 0 for success and GSL_EFAILED if there was a problem reading from the file. The data is assumed to have been written in the native binary format on the same architecture.

 -- Function: int gsl_histogram2d_fprintf (FILE * STREAM, const gsl_histogram2d * H, const char * RANGE_FORMAT, const char * BIN_FORMAT)
     This function writes the ranges and bins of the histogram H line-by-line to the stream STREAM using the format specifiers RANGE_FORMAT and BIN_FORMAT. These should be one of the %g, %e or %f formats for floating point numbers. The function returns 0 for success and GSL_EFAILED if there was a problem writing to the file. The histogram output is formatted in five columns, and the columns are separated by spaces, like this,

          xrange[0] xrange[1] yrange[0] yrange[1] bin(0,0)
          xrange[0] xrange[1] yrange[1] yrange[2] bin(0,1)
          xrange[0] xrange[1] yrange[2] yrange[3] bin(0,2)
          ....
          xrange[0] xrange[1] yrange[ny-1] yrange[ny] bin(0,ny-1)

          xrange[1] xrange[2] yrange[0] yrange[1] bin(1,0)
          xrange[1] xrange[2] yrange[1] yrange[2] bin(1,1)
          xrange[1] xrange[2] yrange[1] yrange[2] bin(1,2)
          ....
          xrange[1] xrange[2] yrange[ny-1] yrange[ny] bin(1,ny-1)

          ....

          xrange[nx-1] xrange[nx] yrange[0] yrange[1] bin(nx-1,0)
          xrange[nx-1] xrange[nx] yrange[1] yrange[2] bin(nx-1,1)
          xrange[nx-1] xrange[nx] yrange[1] yrange[2] bin(nx-1,2)
          ....
          xrange[nx-1] xrange[nx] yrange[ny-1] yrange[ny] bin(nx-1,ny-1)

     Each line contains the lower and upper limits of the bin and the contents of the bin. Since the upper limits of the each bin are the lower limits of the neighbouring bins there is duplication of these values but this allows the histogram to be manipulated with line-oriented tools.

 -- Function: int gsl_histogram2d_fscanf (FILE * STREAM, gsl_histogram2d * H)
     This function reads formatted data from the stream STREAM into the histogram H. The data is assumed to be in the five-column format used by gsl_histogram_fprintf. The histogram H must be preallocated with the correct lengths since the function uses the sizes of H to determine how many numbers to read. The function returns 0 for success and GSL_EFAILED if there was a problem reading from the file.

Resampling from 2D histograms
=============================

As in the one-dimensional case, a two-dimensional histogram made by counting events can be regarded as a measurement of a probability distribution. Allowing for statistical error, the height of each bin represents the probability of an event where (x,y) falls in the range of that bin. For a two-dimensional histogram the probability distribution takes the form p(x,y) dx dy where,

     p(x,y) = n_{ij}/ (N A_{ij})

In this equation $n_ij$ n_{ij} is the number of events in the bin which contains (x,y), $A_ij$ A_{ij} is the area of the bin and N is the total number of events. The distribution of events within each bin is assumed to be uniform.

 -- Data Type: gsl_histogram2d_pdf

     size_t nx, ny
          This is the number of histogram bins used to approximate the probability distribution function in the x and y directions.

     double * xrange
          The ranges of the bins in the x-direction are stored in an array of NX + 1 elements pointed to by XRANGE.

     double * yrange
          The ranges of the bins in the y-direction are stored in an array of NY + 1 pointed to by YRANGE.

     double * sum
          The cumulative probability for the bins is stored in an array of NX*NY elements pointed to by SUM.

The following functions allow you to create a gsl_histogram2d_pdf struct which represents a two dimensional probability distribution and generate random samples from it.

 -- Function: gsl_histogram2d_pdf * gsl_histogram2d_pdf_alloc (const gsl_histogram2d * H)
     This function allocates memory for a two-dimensional probability distribution calculated from the histogram H and returns a pointer to a newly initialized gsl_histogram2d_pdf struct. If any of the bins of H are negative then a null pointer is returned and the error handler is invoked with an error code of GSL_EDOM because a probability distribution cannot contain negative values.

     If insufficient memory is available a null pointer is returned and the error handler is invoked with an error code of GSL_ENOMEM.

 -- Function: void gsl_histogram2d_pdf_free (gsl_histogram2d_pdf * P)
     This function frees the two-dimensional probability distribution function P and all of the memory associated with it.

 -- Function: int gsl_histogram2d_pdf_sample (const gsl_histogram2d_pdf * P, double R1, double R2, double * X, double * Y)
     This function uses two uniform random numbers between zero and one, R1 and R2, to compute a single random sample from the two-dimensional probability distribution P.

Example programs for 2D histograms
==================================

This program demonstrates two features of two-dimensional histograms. First a 10 by 10 2d-histogram is created with x and y running from 0 to 1. Then a few sample points are added to the histogram, at (0.3,0.3) with a height of 1, at (0.8,0.1) with a height of 5 and at (0.7,0.9) with a height of 0.5. This histogram with three events is used to generate a random sample of 1000 simulated events, which are printed out.

     #include <stdio.h>
     #include <gsl_histogram2d.h>

     int
     main ()
     {
       gsl_histogram2d * h = gsl_histogram2d_calloc_uniform (10, 10,
                                                             0, 1, 0, 1) ;

       gsl_histogram2d_accumulate (h, 0.3, 0.3, 1) ;
       gsl_histogram2d_accumulate (h, 0.8, 0.1, 5) ;
       gsl_histogram2d_accumulate (h, 0.7, 0.9, 0.5) ;

       {
         int i ;
         gsl_histogram2d_pdf * p = gsl_histogram2d_pdf_alloc (h) ;

         for (i = 0 ; i < 1000 ; i++) {
           double x, y ;
           double u = ((double) rand ()) / RAND_MAX;
           double v = ((double) rand ()) / RAND_MAX;

           int status = gsl_histogram2d_pdf_sample (p, u, v, &x, &y) ;

           printf("%g %g\n", x, y) ;
         }
       }

      return 0 ;
     }

The following plot shows the distribution of the simulated events. Using a higher resolution grid we can see the original underlying histogram and also the statistical fluctuations caused by the events being uniformly distributed over the the area of these original bins.