make `rand` unbiased for floats

[tpapp]

The current implementation of samplers for eg Float23 and Float64 effectively generates numbers as n*eps(T), where n is an integer between 0 and 1/eps(T)-1 (inclusive). Eg for Float64 this range is 0:(2^52-1).

This makes the random draws (slightly) biased, by -eps(T)/2.

The proposal is to add this term back to the result in the implementation.

Cf this discussion, especially the proposal of @perrutquist

[mschauer]

Just as a data point, to observe this bias, you would need ballpark

julia> Ts = [Float16, Float32, Float64]
julia> [Ts [(std(Uniform())/(eps(T)/2))^2 for T in Ts]]
3×2 Array{Any,2}:
 Float16  3.49525e5 
 Float32  2.34562e13
 Float64  6.7608e30 

number of trials.

[tpapp]

Possibly, but since we can do it right with a very simple correction, why shouldn't we?

[ViralBShah]

cc @rfourquet

[mschauer]

Some more checks (right column after adding eps(T)/2):

julia> T = Float16;[["P(rand() <= 2^(-$i)" for i in 0:16] [mean([(n)*eps(T) <= 2.0^(-i) for n in 0:(1/eps(T)-1)]) for i in 0:16] [mean([(n)*eps(T) + eps(T)/2  <= 2.0^(-i) for n in 0:(1/eps(T)-1)]) for i in 0:16] ]
17×3 Array{Any,2}:
 "P(rand() <= 2^(-0)"   1.0          1.0        
 "P(rand() <= 2^(-1)"   0.500977     0.5        
 "P(rand() <= 2^(-2)"   0.250977     0.25       
 "P(rand() <= 2^(-3)"   0.125977     0.125      
 "P(rand() <= 2^(-4)"   0.0634766    0.0625     
 "P(rand() <= 2^(-5)"   0.0322266    0.03125    
 "P(rand() <= 2^(-6)"   0.0166016    0.015625   
 "P(rand() <= 2^(-7)"   0.00878906   0.0078125  
 "P(rand() <= 2^(-8)"   0.00488281   0.00390625 
 "P(rand() <= 2^(-9)"   0.00292969   0.00195313 
 "P(rand() <= 2^(-10)"  0.00195313   0.000976563
 "P(rand() <= 2^(-11)"  0.000976563  0.000976563
 "P(rand() <= 2^(-12)"  0.000976563  0.0        
 "P(rand() <= 2^(-13)"  0.000976563  0.0  


julia> T = Float16;[["P(rand() < 2^(-$i)" for i in 0:16] [mean([(n)*eps(T) < 2.0^(-i) for n in 0:(1/eps(T)-1)]) for i in 0:16] [mean([(n)*eps(T) + eps(T)/2  < 2.0^(-i) for n in 0:(1/eps(T)-1)]) for i in 0:16] ]
17×3 Array{Any,2}:
 "P(rand() < 2^(-0)"   1.0          1.0        
 "P(rand() < 2^(-1)"   0.5          0.5        
 "P(rand() < 2^(-2)"   0.25         0.25       
 "P(rand() < 2^(-3)"   0.125        0.125      
 "P(rand() < 2^(-4)"   0.0625       0.0625     
 "P(rand() < 2^(-5)"   0.03125      0.03125    
 "P(rand() < 2^(-6)"   0.015625     0.015625   
 "P(rand() < 2^(-7)"   0.0078125    0.0078125  
 "P(rand() < 2^(-8)"   0.00390625   0.00390625 
 "P(rand() < 2^(-9)"   0.00195313   0.00195313 
 "P(rand() < 2^(-10)"  0.000976563  0.000976563
 "P(rand() < 2^(-11)"  0.000976563  0.0        
 "P(rand() < 2^(-12)"  0.000976563  0.0  

This doesn't look too bad at first sight.

[tpapp]

@mschauer: I am sorry, I don't understand the purpose of your execise.

The bias exists and can be described exactly. There is no need to simulate anything.

