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

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

To avoid Math.log(0) requires generating a uniform deviate in (0, 1]. The 
smallest value from the 2^53 dyadic rationals is 2^-53.
h2. BoxMullerGaussianSampler

This create a radius using:
{code:java}
// The extreme tail of the sample is:
// y = 2^-53
// r = 8.57167 
final double r = Math.sqrt(-2 * Math.log(y));
{code}
The tail of the normalized distribution is thus limited to 8.57167.
h2. ZigguratNormalisedGaussianSampler

This creates two samples using Math.log. This is used to create the tail when 
2y >= x^2:
{code:java}
// Note
// R = 3.442619855899
// 1 / R = 0.2905...

double y;
double x;
do {
    // Note: The extreme value y from -Math.log(2^-53) is (to 4 sf):
    // y = 36.74
    // The largest value x where 2y < x^2 is false is sqrt(2*36.74):
    // x = 8.571
    // The extreme tail is:
    // out = +/- 12.01
   double u1 = ...
   double u2 = ...
    y = -Math.log(u1);
    x = -Math.log(u2) * ONE_OVER_R;
} while (y + y < x * x);

final double out = R + x;
{code}
This condition can be generated if u1 = 2^-53 and u2 = 1377 * 2^-53.

The tail of the normalized distribution is thus limited to 12.0141.

 

> 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)

Reply via email to