A better title would be "1 - Random(Random(Random(Random())))", than it's more understandable why the points lean to the circumference.

Yes exactly. I was wondering why the center was becoming more and more empty when the result was moving more to 0.

... That means each point is at a distance of sqrt(f(x)) from the circumference. And nesting random() n times is similar to random()^n, which explains why the first plot has them amassed in the center and the second has a more uniform distribution.

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.

Oh man if my math teachers taught math with application, I would have been so much better off. It's so much easier to understand all of this after writing code for a while.

The cumulative distribution of random(random(...)) repeated k times is 1-(1-x)^k, so 1 - random(random()) should have a cumulative distribution of x^(2).

I think that should even make it exactly uniform, but somehow it doesn't look like it. I may have missed a bit somewhere.

Perhaps it looks nonuniform because we're not using literal points (our dots have nonzero area), so when they get too close they blend together and make the distribution look denser than it is in that area

Could be, but it is also possible I'm missing something somewhere, especially because I can't explain why min(random(), random()) looks different from random(random()) or random()*random().

According to the internet what you wrote is the distribution of min(random(), random()), but the distribution of the product is different.

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.

An attempt to explain the difference without just calculating the things:

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)

(the last line could also read 'return y') Comparing that to min(random(),random()), we can see that if the second call to random is smaller than the first, they will return the same result; otherwise, the program equivalent to random(random()) will return a smaller value, therefore the expected value of random(random()) must be lower that that of min(random(),random()).

Ah I think I see it now, yes given X and Y with uniform distribution then the distribution of X is uniform given Y and X < Y, but the distribution of min(X,Y) is not uniform given max(X,Y).

Nice argument! And to clarify further, the equivalent implementation for `random(random())` could be viewed as a rejection-sampling implementation of random(random()).

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.

1-random(random()) aka 1-random()*random(),looks slightly skewed toward circumference, and 1-random()^2 looks slightly skewed toward center.

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.

Yeah you're wrong :)

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))

That these are exponential random variables comes from the inverse-CDF Sampling Theorem.

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.

[1] https://en.wikipedia.org/wiki/Gamma_distribution .