An unexpected result involving sums of random numbers

[Gerson W. Barbosa]

How many random numbers in the interval [0, 1) are required to add up to 1 or higher, in average?

The following HP-42S and HP-15C programs might give an answer to this question in under one minute on the physical calculators (with some cheating due to the choice of seeds and number of results equal or greater than one), or in reasonable on Free42 and other simulators with a large number of results, regardless of the seeds.

HP-42S:

Code:


00 { 36-Byte Prgm }
01▸LBL "X"
02 STO 00
03 0
04 STO 01
05 SIGN
06▸LBL 00
07 RAN
08▸LBL 01
09 X<>Y
10 STO+ 01
11 X<>Y
12 RAN
13 +
14 X<Y?
15 GTO 01
16 R↓
17 DSE ST Y
18 GTO 00
19 RCL 01
20 RCL+ 00
21 RCL÷ 00
22 END

6 SEED 252 XEQ “X” -> ______________ (55.5 seconds on the real calculator)

Code:


00 { 40-Byte Prgm }
01▸LBL "X"
02 STO 00
03 1
04 SIGN
05▸LBL 00
06 RAN
07▸LBL 01
08 RCL- ST L
09 RAN
10 -
11 1
12 STO+ ST Z
13 +
14 X>0?
15 GTO 01
16 R↓
17 DSE ST Y
18 GTO 00
19 DSE ST X
20 RCL+ 00
21 RCL÷ 00
22 END

67 SEED 252 XEQ “X” -> ______________ (63 seconds on the real calculator)

HP-15C:

Code:


   001-     42 21 11   f LBL A
   002-        44 25     STO I
   003-        44  0     STO 0
   004-            0     0
   005-        44  1     STO 1
   006-        42  0   f x!        ; 1
   007-     42 21  0   f LBL 0
   008-        42 36   f RAN #
   009-     42 21  1   f LBL 1
   010-           34     x<>y
   011-     44 40  1     STO+ 1
   012-           34     x<>y
   013-        42 36   f RAN#
   014-           40     +
   015-     43 30  8   g TEST 8    ; x<y
   016-        22  1     GTO 1
   017-           33     Rv
   018-     42  5 25   f DSE I
   019-        22  0     GTO 0
   020-        45  1     RCL 1
   021-     45 40  0     RCL+ 0
   022-     45 10  0     RCL/ 0
   023-        43 32   g RTN

6 STO RAN# 32 f A -> ______________ (55 seconds on the real HP-15C)

[Valentin Albillo]

.
Hi, Gerson,

(02-28-2023 06:41 PM)Gerson W. Barbosa Wrote:  How many random numbers in the interval [0, 1) are required to add up to 1 or higher, in average?

Excuse me, Gerson, not to rain on your parade or anything but isn't this the same thing I recently discussed in my Concoction the First: Weird limit sub-challenge ?

Regards.
V.

[Gerson W. Barbosa]

(02-28-2023 07:50 PM)Valentin Albillo Wrote:  .
Hi, Gerson,

(02-28-2023 06:41 PM)Gerson W. Barbosa Wrote:  How many random numbers in the interval [0, 1) are required to add up to 1 or higher, in average?

Excuse me, Gerson, not to rain on your parade or anything but isn't this the same thing I recently discussed in my Concoction the First: Weird limit sub-challenge ?

Hello, Valentin,

Indeed it is, albeit this is only the first part of #1. I did participate in that one, but I had completely forgotten about it. I only tried to address #6 and overlooked all the rest, it seems.
I read about this subject one of these days, in a short note in Prime Obsession, a book by John Derbyshire and decided to check it out. I checked also for 2, e and pi, but except for the first one I was unable to identify the constants.
Anyway, we can consider this a (very) belated solution to that particular nice sub-challenge of yours :-)

Best regards,

Gerson.

[Namir]

The answer is a squewed distribution (Gauss, or Poisson, or something else?).

I did some quick calculations in Matlab. For the sum of 1, the most frequent number of random numbers is 2. The distribution dies out fast. For the sum of 2, the most frequent number of random numbers is 4. Also, the distribution dies out fast.

Namir

[Gerson W. Barbosa]

(02-28-2023 11:08 PM)Namir Wrote:  The answer is a squewed distribution (Gauss, or Poisson, or something else?).

