[Bug-gsl] Bug in gsl_ran_exponential?

[Dr. Hans Ekkehard Plesser]

Hi!

I wonder whether there is a slight bug in gsl_ran_exponential():

double
gsl_ran_exponential (const gsl_rng * r, const double mu)
{
  double u = gsl_rng_uniform_pos (r);

  return -mu * log (u);
}

gsl_rng_uniform_pos(r) returns u with 

        0 < u < 1, 

so that  

        0 < -log(u) <= some big number.

This means that gsl_rng_uniform_pos cannot return 0, even though

        p(x) = 1/mu exp(-x/mu)

has its maximum p(x=0) = 1/mu at x=0.

In practice, the difference is probably negligible, since it affects only a 
single possible value of u.

A naive fix would be to use

  double u = 1.0 - gsl_rng_uniform (r);
   
instead. Since gsl_rng_uniform() returns a value in [0, 1), u would be in 
(0, 1]. But the above would map all values of gsl_rng_uniform() < machine 
precision to u = 1, i.e. return value 0, so this fix probably does more 
damage than good.

Best regards,
Hans E Plesser

-- 
Dr. Hans Ekkehard Plesser
Associate Professor

Dept. of Mathematical Sciences and Technology
Norwegian University of Life Sciences

Phone +47 6496 5467
Fax   +47 6496 5401
Email address@hidden
Home  http://arken.umb.no/~plesser