[jira] [Commented] (RNG-154) Avoid infinite samples in the ZigguratNormalisedGaussianSampler

["Alex Herbert (Jira)" <ji...@apache.org>]

    [ https://issues.apache.org/jira/browse/RNG-154?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=17377416#comment-17377416 ] 

Alex Herbert commented on RNG-154:
----------------------------------

h2. Other normalised Gaussian samplers in the library

Infinite samples can also be returned by the BoxMullerNormalizedGaussianSampler due to the use of Math.log(rng.nextDouble()).

The MarsagliaNormalizedGaussianSampler calls Math.log(r2) on a squared radius of a sample (x,y) inside the unit circle. This only occurs when r2 is above zero, thus the sampler eliminates the infinite extreme.

The ZigguratSampler.NormalisedGaussian does not use Math.log. It generates the tail via an exponential distribution that uses recursive calls that add the term 7.56 to repeat samples. It will create a stack overflow before ever reaching infinity.

 

The Box-Muller method could also be fixed to avoid infinite samples. In this case I would favour not doing so as the sampler is a reference implementation and unlikely to be used in practice.

 

> Avoid infinite samples in the ZigguratNormalisedGaussianSampler
> ---------------------------------------------------------------
>
>                 Key: RNG-154
>                 URL: https://issues.apache.org/jira/browse/RNG-154
>             Project: Commons RNG
>          Issue Type: Improvement
>          Components: sampling
>    Affects Versions: 1.3
>            Reporter: Alex Herbert
>            Priority: Minor
>
> The ZigguratNormalisedGaussianSampler uses the following method to sample the tail of the Gaussian:
> {code:java}
> // Note
> // R = 3.442619855899
> // 1 / R = 0.2905...
> double y;
> double x;
> do {
>     y = -Math.log(rng.nextDouble());
>     x = -Math.log(rng.nextDouble()) * ONE_OVER_R;
> } while (y + y < x * x);
> final double out = R + x;
> return hz > 0 ? out : -out;
> {code}
> In the unlikely event the RNG produces 0 then Math.log(0) is infinity. Two zeros in a row can result in an infinite sample for the tail. This is very unlikely but would be unexpected for a user of the library since a sample should be roughly within +/-3.
> Note that if the RNG is sampling from the 2^53 dyadic rationals in [0, 1) then the next value is:
> {code:java}
> Math.log(0x1.0p-53) == -36.7368005696771
> {code}
> Here the returned tail would be 3.44 + 36.7368005696771 / 3.44 = 14.11. This is very far from the extreme of infinity.
> To avoid infinity this can be fixed by:
>  1. Assuming the RNG is returning a value in [0, 1) and using Math.log(1.0 - rng.nextDouble())
>  2. Generating the double u from a long to ensure the value is in [0, 1) and using 1.0 - u.
>  



--
This message was sent by Atlassian Jira
(v8.3.4#803005)