[Top][All Lists]

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[Bug-gsl] Bug in gsl_ran_exponential?

From: Dr. Hans Ekkehard Plesser
Subject: [Bug-gsl] Bug in gsl_ran_exponential?
Date: Wed, 16 Aug 2006 17:08:12 +0200
User-agent: KMail/1.8.2


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

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

reply via email to

[Prev in Thread] Current Thread [Next in Thread]