It's interesting that the second circle is a pretty uniform distribution. `1 - random(random())` must be approximately `sqrt(random())`, iirc that's how you correct for the error demonstrated in the first circle (https://youtube.com/watch?v=4y_nmpv-9lI).
Why are they similar? For a rough idea, random() is a random number between 0 and 1, which on average is 0.5; random(random()) is a random number between 0 and on average 0.5, which on average is 0.25; and so on.
I think that should even make it exactly uniform, but somehow it doesn't look like it. I may have missed a bit somewhere.
Product: https://math.stackexchange.com/questions/659254/product-dist...
Minimum: https://stats.stackexchange.com/questions/220/how-is-the-min...
I can't explain it either, but I also don't think there is a reason they should look the same.
random(random()) has an identical distribution to random()*random() (it may even behave identically for a given rng state), although this is different to Math.pow(random(),2) since in that case there's 100% correlation between both parts which makes the expected value product bigger.
random(random()) is also distributed equivalently to
x=random()
y=random()
while(y>=x) y=random()
return min(x,y)
Another intuitive rephrasing is that while both min(random(), random()) and random(random()) are guaranteed to be at most your first draw, with random(random()) your second draw is constrained to be within the interval of your first, whereas with min(random(), random()) there's no such constraint on the second draw. Thinking about P(x > y); x,y~U should then provide the illumination.
If you plot their difference, 1-random(random()) aka 1-random()*random() is empirically skewed slightly larger than 1-random()^2, but I don't know how to analyze it exactly analytically.
So, random(k) should just be k rand() and so random(random()) should be just random() random().
One really fun way to do this is to rewrite using standard exp/ln rules as
# three exponential random variables each with mean 1
k1 = -ln(rand())
k2 = -ln(rand())
k3 = -ln(rand())
# rand() * rand() * rand() written cryptically
exp(-(k1 + k2 + k3))
The characteristic function E[exp(itX)] is a Fourier transform of the PDF and for a sum of independent variables it becomes a product,
E[exp(it(X + Y))] = E[exp(itX) exp(itY)] = E[exp(itX)] E[exp(itY)]
for an exponential distribution you have
E[exp(itX)] = ∫{0 → ∞} dx exp(-x) exp(itx) = 1/(1 – it).
So the Fourier transform of the PDF of a sum of n geometric random variables is going to be 1/(1 — i t)^n, which is the characteristic function for a gamma-distributed random variable,[1] and so the corresponding PDF is actually
f_n(x) = x^{n-1}/(n-1)! exp(-x)
[this is not too hard to see above when you realize that the exponent can be factored out as a repeated derivative of 1/(1 — i t) and in Fourier space a derivative corresponds to multiplication by the frequency ... derivatives with respect to t become multiplications by x].
So now we know that the distribution of the product U_n = exp(-G_n) where G_n is a Γ(n, 1)-distributed random variable, we work backwards from the definition:
h_n(u) du = the probability that U_n is between u and u + du
in which case G_n is between -ln(u) and -ln(u + du) = -[ ln(u) + ln(1 + du/u) ] ~= -ln(u) - du/u.
so that's f_n(x) dx where x = -ln(u) and dx = du/u, so that's
h_n(u) = (-ln(u))^{n-1}/(n-1)!, for 0 < u < 1.
and if you wanted the CDF it'd be an incomplete gamma function I think. For n=1 you can see that this recapitulates a uniform distribution but even for n=2 you have -ln(u) being the probability density function which I don't think will give you that sqrt() stuff that you need to get exact uniformity.