Is it really possible for rand to throw a zero?

Specific Domains Statistics
61

To quote from the linked issue: The bitfiddling approach would be

randu64tof64(x::UInt64) = reinterpret(Float64, (x>>12) | (1022-trailing_zeros(x))<<52)

That has comparable speed to a single rand, and works because trailing_zeros is already exponentially distributed. (somewhat slower for MersenneTwister that natively supports generation of 52 bit random numbers, and harder to simd)

“Proper” draws are important whenever you really want a random float, instead of a random fixed-point number, i.e. when you care about the cool property that floats have higher precision close to zero (scaling invariance) especially for Float32. E.g. if you have a singularity close to zero like log(x), 1/x, 1/log1p(x), or if you want to start an iteration like x-> mod1(A*x, T) (where zero is a fixed point).

If you just want to add the numbers, then I don’t think the small bias is so relevant (addition destroys the extra precision close to zero anyways), but you get large biases on the tails if you use the obvious approach of sampling from e.g. an exponential distribution (obvious approach: Sample uniform in [0,1] and use log, i.e. inverse cumulative distribution to transform to [0, Inf]; works with “proper” draws, fails with “improper” draws).

1 Like
62

I like this, but I think randu64tof64 would have bias of -ϵ/8 (where ϵ = eps()) because while it would sometimes draw in the range below ϵ/4, it would never draw in the range above 1-ϵ/4, which is rounded to 1.0.

I can see why a “practical” random number generator might avoid returning 1.0, but a “proper” one should return it with probability ϵ/4 imo.

rand(UInt64) * (2^-64) would be almost as fast, and return 1.0 with the right probability.