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 [EMAIL PROTECTED]
Home http://arken.umb.no/~plesser
_______________________________________________
Bug-gsl mailing list
[email protected]
http://lists.gnu.org/mailman/listinfo/bug-gsl
