From: israel@math.ubc.ca (Robert Israel) Subject: Re: adding random digits until digit 0 Date: 13 Apr 1999 17:32:46 GMT Newsgroups: sci.math Keywords: example of probability generating function In article <7etnjm$ah5@sjx-ixn1.ix.netcom.com>, " r e s" writes: |> This looks like a very simple problem, but I can |> find only a partial solution: |> |> Independent uniformly distributed decimal digits |> are generated in sequence until a 0 is obtained. |> Let S be the sum of the digits. |> What is pr(S=k) for k=0,1,2,...? Let f(z) = sum_{k=0}^infinity pr(S=k) z^k be the generating function. You can think of this as the expected payoff in a game. You start off with S = 0. At each round you generate a digit. If the digit is not 0, add it to S. If the digit is 0, the game ends with payoff z^S. It's easy to see by a "first-step analysis" that f(z) = 1/10 (1 + sum_{j=1}^9 z^j f(z)) so that f(z) = 1/(10 - sum_{j=1}^9 z^j). You can extract pr(S=k) from the Taylor series of this. If you want a more explicit formula, expand f(z) in partial fractions (involving the roots of the polynomial 10 - sum_{j=1}^9 z^j, which is irreducible over the rationals). According to Maple: > convert(f,fullparfrac,z); ----- \ / 93206534790699 569851039164221 ) |- ----------------- - ----------------- _alpha / \ 48689722320049889 48689722320049889 ----- _alpha = %1 18924701993010 7 9414801494010 8 - ----------------- _alpha - ----------------- _alpha 48689722320049889 48689722320049889 28530662093010 6 38233652093010 5 - ----------------- _alpha - ----------------- _alpha 48689722320049889 48689722320049889 48034652093010 4 57934652093010 3 - ----------------- _alpha - ----------------- _alpha 48689722320049889 48689722320049889 67934652093010 2\ - ----------------- _alpha |/(z - _alpha) 48689722320049889 / %1 := RootOf( 2 3 4 5 6 7 8 9 -10 + _Z + _Z + _Z + _Z + _Z + _Z + _Z + _Z + _Z ) and the z^k term is thus ----- \ / 93206534790699 569851039164221 ) | ----------------- + ----------------- _alpha / \ 48689722320049889 48689722320049889 ----- _alpha = %1 18924701993010 7 9414801494010 8 + ----------------- _alpha + ----------------- _alpha 48689722320049889 48689722320049889 28530662093010 6 38233652093010 5 + ----------------- _alpha + ----------------- _alpha 48689722320049889 48689722320049889 48034652093010 4 57934652093010 3 + ----------------- _alpha + ----------------- _alpha 48689722320049889 48689722320049889 67934652093010 2\ k / (k+1) + ----------------- _alpha | z / _alpha 48689722320049889 / / %1 := RootOf( 2 3 4 5 6 7 8 9 -10 + _Z + _Z + _Z + _Z + _Z + _Z + _Z + _Z + _Z ) Robert Israel israel@math.ubc.ca Department of Mathematics http://www.math.ubc.ca/~israel University of British Columbia Vancouver, BC, Canada V6T 1Z2 ============================================================================== From: israel@math.ubc.ca (Robert Israel) Subject: Re: adding random digits until digit 0 Date: 14 Apr 1999 21:15:46 GMT Newsgroups: sci.math In article <7f13q8$5g3@dfw-ixnews4.ix.netcom.com>, " r e s" writes: |> Robert Israel wrote ... |> > and the z^k term is thus |> > |> > |> > ----- |> > \ / 93206534790699 569851039164221 |> > ) | ----------------- + ----------------- _alpha |> > / \ 48689722320049889 48689722320049889 |> > ----- |> > _alpha = %1 |> > |> > 18924701993010 7 9414801494010 8 |> > + ----------------- _alpha + ----------------- _alpha |> > 48689722320049889 48689722320049889 |> > |> > 28530662093010 6 38233652093010 5 |> > + ----------------- _alpha + ----------------- _alpha |> > 48689722320049889 48689722320049889 |> > |> > 48034652093010 4 57934652093010 3 |> > + ----------------- _alpha + ----------------- _alpha |> > 48689722320049889 48689722320049889 |> > |> > 67934652093010 2\ k / (k+1) |> > + ----------------- _alpha | z / _alpha |> > 48689722320049889 / / |> > |> > %1 := RootOf( |> > |> > 2 3 4 5 6 7 8 9 |> > -10 + _Z + _Z + _Z + _Z + _Z + _Z + _Z + _Z + _Z ) |> |> Is this to be read as %1 ranging over the roots of this polynomial, |> so that the above summation is over these (9?) roots? (So that the |> factor z^k can be taken outside the summation, making the summation |> part then equal to pr(S=k)?) Yes. |> Is it reasonable to ask for Maple to list these roots, just for |> completeness? (But without something like Maple, it does appear |> that numerical precision will make the formula of limited use |> for computation.) The roots are presumably not expressible in "closed form", but numerical approximations to any desired accuracy are no problem: > allvalues(RootOf(-10+_Z+_Z^2+_Z^3+_Z^4+_Z^5+_Z^6+_Z^7+_Z^8+_Z^9)); -1.302445483-.4524439746*I, -1.302445483+.4524439746*I, -.7320397894-1.149065384*I, -.7320397894+1.149065384*I, .1378726129-1.318816774*I, .1378726129+1.318816774*I, .8861118114-.8906948008*I, .8861118114+.8906948008*I, 1.021001695 If you call these roots r[i], then pr(S=k) = sum(c[i]/r[i]^(k+1), i=1..9) where the c[i] are the following: -.01309690761 - .004561743373 I, -.01309690761 + .004561743373 I, -.007333534801 - .01157953096 I, -.007333534801 + .01157953096 I, .001458726828 - .01326701544 I, .001458726828 + .01326701544 I, .009036853745 - .008883408210 I, .009036853745 + .008883408210 I, .01986972364 The real root 1.021001695 is the smallest in absolute value, so as k -> infinity it dominates the asymptotics: pr(S = k) ~= .01986972364/1.021001695^(k+1) Robert Israel israel@math.ubc.ca Department of Mathematics http://www.math.ubc.ca/~israel University of British Columbia Vancouver, BC, Canada V6T 1Z2 ============================================================================== From: N0$PAM.Bi-Zed@f128.n5015.z2.fidonet.org (Max Alekseyev) Subject: adding random digits until digit 0 Date: 13 Apr 99 22:11:48 +0300 Newsgroups: sci.math Hi, ! Replying to a message of r e s to All: res> Independent uniformly distributed decimal digits res> are generated in sequence until a 0 is obtained. res> Let S be the sum of the digits. res> What is pr(S=k) for k=0,1,2,...? res> For k<10, pr(S=k) has a simple form, but the res> pattern obviously doesn't generalize: res> pr(S=k)=p, k=0 res> pr(S=k)=(p^2)*(1+p)^(k-1), k=1,2,...,9 res> where p=1/10. res> Any suggestions? Sum Pr(S=k)*x^k = (1-x)/(x^10-11x+10) = k = 1/10 + 1/100 x + 11/1000 x^2 + 121/10000 x^3 + 1331/100000 x^4 + 14641/1000000 x^5 + ... We have Pr(S=0)=1/10; Pr(S=1)=1/100; Pr(S=2)=11/1000 and so on... For example, Pr(S=50)=6883985816993102806562986072909868756342081702091/ 10^51 Regards, ø ø Max ~