Using gsl_histogram_pdf_sample to sample from an ad-hoc distribution

Question

I am having some trouble understanding how to use the function gsl_histogram_pdf_sample from the library GSL in C++. The documentation is here,

http://www.gnu.org/software/gsl/manual/html_node/The-histogram-probability-distribution-struct.html#The-histogram-probability-distribution-struct

I am not really an expert yet, so, I was wondering if anyone could tell me what is wrong with this code,

#include <iostream>
#include <gsl/gsl_histogram.h>
#include <gsl/gsl_rng.h>

using namespace std;

int main()
{
    // I am going to use 5 bins
    size_t  Bins = 5;
    // These are the ranges (must be Bins + 1)
    double range[6] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0 };
    // Array with probabilities
    double w[5]  = {0.05, 0.1, 0.3, 0.4, 1};

    // Create the histogram pdf

    gsl_histogram_pdf MyHistPdf;

    MyHistPdf.n     = Bins;
    MyHistPdf.range = range;
    MyHistPdf.sum   = w;

    const gsl_rng_type * T;
    gsl_rng * r;
    T = gsl_rng_default;
    r = gsl_rng_alloc (T);
    double u = gsl_rng_uniform(r);
    cout << u << endl;
    double a = gsl_histogram_pdf_sample(&MyHistPdf, u);

    return 0;
}

The program compiles with no errors, but when I run it, I always get the following error,

gsl: /usr/src/gsl-1.16-1/src/gsl-1.16/histogram/pdf.c:46: ERROR: cannot find r in cumulative pdf

And I have no idea what does it mean.

Answer 1

Before answering your question, you need to understand a fundamental rule about GSL: "Unless you understand GSL source code in detail, never manipulate gsl structs variables directly". If you do that, your program will inevitable fail or your results will make no sense! Why? Because you don't know in detail the internal gsl conventions about how those variables are storage and manipulated!

Read the documentation carefully, and you will see that there are API functions to allocate and manipulate all struct variables for you. If you are a C++ coder, a good "rule of thumb" is to imagine that all variables in GSL structs are private (the only reason this is not the case is because C does not have the capability of making member variables private!).

It is clear in the documentation that a gsl_histogram struct needs to be allocated, some data needs to be inserted and the gsl_histogram_pdf you use needs to be initialized using a specific API function before you can call gsl_histogram_pdf_sample. I saw that you tried to avoid inserting some data by hacking "sum" pointer and GSL failed, as I said it would happen if you don't follow the fundamental rule!

Here is an example that I wrote to guide you

#include <iostream>
#include <cassert>
#include <gsl/gsl_histogram.h>
#include <gsl/gsl_rng.h>

using namespace std;

int main()
{
   size_t  Bins = 5;
   double range[6] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0 };

    // histogram 
    gsl_histogram* hist = gsl_histogram_alloc (Bins);
    assert( hist != NULL );
    gsl_histogram_set_ranges (hist, range, 6);

    // include some elements
    const gsl_rng_type * T = gsl_rng_default;
    gsl_rng *r = gsl_rng_alloc (T);

    for(int j = 0; j < 10000; ++j ) {
      const double u = gsl_rng_uniform(r) * 5.0; 
      gsl_histogram_increment (hist, u);            
    }

    // Create the histogram pdf
    gsl_histogram_pdf*  MyHistPdf = gsl_histogram_pdf_alloc (Bins);
    assert( MyHistPdf != NULL );
    int status = gsl_histogram_pdf_init (MyHistPdf, hist);
    assert( status != GSL_EDOM );

    std::cout << "result you want";
    std::cout << gsl_histogram_pdf_sample(MyHistPdf,gsl_rng_uniform(r)) << std::endl;           
    return 0;
  }