Go to the first, previous, next, last section, table of contents.
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.
A histogram is defined by the following struct,
size_t n
double * range
double * bin
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 <= 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'.
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
.
gsl_histogram
struct.
The bins are uniformly spaced with a total range of
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
.
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].
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.
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.
gsl_histogram_increment
but increases
the value of the appropriate bin in the histogram h by the
floating-point number weight.
GSL_EDOM
and the function returns 0.
GSL_EDOM
.
gsl_histogram
struct directly.
The following functions are used by the access and update routines to locate the bin which corresponds to a given x coordinate.
GSL_SUCCESS
. If x lies outside the valid range of the
histogram then the function returns GSL_EDOM
and the error
handler is invoked.
The library provides functions for reading and writing histograms to a file as binary data or formatted text.
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.
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.
%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.
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.
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 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.
size_t n
double * range
double * sum
The following functions allow you to create a gsl_histogram_pdf
struct which represents this probability distribution and generate
random samples from it.
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
.
s = range[i] + delta * (range[i+1] - range[i])
where i is the index which satisfies sum[i] <= r < sum[i+1] and delta is (r - sum[i])/(sum[i+1] - sum[i]).
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/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 ; double x; while (fscanf(stdin, "%lg", &x) == 1) { gsl_histogram_increment(h, x); } 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 Cauchy distribution with a width of 30 and histogram them over the range -100 to 100, using 200 bins.
gsl-randist 0 10000 cauchy 30 | gsl-histogram -100 100 200 > histogram.dat
A plot of the resulting histogram shows the familiar shape of the Cauchy distribution and the fluctuations caused by the finite sample size.
gnuplot> plot 'histogram.dat' using 1:3 with step
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).
Two dimensional histograms are defined by the following struct,
size_t nx, ny
double * xrange
double * yrange
double * bin
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'.
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
.
gsl_histogram2d
struct. The bins are uniformly spaced with a
total range of
0 <= x < nx in the x-direction and
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
.
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.
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.
gsl_histogram2d_increment
but increases
the value of the appropriate bin in the histogram h by the
floating-point number weight.
GSL_EDOM
and the function returns 0.
GSL_EDOM
.
gsl_histogram2d
struct directly.
The following functions are used by the access and update routines to locate the bin which corresponds to a given (x\,y) coordinate.
GSL_SUCCESS
. If
(x,y) lies outside the valid range of the histogram then the
function returns GSL_EDOM
and the error handler is invoked.
The library provides functions for reading and writing two dimensional histograms to a file as binary data or formatted text.
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.
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.
%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.
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.
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} is the number of events in the bin which contains (x,y), 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.
size_t nx, ny
double * xrange
double * yrange
double * 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.
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
.
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/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.
Go to the first, previous, next, last section, table of contents.