The bias is correctable, rather cheaply. Whether it "looks bad" in some statistical tests that may or may not have power for detecting it (it is a rather small bias, after all) is not the issue here.

The key question is whether there are any objections to correcting this bias on the table. Otherwise, I would just make a PR. @rfourquet, I would appreciate your opinion on this.

[mschauer]

I am sorry, I don't understand the purpose of your exercise.

A good implementation of rand() would for example have the property that the probabilities of rand() <= p and rand() < p should be equal to p. I checked*) this for our Random.rand(T) and for your Random.rand(T) + eps(T)/2 and concluded "it doesn't look to bad" for your idea because e.g. P(rand() <= 0.5) = 0.5 and not 0.500977. I am not sure why you are a bit defensive here. Why not discuss the consequences of this change?

*) Using the information "rand() effectively generates numbers as n*eps(T), where n is an integer between 0 and 1/eps(T)-1 (inclusive)" you gave above.

[tpapp]

Sorry if I came across as defensive, that was not the intention.

The consequences of this change are that rand will be unbiased for floats, including the bias you demonstrated, at no measurable performance loss (seems to be negligible compared to noise):

julia> using Random, BenchmarkTools

julia> rand_corrected(T) = rand(T) + eps(T)/2
rand_corrected (generic function with 1 method)

julia> Random.seed!(1);

julia> @btime foreach(_ -> rand(Float64), 1:1024)
  4.382 μs (0 allocations: 0 bytes)

julia> Random.seed!(1);

julia> @btime foreach(_ -> rand_corrected(Float64), 1:1024)
  4.368 μs (0 allocations: 0 bytes)

[perrutquist]

The performance loss wouldn't just be negligible, it would be exactly zero. The current implementation is something like rand(r, CloseOpen12()) - 1.0 and this could be replaced by rand(r, CloseOpen12()) - prevfloat(1.0)

I.e. the only change is in a compile-time constant.

[tpapp]

My understanding is that this rearrangement is valid because eps is 2^-52 in the [1, 2) range and 1.0-prevfloat(1.0) == 2^-53, but I don't know enough about floating point error to prove it formally. But I can of course make the PR do this.

[mschauer]

I can verify the bias.

[mschauer]

First the Float16 case: With

rnd(::Type{Float16}, i) = Float16(reinterpret(Float32,
                              (UInt32(i) << 13)  | 0x3f800000) - 1)
rndnew(::Type{Float16}, i) = Float16(reinterpret(Float32,
                              (UInt32(i) << 13)  | 0x3f800000) - 0.9995117f0)

being a deterministic model for rand(Float16) generation where i ranges in 0:0x000003ff, I get

julia> mean(rnd(Float16, i) for i in 0:0x000003ff)
Float16(0.4995)

julia> extrema(rnd.(Float16,  0:0x000003ff))
(Float16(0.0), Float16(0.999))

julia> mean(rndnew(Float16, i) for i in 0:0x000003ff)
Float16(0.5)

julia> extrema(rndnew.(Float16,  0:0x000003ff))
(Float16(0.0004883), Float16(0.9995))

[perrutquist]

Although floating-point addition/subtraction is not associative in general, it's quite easy to show that (r - 1.0) + eps()/2 == r - (1.0 - eps()/2) in this case (on systems that adhere to the IEEE standard). It follows directly from these two facts:

1. Addition and subtraction produce correctly rounded results.
2. The intermediate terms (r - 1.0), eps()/2 and (1.0 - eps()/2) can all be expressed exactly, i.e. are not affected by rounding.

[mschauer]

We do not need to add eps(Float16)/2, we need to subtract the right number from r where r in [1,2) such that the result is unbiased in [0,1], which is 0.9995117f0 in my code above. This coincides with julia> Float32((1.0 - eps(Float16)/2))
0.9995117f0, but we do not really need to worry about it.

[andreasnoack]

The problem discussed here appears to be of very small significance in practice. See #8958 (comment) and #8958 (comment)