The underlying thing here is that random(random()) in this case is the same as random() * random(). So random(random(random(...))) is the same as random() * random() * random() and then the analysis goes on. And sure, random() * random() has a mean close to 1/4. But the dynamics change as N becomes large.
Edit - and just in case you doubt whether random() * random() * ... is a valid way of doing this, I also just checked the laborious way of doing it:
>>> def foo(n):
... result = 1.0
... for _ in range(n):
... result = np.random.uniform(0, result)
... return result
...
>>> foo(100)
1.4267531652344414e-46
>>> foo(100)
7.852496730908136e-49
>>> foo(100)
1.3216070221780724e-41
is probably a good marker that it is time for me take my bow in this conversation, however, for an alternative approach I recommend SideQuark's comment on iterative random variables vs chains of multiplied random variables, which have different characteristics when defined as a sequence.
I checked his comment, I think he's incorrect on the approaches being different. Using my function `foo` that does the iterative approach, we can compare the distributions and they're fairly identical.
>>> np.mean([foo(100) for _ in range(100000)])
3.258425093913613e-33
>>> np.mean([np.product(np.random.random(100)) for _ in range(100000)])
8.814732867008917e-33
(There's quite a bit of variance on a log scale, so 3 vs 8 is not a huge difference. I re-ran this and got varying numbers with both approaches. But the iterated code is very slow...)
Note that the mean is actually quite a bit closer to 2^-100 even though the vast majority of numbers from either distribution fall below 2^-100. Even so, the mean for both is approximately a factor of 100 less than 2^-100. Suspicious! Although I think we've both burned enough time on this.
Yeah, the problem with your sampling experiment (as your suspicious observation makes it sound like you’ve started to realize) is that this distribution has a very narrow tail to 1 that you won’t hit by accident but nonetheless contributes very significantly to the computation of the mean.
You correctly computed that the mean of log(foo(n)) is exactly −n, and the median of foo(n) indeed shrinks like e^−n (more specifically like e^(1/3 − n), I think), but the mean of foo(n) is exactly 2^−n.
I don't think the difference will show up for small N. This is an asymptotes thing. Try it for N = 100, that's what I did. For example:
The underlying thing here is that random(random()) in this case is the same as random() * random(). So random(random(random(...))) is the same as random() * random() * random() and then the analysis goes on. And sure, random() * random() has a mean close to 1/4. But the dynamics change as N becomes large.Edit - and just in case you doubt whether random() * random() * ... is a valid way of doing this, I also just checked the laborious way of doing it: