A366995 - OEIS

%I #14 Nov 01 2023 10:01:13

%S 3,39,359,2477,119667,1522705,46419629,6143100517,5472109127035,

%T 790136773603303,278129286200597661,16684426086791338103,

%U 503067648850136040148699,2626565018569118643191009,10920130209346850287269887104735,236686188450953790757840351941895

%N a(n) is the numerator of the expected end time of a game with three gamblers, one of which starts with capital n, the others with capital 1 each, conditional on the event that one of the two poor players wins.

%C In each round of the game, 1 unit is transferred from one randomly chosen player to another. Players play until they are out of money, so when the first player is out the other two continue to play. The winner is the player who ends up with all n+2 units of money.

%H Pontus von Brömssen, <a href="/A366995/b366995.txt">Table of n, a(n) for n = 1..53</a>

%H Persi Diaconis and Stewart N. Ethier, <a href="https://doi.org/10.1214/21-STS826">Gambler’s ruin and the ICM</a>, Statist. Sci. 37 (3) 289-305, August 2022.

%H Persi Diaconis, <a href="https://sites.math.rutgers.edu/~zeilberg/expmath/diaconis23.pdf">Gambler's ruin with k gamblers</a>, slides from talk in the Rutgers Experimental Mathematics Seminar, <a href="https://sites.math.rutgers.edu/~zeilberg/expmath/archive23.html">Fall 2023 Semester</a>, Oct 12, 2023.

%H Experimental Mathematics, <a href="https://vimeo.com/876882233">Gambler’s ruin with k gamblers</a>, recording of talk, Vimeo video, Oct 22, 2023.

%o (Sage)

%o from itertools import permutations

%o def T(n):

%o     nodes = [(i,j) for i in range(n+2) for j in range((n+2-i)//2+1)]

%o     m = len(nodes)

%o     Q0 = {x:{y:0 for y in nodes} for x in nodes}

%o     for x in nodes:

%o         c1 = x+(n+2-sum(x),)

%o         for i,j in permutations(range(3),int(2)):

%o             if c1[i] and c1[j]:

%o                 c2 = list(c1)

%o                 c2[i] -= 1

%o                 c2[j] += 1

%o                 y = (c2[0],min(c2[1:]))

%o                 if c2[0] != n+2:

%o                     Q0[x][y] += n+2-c2[0]

%o     Q0 = matrix(QQ,[list(R.values()) for R in Q0.values()])

%o     s = sum(Q0.columns())

%o     Q = identity_matrix(QQ,m-1)

%o     for i in range(1,m):

%o         for j in range(1,m):

%o             if s[i] != 0: Q[i-1,j-1] -= Q0[i,j]/s[i]

%o     return (Q**(-1)*ones_matrix(QQ,m-1))[-2,0]

%o def A366995(n):

%o     return T(n).numerator()

%o def A366996(n):

%o     return T(n).denominator()

%Y Cf. A366566 (a(n)/A366996(n) rounded to nearest integer), A366996 (denominators).

%K nonn,frac

%O 1,1

%A _Pontus von Brömssen_, Oct 31 2023