Hi, Namir,

Perhaps my phrasing was not so clear. Let’s do an example by hand on the HP-15C, step by step:

0 STO RAN#; start the sequence of random numbers with zero

1)
f RAN# -> 0.1018
f RAN# -> 0.7365
+ -> 0.8383 ; sum < 1, go on:
f RAN# -> 0.3248
+ -> 1.1631 ; sum >= 1 obtained with 3 random numbers.

2)
f RAN# -> 0.3485
f RAN# -> 0.8685
+ -> 1.2171 ; sum >= 1 obtained with 2 random numbers.

3)
f RAN# -> 0.2789
f RAN# -> 0.3615
+ -> 0.6404 ; sum < 1, go on:
f RAN# -> 0.1074
+ -> 0.7448 ; sum still < 1, go on:
f RAN# -> 0.2629
+ -> 1.0106 ; sum >= 1 obtained with 4 random numbers.

4)
f RAN# -> 0.9252
f RAN# -> 0.5599
+ -> 1.4851 ; sum >= 1 obtained with 2 random numbers.

Average: (3 + 2 + 4 + 2)/4 = 2.75

That’s exactly what we get when running the program:

0 STO RAN# 4 f A -> 2.7500

The result, 2.75, is close to the constant we’re looking for.

I thought of providing an HP-41C program, but then I remembered it lacks a random number generator.

Gerson.

[KeithB]

The answer I got was about 3.5. However that was with a normal distribution with mean of 0.5 and standard deviation of 1. 8^)

While I did it in Matlab, it would be easy on the Prime.

[Gil]

Could anybody present a nice, understandable proof for one of the results, for example for the average "—> e" in case "sum of necessary numbers >=1"?

In a far-fetched way, it reminds me of Benford law, with % of number i(i=1,2... 9) at first place = log (1+1/number i).

Regards,
Gil

[robve]

This post reminded me of something I had seen before posted on Reddit r/programming a couple of years back.

The mean-sum-rand-to-one turns out to be a kinda-sorta-known result, if you're "obsessed with primes"

Spoiler alert, do not follow these links: a wonderful article at Mostlymaths with an explanation at Math.Stackexchange

Quite intriguing.

- Rob

[Gerson W. Barbosa]

(03-01-2023 03:36 PM)robve Wrote:  Spoiler alert, do not follow these links: a wonderful article at Mostlymaths with an explanation at Math.Stackexchange

Thanks for the links! These alone have made worthwhile starting this thread.

Gerson.

[Gerson W. Barbosa]

(03-01-2023 03:54 AM)Gerson W. Barbosa Wrote:  Average: (3 + 2 + 4 + 2)/4 = 2.75

Actually the programs compute the average as

[(2 + 1 + 3 + 1) + 4]/4 = 2.75,

thus doing only six additions instead of ten.

Edited to fix a mistake

[Namir]

(03-01-2023 03:54 AM)Gerson W. Barbosa Wrote:  

(02-28-2023 11:08 PM)Namir Wrote:  The answer is a squewed distribution (Gauss, or Poisson, or something else?).

Hi, Namir,

Perhaps my phrasing was not so clear. Let’s do an example by hand on the HP-15C, step by step:

0 STO RAN#; start the sequence of random numbers with zero

1)
f RAN# -> 0.1018
f RAN# -> 0.7365
+ -> 0.8383 ; sum < 1, go on:
f RAN# -> 0.3248
+ -> 1.1631 ; sum >= 1 obtained with 3 random numbers.

2)
f RAN# -> 0.3485
f RAN# -> 0.8685
+ -> 1.2171 ; sum >= 1 obtained with 2 random numbers.

3)
f RAN# -> 0.2789
f RAN# -> 0.3615
+ -> 0.6404 ; sum < 1, go on:
f RAN# -> 0.1074
+ -> 0.7448 ; sum still < 1, go on:
f RAN# -> 0.2629
+ -> 1.0106 ; sum >= 1 obtained with 4 random numbers.

4)
f RAN# -> 0.9252
f RAN# -> 0.5599
+ -> 1.4851 ; sum >= 1 obtained with 2 random numbers.

Average: (3 + 2 + 4 + 2)/4 = 2.75

