26 ms·

The claim here is not about expectations. random(X) is identically distributed to random() * X, even if X itself is a random variable, because the randomness i

The claim here is not about expectations. random(X) is identically distributed to random() * X, even if X itself is a random variable, because the randomness in random() is independent from the randomness in X. Indeed, the way you would implement random(X), if you only had random() and X available, would be to compute random() * X. So in particular, random(random()) is identically distributed to random() * random(). (To be clear, this is different from random()^2, as pointed out elsewhere.)

jasomill 3y ago
As a quick sanity check, Mathematica 13.2.1 Kernel for Linux ARM (64-bit) Copyright 1988-2023 Wolfram Research, Inc. In[1]:= DI = Distributed; In[2]:= EV = Expectation; In[3]:= SS = Subscript; In[4]:= UD = UniformDistribution; In[5]:= f[1] := EV[x, DI[x, UD[{0, 1}]]] In[6]:= f[n_Integer] := EV[x, DI[x, UD[{0, f[n - 1]}]]] In[7]:= Block[{n = 42}, EV[Product[SS[x, k], {k, 1, n}], Map[DI[SS[x, #], UD[{0, 1}]] &, Range[1, n]]] == Product[EV[x, DI[x, UD[{0, 1}]]], {k, 1, n}] == f[n]] Out[7]= True where n = 42 ∈ ℕ is arbitrary. Alas, Mathematica doesn't appear to support actual nested distributions of the form UD[{1,UD[{1,...}]}], so, e.g., f[2] == EV[x_2, DI[x_2, UD[{0, EV[x_1, DI[x_1, UD[{0,1}]]]}]]] == EV[x_2, DI[x_2, UD[{0,1/2}] == 1/4 which is trivial.
- jasomill 3y ago
  (As is the expectation of a uniform distribution in the first place. From the PDF for the uniform distribution (1/(b - a) for x ∈ [a,b], zero elsewhere) and the definition of EV as the integral of the product of the identity map and the PDF over the reals, EV[x, DI[x, UD[{a,b}]]] == Integrate[x * Piecewise[{ {1/(b - a), a <= x <= b}, {0, x < a || x > b}}], {x, -∞, ∞}] == Integrate[x * 1/(b - a), {x, a, b}] == ((b^2/2) - (a^2/2)) / (b - a) == (a + b)/2 so for a = 0, == b/2 no matter how many nested integrals it takes to define the real number b.)
SideQuark 3y ago
Agreed, in this case they agree, but it's completely a numerical coincidence and not general. So people with a rough feeling that they compute rand * 5 should equate to rand(rand())=rand * rand should be admonished to check it extremely carefully. And the result is not true at the level of code. In practice it's not true, unless extreme care is taken to ensure perfect uniformity, which is rare due to nonuniformity of IEEE floats. I've not seen a codebase in the wild that makes all the required pieces work. It's doable, I've written articles on it, but it's costly to do. Another error in practice is that the method of taking a rand * 5 is not uniform due to numerical representations. So these are not identical in practice. This the best way to analyze, which is trivially analyzable, is the nested form. It's also why I wrote "may" in the previous comment.
- anderskaseorg 3y ago
  It is true, both in the ideal mathematical sense, and also in the IEEE sense under the totally obvious implementation that most programmers would write without any special care: function rand(endpoint = 1) { return endpoint * Math.random(); } Here rand(rand()) would evaluate to literally the same number as rand() * rand() given the same state of the underlying generator. Someone who follows this chain of reasoning is perfectly correct to consider it intuitively obvious. And rand() * rand() is a great form for analysis, since it allows us to take advantage of E[XY] = E[X]E[Y] for independent X, Y to compute the mean very quickly. Clearly you followed a different chain of reasoning, and that’s fine; maybe you have a much more complicated generator in mind that takes more care to ensure ideal rounding at the limits of floating-point precision, and that’s fine too. But let’s not admonish anyone for making correct statements.
  - SideQuark 3y ago
    >It is true, both in the ideal mathematical sense, and also in the IEEE sense under the totally obvious implementation ..... That implementation is not uniform under IEEE 754. Here [1] is but one of hundreds of papers on the subject. To quote "Drawing a floating-point number uniformly at random from an interval [a, b) is usually performed by a location-scale transformation of some floating-point number drawn uniformly from [0, 1). Due to the weak properties of floating-point arithmetic, such a transformation cannot ensure respect of the bounds, uniformity or spatial equidistributivity.". This is trivial to understand. Suppose you truly had a IEEE uniform [0,1) in the sense every possible float representation could occur and an interval of [a,b) (for representable a,b) is represented with probability b-a. Then if you multiply, say, by 0.99, then by the pigenhole principle some of the original values would coalesce, and some would not, so now you are no longer uniform. Some values are more likely to occur. This also happens when you multiply by 5, or 7, or any number, except the ranges, as they stretch, due to finite precision, do not map to the correct length ranges. The new range [0,5) has a different set of representable floats, and the precision of them moves in predictable but non-uniform ways, so the result is no longer uniform. So it is not true that rand() * N is uniform for IEEE 754 numbers. Another simple way to see it: when you multiply by floats > 1, then some of the low bits get lost for some floats, so you cannot begin to hope the results are uniform. You're truncating/rounding, which loses information, so you by necessity only have an approximation. And note that very few (if any) popular languages even manage to get [0,1) uniform (more details on 15 languages are in [1]). Their conclusion, Table II, to the question of whether the language rng provides uniformity (and spatial equidistributivity, needed to scale as the above), is no for all 15 languages. And this is because most do what you claim works. Even getting the [0,1) uniform as a start is extremely tricky and nuanced. The usual "uniform random int in [0,2^p-1] then divide by 2^p (or 2^p-1)" fails to be uniform most of the time and (as the paper shows) for pretty much every language implementation out there. This issue is important for places needing super careful accuracy like physics sims, weather sims, many monte carlo algorithms where you don't want bad numerics screwing with your results, and in all of those places, they use custom, properly written methods to create such numbers. The general idea is that for getting a uniform in [a,b) you need to sample enough bits for the range, then carefully bit construct the floating point value - it cannot be transformed (easily) from smaller or larger ranges. Table V lists various methods to draw floats in [a,b). Your method is listed (as I claimed above) as not being uniform (in fact, it fails all 5 criteria in the table, the only method of the 7 listed to score so poorly). Here's [2] a paper showing the common method of using a PRNG then dividing by some number to get "uniform" in [0,1) is not uniform, which is well known [2]. So most implementations don't even start with uniform. This is a well known problem in random number generation. And the rand(rand()) == rand() * rand() relies on uniform distributions. So please stop claiming things until you have checked them carefully. Fuzzy feelings about what you guess is true is no substitute for actually doing the work to check, and doubling down when someone points you on the right path is bizarre to me. [1] https://hal.science/hal-03282794v4/document https://hal.science/hal-03282794v4/document [2] https://hal.science/hal-02427338/document\\ https://hal.science/hal-02427338/document\\*

jasomill 3y ago

As a quick sanity check, Mathematica 13.2.1 Kernel for Linux ARM (64-bit) Copyright 1988-2023 Wolfram Research, Inc. In[1]:= DI = Distributed; In[2]:= EV = Expectation; In[3]:= SS = Subscript; In[4]:= UD = UniformDistribution; In[5]:= f[1] := EV[x, DI[x, UD[{0, 1}]]] In[6]:= f[n_Integer] := EV[x, DI[x, UD[{0, f[n - 1]}]]] In[7]:= Block[{n = 42}, EV[Product[SS[x, k], {k, 1, n}], Map[DI[SS[x, #], UD[{0, 1}]] &, Range[1, n]]] == Product[EV[x, DI[x, UD[{0, 1}]]], {k, 1, n}] == f[n]] Out[7]= True where n = 42 ∈ ℕ is arbitrary. Alas, Mathematica doesn't appear to support actual nested distributions of the form UD[{1,UD[{1,...}]}], so, e.g., f[2] == EV[x_2, DI[x_2, UD[{0, EV[x_1, DI[x_1, UD[{0,1}]]]}]]] == EV[x_2, DI[x_2, UD[{0,1/2}] == 1/4 which is trivial.

jasomill 3y ago
(As is the expectation of a uniform distribution in the first place. From the PDF for the uniform distribution (1/(b - a) for x ∈ [a,b], zero elsewhere) and the definition of EV as the integral of the product of the identity map and the PDF over the reals, EV[x, DI[x, UD[{a,b}]]] == Integrate[x * Piecewise[{ {1/(b - a), a <= x <= b}, {0, x < a || x > b}}], {x, -∞, ∞}] == Integrate[x * 1/(b - a), {x, a, b}] == ((b^2/2) - (a^2/2)) / (b - a) == (a + b)/2 so for a = 0, == b/2 no matter how many nested integrals it takes to define the real number b.)

SideQuark 3y ago

Agreed, in this case they agree, but it's completely a numerical coincidence and not general. So people with a rough feeling that they compute rand * 5 should equate to rand(rand())=rand * rand should be admonished to check it extremely carefully. And the result is not true at the level of code. In practice it's not true, unless extreme care is taken to ensure perfect uniformity, which is rare due to nonuniformity of IEEE floats. I've not seen a codebase in the wild that makes all the required pieces work. It's doable, I've written articles on it, but it's costly to do. Another error in practice is that the method of taking a rand * 5 is not uniform due to numerical representations. So these are not identical in practice. This the best way to analyze, which is trivially analyzable, is the nested form. It's also why I wrote "may" in the previous comment.

anderskaseorg 3y ago
It is true, both in the ideal mathematical sense, and also in the IEEE sense under the totally obvious implementation that most programmers would write without any special care: function rand(endpoint = 1) { return endpoint * Math.random(); } Here rand(rand()) would evaluate to literally the same number as rand() * rand() given the same state of the underlying generator. Someone who follows this chain of reasoning is perfectly correct to consider it intuitively obvious. And rand() * rand() is a great form for analysis, since it allows us to take advantage of E[XY] = E[X]E[Y] for independent X, Y to compute the mean very quickly. Clearly you followed a different chain of reasoning, and that’s fine; maybe you have a much more complicated generator in mind that takes more care to ensure ideal rounding at the limits of floating-point precision, and that’s fine too. But let’s not admonish anyone for making correct statements.
- SideQuark 3y ago
  >It is true, both in the ideal mathematical sense, and also in the IEEE sense under the totally obvious implementation ..... That implementation is not uniform under IEEE 754. Here [1] is but one of hundreds of papers on the subject. To quote "Drawing a floating-point number uniformly at random from an interval [a, b) is usually performed by a location-scale transformation of some floating-point number drawn uniformly from [0, 1). Due to the weak properties of floating-point arithmetic, such a transformation cannot ensure respect of the bounds, uniformity or spatial equidistributivity.". This is trivial to understand. Suppose you truly had a IEEE uniform [0,1) in the sense every possible float representation could occur and an interval of [a,b) (for representable a,b) is represented with probability b-a. Then if you multiply, say, by 0.99, then by the pigenhole principle some of the original values would coalesce, and some would not, so now you are no longer uniform. Some values are more likely to occur. This also happens when you multiply by 5, or 7, or any number, except the ranges, as they stretch, due to finite precision, do not map to the correct length ranges. The new range [0,5) has a different set of representable floats, and the precision of them moves in predictable but non-uniform ways, so the result is no longer uniform. So it is not true that rand() * N is uniform for IEEE 754 numbers. Another simple way to see it: when you multiply by floats > 1, then some of the low bits get lost for some floats, so you cannot begin to hope the results are uniform. You're truncating/rounding, which loses information, so you by necessity only have an approximation. And note that very few (if any) popular languages even manage to get [0,1) uniform (more details on 15 languages are in [1]). Their conclusion, Table II, to the question of whether the language rng provides uniformity (and spatial equidistributivity, needed to scale as the above), is no for all 15 languages. And this is because most do what you claim works. Even getting the [0,1) uniform as a start is extremely tricky and nuanced. The usual "uniform random int in [0,2^p-1] then divide by 2^p (or 2^p-1)" fails to be uniform most of the time and (as the paper shows) for pretty much every language implementation out there. This issue is important for places needing super careful accuracy like physics sims, weather sims, many monte carlo algorithms where you don't want bad numerics screwing with your results, and in all of those places, they use custom, properly written methods to create such numbers. The general idea is that for getting a uniform in [a,b) you need to sample enough bits for the range, then carefully bit construct the floating point value - it cannot be transformed (easily) from smaller or larger ranges. Table V lists various methods to draw floats in [a,b). Your method is listed (as I claimed above) as not being uniform (in fact, it fails all 5 criteria in the table, the only method of the 7 listed to score so poorly). Here's [2] a paper showing the common method of using a PRNG then dividing by some number to get "uniform" in [0,1) is not uniform, which is well known [2]. So most implementations don't even start with uniform. This is a well known problem in random number generation. And the rand(rand()) == rand() * rand() relies on uniform distributions. So please stop claiming things until you have checked them carefully. Fuzzy feelings about what you guess is true is no substitute for actually doing the work to check, and doubling down when someone points you on the right path is bizarre to me. [1] https://hal.science/hal-03282794v4/document https://hal.science/hal-03282794v4/document [2] https://hal.science/hal-02427338/document\\ https://hal.science/hal-02427338/document\\*