Random number in (0,1)

Specific Domains Statistics
1

In Random.jl, there is the possibility to generate a random number in the interval [0,1) through CloseOpen01, or even in [1,2) (I can see the reason for this for implementation, less so for a user) through CloseOpen12, but there is no way to specify OpenOpen01.

Generating in (0,1) is particularly important to simulate various distributions through the inverse cumulative distribution function, a very common technique in Monte-Carlo methods. The dSFMT code provides such methods. Of course, it is always possible to add an if statement around each random number generation, but this is not particularly clean or efficient, plus it will not match the original dSFMT numbers in (0,1).

1 Like
2

I just want to add question, if someone is able to answer you: How important this can be in practice? With a period of \sim 2^{20000}, an implementation with open or closed intervals makes any difference in any actual application?

(I am assuming that this is a question addressed to the generation of a sequence of “continuous” floats, if one is sampling a discrete distribution that might be a problem, but then we only need to adjust the interval of the sequence generated).

Anyway, a simple workaround would be rnd = nextfloat(0.) + rand() (which actually does not affect the upper limit, because 1. + nextfloat(0.) == 1.)

6 Likes
3

There is a solution by @Elrod here to get the random numbers distribution over (0,1):

prevfloat(1.0)*(1-rand())
5 Likes
4

That will still skew the result slightly, as it unproportionally increases the likelyhood of 2^{-2021}, but that’s probably fine for most practical uses.

2 Likes
5

I think that other thread exhaust the subject. I was thinking only that one might want to avoid 0. because of numerical instabilities (even if very, very rare), but even for that reason changing that does not make any difference, because

julia> 1. / nextfloat(0.) == 1. / 0. == Inf
true
6

Currently, I’m subtracting the following from the [1,2) random number in VectorizedRNG.jl:

oneopenconst(::Type{Float64}) = 0.9999999999999999
oneopenconst(::Type{Float32}) = 0.99999994f0

This provides:

julia> 1 - oneopenconst(Float64)
1.1102230246251565e-16

julia> 1 - oneopenconst(Float32)
5.9604645f-8

julia> prevfloat(2.0) - oneopenconst(Float64)
0.9999999999999999

julia> prevfloat(2f0) - oneopenconst(Float32)
0.99999994f0
5 Likes
7

I could see many use cases where this nextfloat(0) will be problematic. The Inverse CDF of this number is just typically huge and will likely strongly skew results of a Monte-Carlo simulation

8

This will clearly break any Monte-Carlo simulation relying on an inverse cdf as invcdf(0) = -Inf.
The result of the simulation will then be Inf or NaN.

9

The first float after 0. is 5.0e-324, while you get 1/x = Inf for numbers up to about 1e-309. So the problem, if it is there, is not only about the open or closed intervals, but defining the inverse of these number in a safe way.

3 Likes
10

@Elrod I like your idea best. I suppose we loose a tiny bit of accuracy.
It also does not break any symmetry and is fast:

1-((2-eps(Float64)) - oneopenconst(Float64)) = 1 - oneopenconst(Float64)