That’s exactly what we get when running the program:

0 STO RAN# 4 f A -> 2.7500

The result, 2.75, is close to the constant we’re looking for.

I thought of providing an HP-41C program, but then I remembered it lacks a random number generator.

Gerson.

Your results are good. With one million iteration I get an average of 2.71. With a sum of 2 I get an average of 4.6. And with an average of 3 I get an average of 6.5.

Namir

[Albert Chan]

If random values are discrete, say coin flips, (head=1, tail=0), the expected value (sum ≥ 1) = 1/p = 2

Code:

        n * P(n)
1       1 * 1/2
01      2 * 1/4
001     3 * 1/8
0001    4 * 1/16
...

This gives a neat prove for E = Σ(n/2^n, n=1 .. Inf) = 1/p = 2
see https://math.stackexchange.com/questions...onverge-to

At the end of each simulation, sum is exactly 1, without "waste".

---

https://www.fourmilab.ch/documents/random_sum/ examples ("waste" in bold)

0.4770 + 0.0516 + 0.1793 + 0.2250 + 0.0438 + 0.3006 = 1.2773
0.8474 + 0.3005 = 1.1479
0.5535 + 0.2721 + 0.5390 = 1.3646
...

For continuous uniform distribution of [0, 1), we do have wastage (almost always)

Last random number, on average = 0.5
Balance (before adding last random number), on average, between 0.5 to 1.0

1.0 < E*0.5 < 1.5      → 2.0 < E < 3.0

FYI, above link has the prove for E = e ≈ 2.71828

[Gil]

On HP50G, after a loop of 100 000,
I get an average of 2.80983 for "sum > 1"
and 2.58946 for "sum >= 1"
with seed (RDZ) =1.

It does not seem to converge to e, or very, very slowly.

Is it normal?

[Gerson W. Barbosa]

(03-01-2023 07:14 PM)Namir Wrote:  Your results are good. With one million iteration I get an average of 2.71. With a sum of 2 I get an average of 4.6. And with an average of 3 I get an average of 6.5.

I have modified two programs to handle arbitrary sums. Even 5000 iterations, which is quite feasible on the HP-15C LE, will give results close to the theoretical values:

1 -> 2.7182818; e

2 -> 4.6707743; e(e - 1)

3 -> 6.6665656; e(2e² - 4e + 1)/2

Code:


   001-    42 21 11   f LBL A
   002-       44 25     STO I
   003-          34     x<>y
   004-          20     ×
   005-       44  0     STO 0
   006-           0     0
   007-       44  1     STO 1
   008-       43 36   g LSTΧ
   009-    42 21  0   f LBL 0
   010-       42 36   f RAN#
   011-    42 21  1   f LBL 1
   012-          34     x<>y
   013-    44 40  1     STO+ 1
   014-          34     x<>y
   015-       42 36   f RAN#
   016-          40     +
   017-    43 30  8   g TEST 8    ; x<y
   018-       22  1     GTO 1
   019-          33     Rv
   020-    42  5 25   f DSE I
   021-       22  0     GTO 0
   022-       45  1     RCL 1
   023-    45 40  0     RCL+ 0
   024-    45 10  0     RCL/ 0
   025-       43 32   g RTN

0 STO RAN# 1 ENTER 5000 f A -> 2.7152 ( 55.5 seconds )

0 STO RAN# 2 ENTER 5000 f A -> 4.6656 ( 100.6 seconds )

0 STO RAN# 3 ENTER 5000 f A -> 6.6754 ( 146.7 seconds )


Of course, on a computer or smartphone running Free42 we can go even further. For instance,

Code:

00 { 40-Byte Prgm }
01▸LBL "X"
02 X<>Y
03 RCL× ST Y
04 STO 00
05 CLX
06 STO 01
07 X<> ST L
08▸LBL 00
09 RAN
10▸LBL 01
11 X<>Y
12 STO+ 01
13 X<>Y
14 RAN
15 +
16 X<Y?
17 GTO 01
18 R↓
19 DSE ST Y
20 GTO 00
21 RCL 01
22 RCL+ 00
23 RCL÷ 00
24 END

4 SEED 1 ENTER 1E7 XEQ “X” -> 2.7181530 ( 29.4 seconds )

