Re: Bug: Laplace RNG formula seems to be incorrect

[Peter Johansson]

Hi Elias,

I'm just a user of GSL, not maintainer or contributor of this code, but

On 27/11/20 9:18 am, Elias Youssef H Netto wrote:
> Dear all,
>
> Thanks for your effort on this awesome package.
> It is great to have a fast, free implementation of several useful
> algorithms.
>
> Bug Report:
>
> Apparently, the formula for the Laplace RNG in the function gsl_ran_laplace
> (file /randist/laplace.c) seems to be wrong.
>
> See details below:
> - GSL version: 2.6.2
> - OS - Arch on Lenovo desktop
> - The compiler used - not relevant
> - A description of the bug behavior - see below
> - A short program which exercises the bug - see below
>
> My own calculations point to a formula that should be (see the attached
> pdf):
>
> $F^{-1}(y) &= u -sgn(y - 0.5) \cdot b \cdot log(1 - sgn(y - 0.5)(2y - 1))$
>
> Which agrees with what Wikipedia presents on its page (
> https://en.wikipedia.org/wiki/Laplace_distribution).
>
> The gsl_ran_laplace function is currently implemented as:
>
> ```C
> double
> gsl_ran_laplace (const gsl_rng * r, const double a)
> {
>    double u;
>    do
>      {
>        u = 2 * gsl_rng_uniform (r) - 1.0;
>      }
>    while (u == 0.0);
>
>    if (u < 0)
>      {
>        return a * log (-u);
>      }
>    else
>      {
>        return -a * log (u);
>      }
> }
> ```
> The problem is in the returns inside the if...else statement.
> It should be:
>
> ```C
> double
> gsl_ran_laplace (const gsl_rng * r, const double a)
> {
>    double u;
>    do
>      {
>        u = 2 * gsl_rng_uniform (r) - 1.0;
>      }
>    while (u == 0.0);
>
>    if (u < 0)
>      {
>        return a * log (1+u);
>      }
>    else
>      {
>        return -a * log (1-u);
>      }
> }
> ```


I don't see what the problem is. In the if (u<0.0) clause we have u 
uniformly distributed on [-1.0, 0.0) and the GSL implementation returns

return a * log (-u);

and you suggest

return a * log (1+u)

that is very similar because -u belongs to (0.0, 1.0] and 1+u belongs to 
[0.0, 1.0) so only difference is that you include 0.0 instead of 1.0. 
Returning log(0.0) seems problematic, while including log(1.0) = 0.0 in 
the possible return value seems perfectly legit.

Cheers,

Peter