Barbara writes a sequence of integers starting with the number 12. Each subsequent
integer she writes is chosen randomly with equal probability from amongst the positive divisors of the previous integer (including the possibility of the integer itself). She keeps writing integers until she writes the integer 1 for the first time, and then she stops.
An example of one such sequence is 12, 6, 6, 3, 3, 3, 1.
What is the expected value of the number of terms in Barbara’s sequence?
If n is the starting number (where n=12 for the given problem), then the expected number of terms in the sequence, x, can be seen to be:
x = 1 + (c/(c+1))*e + (1/(c+1))*x
where e is the average expectation for the other possible chosen divisors of n, and c is the number of other such divisors, including 1, so that c+1 is the total number of divisors including n itself. The reason for this is:
The 1 accounts for this number in the sequence itself.
The (c/(c+1))*e represents the c/(c+1) probability that a divisor other than n itself will be chosen.
The (1/(c+1))*x represents the 1/(c+1) probability that n itself will be chosen as the divisor, with of course the same expectation again.
Solving for x, this comes out to
x = 1 + e + 1/c
This allows the recursive calculation of the expectation:
DECLARE FUNCTION expct# (n#)
DEFDBL A-Z
PRINT expct(12)
END
FUNCTION expct (n)
IF n = 1 THEN expct = 1: EXIT FUNCTION
tot = 0: ct = 0
FOR i = 1 TO n - 1
IF n MOD i = 0 THEN
tot = tot + expct(i)
ct = ct + 1
END IF
NEXT
e = tot / ct
PRINT n; 1 + e + 1 / ct
expct = 1 + e + 1 / ct
END FUNCTION
It finds the expectation for n=12 to be 4.033333333333333, or 4 + 1/30, or 121/30.
Along the way it reports the expectations for the other divisors of 12:
2: 3
3: 3
4: 3.5
6: 3 + 2/3 = 11/3
with of course the known expected value for 1, as 1.
As a check, 121/30 is indeed 1 more than the average of 1, 3, 3, 3.5, 11/3 and 121/30.
A simulation, likewise checks:
DO
n = 12: ct = 1
DO
PRINT n;
ct = ct + 1
num = 0
FOR i = 1 TO n
IF n MOD i = 0 THEN
num = num + 1
d(num) = i
END IF
NEXT
r = INT(RND(1) * num + 1)
n = d(r)
LOOP UNTIL n = 1
PRINT 1, ct
totCt = totCt + ct: trials = trials + 1
PRINT totCt; trials, totCt / trials
LOOP
which finds, at 55,325 trials, 222,991 total count of terms used in the sequences, for an average of 4.030565.
For n other than 12, here is a table of expected sequence lengths:
2 3
3 3
4 3.5
5 3
6 3.666666666666667
7 3
8 3.833333333333333
9 3.5
10 3.666666666666667
11 3
12 4.033333333333333
13 3
14 3.666666666666667
15 3.666666666666667
16 4.083333333333334
17 3
18 4.033333333333333
19 3
20 4.033333333333333
For a prime, it's 3; for the square of a prime it's 3.5; for a product of two primes, it's 3+2/3; for the square of a prime multiplied by another prime, it's the current answer, etc.
|
|
Posted by Charlie
on 2010-04-09 16:27:06 |
solution