4 SEED 2 ENTER 1E7 XEQ “X” -> 4.6711803 ( 50.1 seconds )

4 SEED 3 ENTER 1E7 XEQ “X” -> 6.6667329 ( 71.6 seconds )


As mentioned earlier in this thread, I was motivated by a short note on a book I read some years ago and ran into it again last week, which read

Quote:Here is an example of e turning up unexpectedly. Select a random number between 0 and 1. Now select another and add it to the first. Keep doing this, piling on random numbers. How many random numbers, on average, do you need to make the total greater than 1? The answer is e.

Not forgetting about Valentin's Math Challenge linked above would have spared me some time, but then again I would have lost the chance to dust off my calculators and check that out for myself. Detailed information and a shorter and possibly faster RPN program also available in his original solution.

Thanks you all for your interest and contributions.

Gerson.

[KeithB]

In Matlab with the following script:
sumcount = 0;
for i = 1:100000
s = 0;
count = 0;
while s <= 1
s = s+random('Uniform', 0,1);
count = count + 1;
end
sumcount = sumcount + count;
end

avgcount = sumcount/100000;

I get:
2.7195
2.7196
2.7190
2.7175
2.7183

Or, much better than yours. Of course this does not really "converge" to e.

[Gil]

On HP50G, after a loop of 100 000,
I get an average of 2.80983 for "sum > 1"
and 2.58946 for "sum >= 1"
with seed (RDZ) =1.

And, after a loop of 1000 000,
I get an average of 2.80983 for "sum > 1"
and 2.590195 for "sum >= 1" for "sum >= 1"
with seed (RDZ) =1.

It does not seem to converge to e, or very, very slowly.

The program is the following:

\<< TIME 0. 'N' STO 1. RDZ 1. 1000000.
START 0. 'S' STO
WHILE S 1. \<=
REPEAT RAND 1. RND 'S' STO+ 1. 'N' STO+
END
NEXT N 1000000. / TIME ROT HMS-
\>>
Is it normal?

[Gil]

Sorry.

I added a "1 RDZ" unduly in the second part of the program.

The found values indeed —> clearly to 2.7 or 2.8, ie e.

"\\<< TIME 0. 'N' STO 1. RDZ 1. 100000.
START 0. 'S' STO
WHILE S 1. \\<=
REPEAT RAND 'S' STO+ 1. 'N' STO+
END
NEXT N 100000. / TIME ROT HMS-
\\>>"

Gives 2.71959.

With my apologies.

[Gerson W. Barbosa]

Same result on the HP 50g. It should take longer on the HP-48G/GX. Even longer on the HP-28S.

1 RDZ 100000

« 0 DUP2 NOT SWAP
START RAND
DO RAND + SWAP 1 + SWAP DUP
UNTIL 1 >
END DROP
NEXT SWAP / 1 +
»

TEVAL ->

2.71959
s:1365.9755 (HP 50g)

s:57.9261 (iHP48)

[Albert Chan]

(03-02-2023 02:06 PM)Albert Chan Wrote:  Last random number, on average = 0.5
Balance (before adding last random number), on average, between 0.5 to 1.0

1.0 < E*0.5 < 1.5      → 2.0 < E < 3.0

I was wrong about last random number average = 0.5
Balance likely contains many small numbers. Last random thus pushed higher (on average).

lua> N = 1e7 -- simulations
lua> S, R, E = 0, 0, 0
lua> for i = 1, N do
:            s = 0; repeat r=random(); s=s+r; E=E+1 until s >= 1
:            S, R = S+s, R+r
:      end

lua> S, R, E = S/N, R/N, E/N
lua> B = S-R
lua> B, R, S
0.7182663624489083      0.6408642227941028      1.3591305852430111

lua> e = exp(1) -- simulation suggested (B+R = S) values, for huge n
lua> e-2, 2-e/2, e/2
0.7182818284590451      0.6408590857704775      1.3591409142295225

We expected E*0.5 ≈ S. Here, expected counts estimated with 2S is better than E.

lua> E, 2*S
2.7186766      2.7182611704860222

[Gil]

And Barbosa's program on EMU48, with Samsung A53:
Less than 13s in approximate mode [/php]for average "1 to 100 000 and seed =1".

Parabéns, Gérson!