# Average number of uniform variates needed for their sum to be $> 1$ should be $\approx e$

An interesting tweet from Fermat's Library proves that the average number of random values [0,1] that are needed to Total > 1 is about $e$ (2.718...).

I don't get the proof, but I thought that a million trials in Mathematica should show something close to $e$ but I keep getting about 2.58 - I was hoping it would be an example of how you could get a clue about the proof by experimenting in Mathematica.

I suspect my coding skills are at fault, not the proof...

f[n_] := If[Total@RandomReal[1, n] > 1, n, f[n + 1]]
Table[f[2], 1000000] // Mean // N

• When f is invoked recursively, it doesn't know what random numbers have been chosen and summed so far. For example, try something like N[Mean[Table[count = 0; sum = 0; While[sum <= 1, sum += RandomReal[]; count++]; count, 1000000]]] – ilian Oct 30 '17 at 3:49
• Here is the CV thread that explains the result. – J. M.'s technical difficulties Oct 31 '17 at 6:58

## 3 Answers

Mean[Length/@Table[NestWhile[Append[#,RandomReal[]]&,{RandomReal[]},Total@# < 1 &], 1000]]//N

• You could also use Length@NestWhileList[# + RandomReal[] &, RandomReal[], # < 1 &]. It's more efficient (although in this case it really doesn't matter). – b3m2a1 Oct 30 '17 at 4:34

Here is a very simple recursive solution.

counter[cnt_, total_] /; total > 1. := cnt
counter[cnt_, total_] := counter[cnt + 1, total + RandomReal[]]
Mean[Table[counter[0, 0], 1000000]] // N


2.71805

• Far more readable answer, +1 – LLlAMnYP Oct 30 '17 at 10:38

This question has popped up at stats.stackexchange.com ... where alternative algorithms (including solved in mma) are discussed. See:

http://stats.stackexchange.com/questions/193990/approximate-e-using-monte-carlo-simulation

The code I suggested there is:

Mean[Table[Module[{u=RandomReal[], t=1},  While[u<1, u=RandomReal[]+u; t++]; t] , {10^6}]]


{0.208377, 2.71887}

... which is about 50 times faster than the NestWhile approach:

Mean[ Length /@
Table[NestWhile[Append[#, RandomReal[]] &, {RandomReal[]},
Total@# < 1 &], 10^6]] // N // AbsoluteTiming


{9.45841, 2.71663}

And:

counter[cnt_, total_] /; total > 1. := cnt
counter[cnt_, total_] := counter[cnt + 1, total + RandomReal[]]
Mean[Table[counter[0, 0], 1000000]] // N // AbsoluteTiming


{5.30918, 2.71847}

Notes

1. For an explanation of why the relation holds, see:

http://stats.stackexchange.com/questions/194352/

2. For larger samples, the above can be further improved replacing Table with ParallelTable.

• I would advocate against using Random since the underlying method is not necessarily good as an RNG. – Daniel Lichtblau Oct 30 '17 at 15:37
• Thanks Daniel - updated to RandomReal – wolfies Oct 30 '17 at 15:43