I want a random Float64 in the interval (0,1].
I understand rand() returns a number in [0,1).
Should I just use 1-rand() ?
1 Like
is rand() not more like (0,1)? And how much does [0,1) matter vs (0,1]? For most numerical problems, I assume the distinct makes no different what-so-ever in practice.
That’s a nice and quick solution. You can also check for 1 and make it into a 0, eg
rand(Float64) |> x -> iszero(x) ? one(x) : x
You should not use 1.0 - rand().
The density of Float64s (how many occur within a subspan) varies within [0.0,1.0].
(source: Float)
(source: Luis R. Izquierdo and J. Gary Polhill: Is Your Model Susceptible to Floating-Point Errors?)
for a single random Float64 in (0.0, 1.0]
randfloat = rand()
if randfloat === 0.0
randfloat = 1.0
end
for n of them
zerosasones(x) = ifelse(x === 0.0, 1.0, x)
randfloats = zerosasones.(rand(n))
2 Likes
I am not sure how that is relevant, AFAIK rand accounts for that and the result is uniform.
The density is varying but rand() does sample from it uniformly as it actually samples from [1,2[ and subtracts 1.0.
This is equivalent to running 2.0-rand(Random.CloseOpen12()) which avoids this subtraction.
You can check it by running:
julia> function t(m1, m2)
x = 1.0 - rand(m1)
y = 2.0 - rand(m2, Random.CloseOpen12())
x == y
end
t (generic function with 2 methods)
julia> m1, m2 = MersenneTwister(100), MersenneTwister(100);
julia> all(t(m1,m2) for i in 1:10^8)
true
The key thing here is that on [1,2[ interval all Float64 are uniform. You can check that eps for any value in [1,2[ range is the same. as well as the difference between x and nextfloat(x) in this range.
In other words if you look at Luis R. Izquierdo and J. Gary Polhill: Is Your Model Susceptible to Floating-Point Errors? image [1,2[ just happens to be one subsegments that is uniform on it.
Actually on Float it is shown on this picture:
10 Likes
using StatsBase: skewness
rs01 = rand(n)
rs10 = 1.0 .- rs01
skewness(rs01) == -skewness(rs10)
You are demonstrating numerical error (after correcting to rs10 = 1.0 .- rs01, which makes it run).
I am demonstrating skewness in small random samples and that it swaps signs using 1.0 .- rand()
More to the point, a uniform random number generator returning values in [0.0, 1.0) is used internally in simulation code that requires random variates from some specific non-unform distribution, often also in [0.0, 1.0). Once that happens, using one minus random variate probably maps some distinct smaller numbers together as they cross power of 2 bounds.
What you show is important, but it is unrelated to random number generation. If you take any vector x (random or not) you will get the same result. What I understand @Tamas_Papp says is that skewness(1.0 .- x) == skewness(x) will be often false due to roundoff errors.
The swapping of signs in skewness is not a function of roundoff errors, it is a function of the mirroring.
This is correct - if you have a non-uniform distribution on [0,1[ then taking 1 minus it might hit the non uniformity problem.
No, as s(x) = -s(A-x) is an identity for s being skewness, and A a constant. What you are showing is that this does not hold in practice, because of numerical error. That is orthogonal to this discussion. Cf eg
julia> x = rand(10^9);
julia> A = 100
100
julia> sum(x) + A == sum(x .+ A)
false
and various similar examples.
I am referring to:
julia> using Random
julia> Random.seed!(100);
julia> x = rand(100000);
julia> skewness(1.0 .- x) == -skewness(x)
false
and you see that the sign is swapped but the values are not exactly opposite - they are only approximately opposite.
this is getting away from itself – try the same thing using Float16 rands, Float16(skewness(x)) and 1000 samples. I know that floating point errors accrue – no argument there. 
coding to the internal workings of rand() … kinda iffy for FAA system simulations
So you’re saying that if you take 1.0 - rand() then many of the small floats, close to zero, can never be selected, since eps is larger close to 1.0? Those small numbers just become inaccessible to the algorithm.
1 Like
Yes. You can theoretically get 0 and the next theoretically possible value is eps().