On reviewing the code of numpy/random/mtrand/distributions.c, we noticed that in rk_standard_exponential, a trick is used to mitigate the fact that rk_double returns [0, 1).
|
double rk_standard_exponential(rk_state *state) |
|
{ |
|
/* We use -log(1-U) since U is [0, 1) */ |
|
return -log(1.0 - rk_double(state)); |
but why isn't the same trick used in rk_laplace?
|
double rk_laplace(rk_state *state, double loc, double scale) |
|
{ |
|
double U; |
|
|
|
U = rk_double(state); |
|
if (U < 0.5) |
|
{ |
|
U = loc + scale * log(U + U); |
|
} else |
|
{ |
|
U = loc - scale * log(2.0 - U - U); |
|
} |
|
return U; |
|
} |
This may lead to problem if U == 0 since log(U + U) will give -inf.
CC: @dingsquared
Numpy/Python version information:
Current master (8f31f95)
On reviewing the code of
numpy/random/mtrand/distributions.c, we noticed that inrk_standard_exponential, a trick is used to mitigate the fact thatrk_doublereturns [0, 1).numpy/numpy/random/mtrand/distributions.c
Lines 110 to 113 in 8f31f95
but why isn't the same trick used in rk_laplace?
numpy/numpy/random/mtrand/distributions.c
Lines 664 to 677 in 8f31f95
This may lead to problem if
U == 0sincelog(U + U)will give-inf.CC: @dingsquared
Numpy/Python version information:
Current
master(8f31f95)