How do I compute hitting time in a finite Markov chain? (to quantify my Machi Koro pwnage)

Question

In a game of Machi Koro, I had the Sushi Bar (active) while my opponent had a Wheat Field and Bakery with no low-roll buildings available. I want to know how long they have to roll before they can afford the Train Station.

That is, they're engaged in a Markov chain where the state space is "0 to 4 coins" ($\{0, 1, 2, 3, 4\}$) and I want to know the probability, for each $n$, that they reach $4$ for the first time after exactly $n$ transitions (preferably for all starting states, but just for zero is interesting enough). I might also be interested in the cumulative version: for each $n$, what's the likelihood that they'll have reached $4$ at least once after $n$ transitions.

The transition probabilities are as follows: from $k$ coins, $k < 4$: to $1$ with probability $1/6$; otherwise, to $k$ with probability $3/6$ and to $k+1$ with probability $2/6$. Since the transitions from $4$ to anything are not interesting, I'll just define them arbitrarily as $4$ to $4$ with probability $1$ (really they'll spend all their coins, but the questions don't care what happens after $4$ is reached).

Spelled out:

$ \left[ \begin{matrix} 3/6 & 3/6 & & & \\ & 4/6 & 2/6 & & \\ & 1/6 & 3/6 & 2/6 & \\ & 1/6 & & 3/6 & 2/6 \\ & & & & 1 \end{matrix} \right] $

(The probability in row $i$ column $j$ is the probability of transition from $i$ to $j$ coins in a single step.)

Some analysis: in the transition graph, there are strongly connected components $\{0\}, \{1, 2, 3\}$ and $\{4\}$ with transitions $0 \rightarrow 1$ and $3 \rightarrow 4$. Maybe it helps to analyse each scc separately. You go from zero coins to one in $k$ steps, $k > 0$, with probability $2^{-k}$; the big one essentially boils down to "how soon will you increment your state three times in a row (each w. prob. $2/6$) without falling back to $1$ (w. p. $1/6$ each step)". I'm not sure how to analyze this, though.

I could probably spell my way through https://en.wikipedia.org/wiki/Markov_chain but given that my matrix-fu is a little rusty, some hand-holding would be appreciated.

Edit: Come to think of it, given the City Hall the transition from $0$ is to $1$ with probability $1$. The 50/50-out-of-0 would occur if I had three Cafés instead of the Sushi Bar and we were playing without the Harbor expansion. Feel free to analyze either—or both! Also, for purposes of simplifying the math, in both cases I'm ignoring the fact that my opponent earns a coin whenever I roll a 1.

Answer 1

A bit of simulation can say something about the shape of the distribution. See below for the full simulation code.

Based on 100,000 games, the quartile (lower) boundaries are at 8, 12 and 20; that is, in fewer than 25% of games four coins are reached by turn 8 or earlier, and in more than 25% of games four coins are reached by turn 9 (or earlier).

The decile (lower) boundaries are 5, 7, 8, 10, 12, 14, 18, 22, 29.

Based on a slightly simplistic simulation of two competing strategies which focus on the Convenience Store and Cheese Factory, respectively, the median and mode game length are between 25 and 30 (in 1000 games), ranging from 20 to 60, falling off quickly from 40 and beyond.

The cumulative likelihood of reaching 4 for the first time by turn $n$, starting with 0 coins, is given below, up to 50%. Note that it can't be done before turn 4. The right tail goes towards 100% as $n$ goes to infinity, of course; it crosses 99% in turns 54-55, 99.9% in turns 79-80 and 99.99 in turns 108-109.

$ \begin{array}{r|r} \textrm{turn} & \textrm{cumul. prob.} \\ \hline 4 & 1.825\% \\ 5 & 5.787\% \\ 6 & 11.450\% \\ 7 & 17.811\% \\ 8 & 24.228\% \\ 9 & 30.599\% \\ 10 & 36.591\% \\ 11 & 42.104\% \\ 12 & 47.358\% \\ 13 & 52.032\% \\ \end{array} $

Here's the simulation code:

#!/somewhere/over/the/rainbow/python
from random import choice
from collections import Counter

probs = [(0, 1),
         (1, 1, 2),
         (1, 2, 2, 2, 3, 3),
         (1, 3, 3, 3, 4, 4)]

def simulation(state=0, target=4):
    k = 0
    while state != target:
        state = choice(probs[state])
        k += 1
    return k

def many(n, state=0, target=4):
    return Counter(simulation(state, target) for _ in range(n))

n = 100000
histogram = many(n)
cumulative = 0.0
pct = lambda x: 100 * x / float(n)
for i in range(max(histogram.keys())):
    cumulative += histogram[i]
    print '%5d: % 7.3f%% (cumulative % 7.3f%%)' % \
        (i, pct(histogram[i]), pct(cumulative))

Have fun running your own simulations. An analytical answer, and an explanation of the techniques for finding it, will still be much appreciated :-)

Answer 2

I found the following answer to a related question: https://math.stackexchange.com/a/1317088/234337.

I also derived the main formula myself. Let $t_i$ be the expected number of transitions before 4 coins is reached, when starting with $i$ coins. Then $t_0 = 1 + (\frac12 t_0 + \frac12 t_1)$, implying that $t_0 = 2 + t_1$. By similar substitutions and a bit of fiddly algebra, $t_1 = 3 + t_2$ and $t_2 = \frac92 + t_3$ and $t_3 = \frac{81}{12}$. Putting this together, $t_0 = \frac{65}4 = 16.25$.

I also read a bit of the book linked in the answer, which helped me understand what linear algebra to do. Here's some code to compute an exact answer to the question "what's the probability of having hit 4 coins at least once in turn $n$ or earlier":

import numpy
from fractions import Fraction as frac

p = [[frac(1, 2), frac(1, 2),          0,          0,          0],
     [         0, frac(4, 6), frac(2, 6),          0,          0],
     [         0, frac(1, 6), frac(3, 6), frac(2, 6),          0],
     [         0, frac(1, 6),          0, frac(3, 6), frac(2, 6)],
     [         0,          0,          0,          0,          1]]

p = numpy.matrix(p)
cumulative = p.copy() * p * p

pct = lambda frac: float(100 * frac)
for i in range(4, 60):
    before = cumulative[0, 4]
    cumulative *= p
    after = cumulative[0, 4]
    print '%2d %7.3f (%5.3f)' % (i, pct(after), pct(after - before))

The distribution it prints approximately matches the result of my simulations in another answer; the deviations are small enough that I don't suspect a wrong implementation or a correct implementation of the wrong thing.

I'm stuck with one observation: the second quartile cutoff is between 12 and 13, yet the expected hitting time is $\not \in [12, 13]$, which I would naively expect. I suspect that this is because the cutoff is the median hitting time whereas the expectation is the mean; in the land of the blind, everyone but the one-eyed king has below average vision.

The following code may be helpful if you want to edit the above to take the other player's Wheat Field into account; each state is then A's coins at the start of B's turn, where B (me) has the three Cafés.

from collections import Counter
for k in range(4):
    histogram = Counter()
    for me in range(1, 7):
        for you in range(1, 7):
            coins = k + (me == 1)
            if you == 1:
                coins = max(coins - 3, 0) + 1
            elif you in (2, 3):
                coins += 1
            histogram[coins] += 1
    print k, sorted(histogram.items())