[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 |
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
- [Bug-gsl] Bug in gsl_ran_exponential?,
Dr. Hans Ekkehard Plesser <=