[Bug-gsl] gsl_ran_beta can return 0 or 1

[Toby Johnson]

I'm not sure whether this is a bug.

gsl_ran_beta() can return values that are 0 or 1.  This can cause
problems if I then want to look up the pdf for the value just
generated by calling gsl_ran_beta_pdf, because the pdf can be infinite
at 0 or 1.  This sequence of events is often used in a
Metropolis--Hastings algorithm.  The following code generates this
`bug' on the 10th iteration, using the default generator and seed.  I
haven't tried to fix this since I'm not sure whether others will
consider it a bug.

#include <iostream>
using namespace std; 

#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

int main() {

  const gsl_rng_type *rng_type;
  long int rng_seed;
  gsl_rng *rng;
  
  gsl_rng_env_setup();
  rng_type = gsl_rng_default;
  rng_seed = gsl_rng_default_seed;
  rng = gsl_rng_alloc (rng_type);

  gsl_rng_set(rng,rng_seed);

  double x,newx,probx;

  x = 0.5;
  
  do {
    newx = gsl_ran_beta(rng,5.*x,5.*(1.-x));
    cout << x << " -> " << newx << endl;
    probx = gsl_ran_beta_pdf(newx,5.*x,5.*(1.-x));
    cout << "p = " << probx << endl;

    x = newx;
  }   while (1);

  return(0);
}

signature.asc
Description: Digital signature