# Output distribution of rand(Float32) and rand(Float64), thread 2

**URL:** <https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184>\
**Category:** Internals & Design\
**Tags:** float, random\
**Created:** [October 19, 2023, 1:10pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184 "2023-10-19T13:10:07Z")\
**Posts on this page:** 20\
**Page:** 10

<div class="post-metadata">

**Author:** ![foobar\_lv2](https://avatars.discourse-cdn.com/v4/letter/f/ee59a6/32.png) [@foobar\_lv2](https://discourse.julialang.org/u/foobar_lv2)\
**Post date:** [October 28, 2023, 8:57pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/182 "2023-10-28T20:57:54Z")

</div>

> [@StefanKarpinski](#):
>
> For those advocating sampling from (0, 1), please provide a proposed support set along with a justifiable uniform distribution criterion on that support set.

If there was no consideration of performance, and this was pure math, I would do the following:

Sample a real number in `[0,1]`, uniformly according to standard measure. Then round it down to the next float.

In other words: `P(rand() == x) = |[x, nextfloat(x)]|` for `x: x == round(x)`, and `P(rand() < p) = p` for `p: p == round(p)`.

We now have a distribution that includes zero and subnormals. However, for Float64, Float32 and BFloat16, the probability of hitting exact zero or a subnormal number is smaller than `ldexp(1.0, -120)`. This probability means “never in our finite universe”.

More precisely, it means “only if the random generator is badly seeded”. On the other hand, exact zeros and subnormals are liable to trigger bugs that can mess up worse than bad statistics – e.g. security bugs.

So my principled distribution-with-safety-belt has `P(rand() < p) = p` for `p: p == round(p) and p > minnormal`, where minnormal is the smallest positive non-subnormal float. Then we clamp `P(rand() < minnormal) = 0`. This safety-belt is statistically invisible in our finite universe, but saves the bacon if somebody is bad at seeding the RNG.

======================  
Ok, we now need to do compromises with performance. E.g. your @StefanKarpinski proposed `randf`.

I am saying that we should retain the clamping close to zero from the ideal-with-safety-belt distribution. Clamp at something small like e.g. `ldexp(1.0, -60)`. This clamping is no more or less principled than the choice how to distribute our subset of floats where we have the desired equality `P(rand() < p) = p`

======================

All these clamping arguments do not apply to Float16. This is because `nextfloat(Float16(0)) == Float16(6.0e-8)` is a pretty large number with a probability that is significant! Clamping that introduces real bias that can mess up real statistics!

Hence, I have no opinion on Float16.

======================

PS. I think rounding `Bernoulli(1.0e-100)` to `Bernoulli(1.0e-20)` is a much worse mistake to make than e.g. rounding `Bernoulli(1.0e-18)` to `Bernoulli(0.0)`.

---

<div class="post-metadata">

**Author:** ![adienes](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/adienes/32/37459_2.png) [@adienes](https://discourse.julialang.org/u/adienes)\
**Post date:** [October 28, 2023, 8:59pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/183 "2023-10-28T20:59:13Z")

</div>

> [@foobar\_lv2](#):
>
> Sample a real number in `[0,1]`, uniformly according to standard measure.

I have seen you say this multiple times, but like… how?? isn’t that the whole point of this thread, to come up with algorithms to sample uniformly from unit interval?

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [October 28, 2023, 9:03pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/184 "2023-10-28T21:03:02Z")

</div>

> [@adienes](#):
>
> but like… how??

He’s giving the thing we’re supposed to model (he says “and this was pure math”) and then describes the suggested compromises…

> [@foobar\_lv2](#):
>
> the probability of hitting exact zero or a subnormal number is smaller than `ldexp(1.0, -120)`. This probability means “never in our finite universe”.

I think this is the key. Yes, excluding 0 induces a bias, but the bias is “essentially unobservable” and I’m with @foobar_lv2 and also with him about Float16. You really can’t expect much from Float16 when it comes to accurately modeling low probability events.

---

<div class="post-metadata">

**Author:** ![mbauman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mbauman/32/31082_2.png) [@mbauman](https://discourse.julialang.org/u/mbauman)\
**Post date:** [October 28, 2023, 9:03pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/185 "2023-10-28T21:03:33Z")

</div>

> [@dlakelan](#):
>
> Computer floats are only an approximation to real numbers. That they break down when very fine resolution is important is well known.

This is a fundamental disagreement between our world views. Floats are not approximate. They exactly identify a real value.

The current `rand() < p` is **exactly correct** for all `p`s that are multiples of `eps(0.5)`.

---

<div class="post-metadata">

**Author:** ![StefanKarpinski](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stefankarpinski/32/24_2.png) [@StefanKarpinski](https://discourse.julialang.org/u/StefanKarpinski)\
**Post date:** [October 28, 2023, 9:03pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/186 "2023-10-28T21:03:36Z")

</div>

That is neither uniform under the P(rand() \< p) = p criterion since P(rand() \< min(S)) = 0 \< min(S) nor under the P(rand() ≤ p) = p criterion since P(rand() \< max(S)) = 1 \> max(S). More generally, if `_rand()` is uniform in $[0, 1)& then your rand() will have probability too small at the smallest value in the support and too big everywhere else. This means that doing `rand() < p` to get an event with probability p is always wrong and not fixable in any obvious way.

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [October 28, 2023, 9:05pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/187 "2023-10-28T21:05:14Z")

</div>

> [@foobar\_lv2](#):
>
> PS. I think rounding `Bernoulli(1.0e-100)` to `Bernoulli(1.0e-20)` is a much worse mistake to make than e.g. rounding `Bernoulli(1.0e-18)` to `Bernoulli(0.0)`.

> [@mbauman](#):
>
> The current `rand() < p` is **exactly correct**

Except that no computer RNGs have that level of accuracy/uniformity in bit generation or even close, and as @foobar_lv2 says, if you ask for a Bernoulli(1f-50) you’ll get a Bernoulli(0.0)

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [October 28, 2023, 9:17pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/188 "2023-10-28T21:17:59Z")

</div>

More to the point, what’s the observable consequences of cutting off an interval between [0.0,nextfloat(0.0)] or [0.0f0,nextfloat(0.0f0)]?

Those are of size:

```julia
julia> nextfloat(0.0)
5.0e-324

julia> nextfloat(0.0f0)
1.0f-45

```

The first one is unobservable in our finite universe as @foobar_lv2 says, the second one might be observable… lemme think, a cluster of 1 million cores each operating at 1 billion operations per second would take on average …

```julia
julia> 1.0/nextfloat(0.0f0)/(1e9*1e6*3600*24*365)
2.2628863722506974e22

```

So 10^22 years to see the first consequence.

Again, unobservable in our universe.

However, computers aren’t perfect, bugs in RNG implementations or seeding could induce bias towards 0.0 dramatically more than they should have, so we’re talking about protecting against a known special case which will only be hit if the RNG is already “broken” in some sense. Which is the more important source of bias? the risk of broken software or the risk that you’ll miss out on 1 event in 10^22 years of massive supercomputing?

---

<div class="post-metadata">

**Author:** ![foobar\_lv2](https://avatars.discourse-cdn.com/v4/letter/f/ee59a6/32.png) [@foobar\_lv2](https://discourse.julialang.org/u/foobar_lv2)\
**Post date:** [October 28, 2023, 9:20pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/189 "2023-10-28T21:20:16Z")

</div>

> [@adienes](#):
>
> I have seen you say this multiple times, but like… how?? isn’t that the whole point of this thread, to come up with algorithms to sample uniformly from unit interval?

This thread has multiple points. We at first need to reach consensus what “sample uniformly from unit interval” means, mathematically speaking, in the context of floating point numbers. I had hoped that this consensus was reached by now?

Hence I answered with a mathematical definition, not with an algorithm.

Then we need to come up with fast algorithms to do that, approximately (there will be trade-offs between perf and precision).

And we may want appropriate safety rails. An important point I tried to make is that users are incapable of dealing with precise zeros in the rand output. I brought empirical evidence to the table.

This strongly suggests that the following property of the “real” distribution should not be traded away for performance: Exact zeros don’t happen with the computational resources available to mankind (unless the RNG is ill-seeded).

I was a little hyperbolic on that – `ldexp(1.0, -1000)` is cosmologically impossible, but `ldexp(1.0, -120)` is merely humanly impossible.

---

<div class="post-metadata">

**Author:** ![StefanKarpinski](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stefankarpinski/32/24_2.png) [@StefanKarpinski](https://discourse.julialang.org/u/StefanKarpinski)\
**Post date:** [October 28, 2023, 9:21pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/190 "2023-10-28T21:21:53Z")

</div>

> [@foobar\_lv2](#):
>
> All these clamping arguments do not apply to Float16. This is because `nextfloat(Float16(0)) == Float16(6.0e-8)` is a pretty large number with a probability that is significant! Clamping that introduces real bias that can mess up real statistics!
> 
> Hence, I have no opinion on Float16.

Wait, so the stance is that we should exclude zero for some float types but not others? That’s so much worse.

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [October 28, 2023, 9:23pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/191 "2023-10-28T21:23:31Z")

</div>

> [@StefanKarpinski](#):
>
> Wait, so the stance is that we should exclude zero for some float types but not others? That’s so much worse.

IMHO we shouldn’t by default provide rand for Float16’s but instead provide a function that errors with a message about please define rand(Float16) yourself, according to whichever types of biases you prefer and here’s a palette of several implementations you could choose from.

---

<div class="post-metadata">

**Author:** ![StefanKarpinski](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stefankarpinski/32/24_2.png) [@StefanKarpinski](https://discourse.julialang.org/u/StefanKarpinski)\
**Post date:** [October 28, 2023, 9:28pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/192 "2023-10-28T21:28:05Z")

</div>

Are you trolling or serious?

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [October 28, 2023, 9:28pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/193 "2023-10-28T21:28:57Z")

</div>

No I’m serious, the tradeoffs should be explicit and only a small number of people should be using Float16 I think? Or maybe I am not aware of people using Float16 like crazy in some current applications?

Or let me put it this way, the probability to get Inf in float16 sampling of Exponentials is… 6e-8 so “10 times a second”

---

<div class="post-metadata">

**Author:** ![mbauman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mbauman/32/31082_2.png) [@mbauman](https://discourse.julialang.org/u/mbauman)\
**Post date:** [October 28, 2023, 9:37pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/194 "2023-10-28T21:37:07Z")

</div>

You and @foobar_lv2 are continuing to conflate the number of available bit patterns with the values themselves. Ironically, those two things are only the same in the [1, 2) interval.

Once you have floating point numbers that have varying epsilons, you must be more considered in what it means to draw one such number.

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [October 28, 2023, 9:41pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/195 "2023-10-28T21:41:21Z")

</div>

> [@mbauman](#):
>
> You and @foobar_lv2 are continuing to conflate the number of available bit patterns with the values themselves

No, the size of the region that rounds to 0 is | [0,6e-8] | = 6e-8 (or within a factor of 2 of that, depending on rounding mode I guess)

---

<div class="post-metadata">

**Author:** ![Sukera](https://avatars.discourse-cdn.com/v4/letter/s/ce7236/32.png) [@Sukera](https://discourse.julialang.org/u/Sukera)\
**Post date:** [October 28, 2023, 9:42pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/196 "2023-10-28T21:42:12Z")

</div>

> [@dlakelan](#):
>
> Or let me put it this way, the probability to get Inf in float16 sampling of Exponentials is… 6e-8 so “10 times a second”

I can’t seem to reproduce that:

```julia-repl
julia> using Dates, Random

julia> try
       println(now())
       while true
           isinf(randexp(Float16)) && println("gotcha")
       end
       finally
       println(now())
       end
2023-10-28T23:41:26.450
^C2023-10-28T23:41:42.422
ERROR: InterruptException:
Stacktrace:
 [1] _randexp(rng::TaskLocalRNG, ri::UInt64)
   @ Random ~/julia/usr/share/julia/stdlib/v1.11/Random/src/normal.jl:129
 [2] randexp
   @ ~/julia/usr/share/julia/stdlib/v1.11/Random/src/normal.jl:127 [inlined]
 [3] randexp
   @ ~/julia/usr/share/julia/stdlib/v1.11/Random/src/normal.jl:198 [inlined]
 [4] randexp
   @ ~/julia/usr/share/julia/stdlib/v1.11/Random/src/normal.jl:199 [inlined]
 [5] top-level scope
   @ ./REPL[7]:4

```

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [October 28, 2023, 9:45pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/197 "2023-10-28T21:45:57Z")

</div>

> [@Sukera](#):
>
> I can’t seem to reproduce that:

Interesting. On the other hand if you ask it to print gotcha when rand(Float16) gives Float16(0.0) you get Gotcha 10 times a second. So evidently there’s a difference in the implementation for Exponential?

---

<div class="post-metadata">

**Author:** ![foobar\_lv2](https://avatars.discourse-cdn.com/v4/letter/f/ee59a6/32.png) [@foobar\_lv2](https://discourse.julialang.org/u/foobar_lv2)\
**Post date:** [October 28, 2023, 9:46pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/198 "2023-10-28T21:46:48Z")

</div>

> [@mbauman](#):
>
> The current `rand() < p` is **exactly correct** for all `p`s that are multiples of `eps(0.5)`.

> [@mbauman](#):
>
> You and @foobar_lv2 are continuing to conflate the number of available bit patterns with the values themselves. Ironically, those two things are only the same in the [1, 2) interval.

Are we posting in the same thread? Am I that bad at communicating?

```julia
julia> function foo(p, N)
       smaller = 0
       smallerEqual = 0
       for i=1:N
       r = rand(typeof(p))
       r <= p && (smallerEqual += 1)
       r < p && (smaller += 1)
       end
       @show N, p, smaller, smallerEqual, smaller/N, smallerEqual/N
       end
foo (generic function with 1 method)

julia> foo(ldexp(eps(0.5f0), -30), 1<<32)
(N, p, smaller, smallerEqual, smaller / N, smallerEqual / N) = (4294967296, 5.551115f-17, 276, 276, 6.426125764846802e-8, 6.426125764846802e-8)
(4294967296, 5.551115f-17, 276, 276, 6.426125764846802e-8, 6.426125764846802e-8)

```

This is what I’m complaining about. This result is wrong by 9 orders of magnitude.

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [October 28, 2023, 9:50pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/199 "2023-10-28T21:50:40Z")

</div>

> [@foobar\_lv2](#):
>
> This result is wrong by 9 orders of magnitude.

This literally made me laugh out loud.

---

<div class="post-metadata">

**Author:** ![Sukera](https://avatars.discourse-cdn.com/v4/letter/s/ce7236/32.png) [@Sukera](https://discourse.julialang.org/u/Sukera)\
**Post date:** [October 28, 2023, 9:53pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/200 "2023-10-28T21:53:24Z")

</div>

> [@foobar\_lv2](#):
>
> This is what I’m complaining about. This result is wrong by 9 orders of magnitude.

Yes, and as I pointed out in [message #37](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/38) of this thread, this is _expected_, because the current generation scheme does not allow for the precision you’re asking for.

Could it be improved? Maybe, if performance of such a scheme was good AND followed the currently guaranteed distribution of possibly generating `0.0`. So far, you’ve either generated a different distribution, generated a biased one, or had much worse performance.

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [October 28, 2023, 9:56pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/201 "2023-10-28T21:56:13Z")

</div>

> [@Sukera](#):
>
> generated a biased one

My point above is that generating a biased distribution when the bias is at the level of “global scale supercomputing can’t observe the bias in billions of times longer than the duration of the universe” is the right thing to do because _errors in computer code_ are much much more of a problem than biases of this size. From that perspective, excluding 0.0 is the right thing to do for Float64 and Float32 because it’s unobservable bias. (and note, it’s what R does when generating rexp)

Float16 is problematic as a model for probability.

[Previous page](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184.md?page=9)

[Next page](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184.md?page=11)