Riddler: Nonconformist Dice

Last week’s Riddler Classic is a question about rolling tetrahedral dice:

You have four fair tetrahedral dice whose four sides are numbered 1 through 4.
You play a game in which you roll them all and divide them into two groups: those whose values are unique, and those which are duplicates. For example, if you roll a 1, 2, 2 and 4, then the 1 and 4 will go into the “unique” group, while the 2s will go into the “duplicate” group.
Next, you reroll all the dice in the duplicate pool and sort all the dice again. Continuing the previous example, that would mean you reroll the 2s. If the result happens to be 1 and 3, then the “unique” group will now consist of 3 and 4, while the “duplicate” group will have two 1s.
You continue rerolling the duplicate pool and sorting all the dice until all the dice are members of the same group. If all four dice are in the “unique” group, you win. If all four are in the “duplicate” group, you lose. What is your probability of winning the game?

We will answer this question using a Markov Chain in python. For that, we will use numpy (for matrix multiplication), collections.Counter (for easy counting) and itertools.product (for throwing dice).

from collections import Counter
from itertools import product
import numpy as np

First, let’s calculate the probability of every possible throw with four tetrahedral dice. We will represent a dice throw as an ordered string.

# all possible throws
throws = ["".join(str(i) for i in state) for state in product([1,2,3,4], repeat=4)]
# aggregate identical throws and calculate probabilities
all_states = ["".join(sorted(throw)) for throw in throws]
freqs = Counter(all_states).items()
throw = {k: v/(4**4) for k, v in freqs}
print("\n".join([f"{k}: {v}" for k, v in throw.items()][:12]) + "\n...")

## 1111: 0.00390625
## 1112: 0.015625
## 1113: 0.015625
## 1114: 0.015625
## 1122: 0.0234375
## 1123: 0.046875
## 1124: 0.046875
## 1133: 0.0234375
## 1134: 0.046875
## 1144: 0.0234375
## 1222: 0.015625
## 1223: 0.046875
## ...

In addition to these throw results, we have two states the chain can be in: "loss" and "win".

states = list(throw.keys()) + ["loss", "win"]
states

## ['1111', '1112', '1113', '1114', '1122', '1123', '1124', '1133', '1134', '1144', '1222', '1223', '1224', '1233', '1234', '1244', '1333', '1334', '1344', '1444', '2222', '2223', '2224', '2233', '2234', '2244', '2333', '2334', '2344', '2444', '3333', '3334', '3344', '3444', '4444', 'loss', 'win']

Next, we calculate the probability of going from one state to another. The following function will calculate all transition probabilities starting from state from_state. Remember we can roll every die in the “unique” group, but cannot change the die in the “duplicate” group. The states for "loss" and "win" are absorbing states.

def transition_probs(from_state):
  # Return all transition probabilities from state `from_state`
  
  # initially fill all probabilities with zero
  probs = {state: 0 for state in states}
  
  counts = Counter(from_state)
  uniques = [k for k, v in counts.items() if v == 1]
  
  if len(uniques) == 0 or from_state == "loss":
    probs["loss"] = 1
    return probs
  elif len(uniques) == 4 or from_state == "win":
    probs["win"] = 1
    return probs
  
  # number of dice to throw again
  n_dice = 4 - len(uniques) 
  # calculate possible new states
  draws = ["".join(str(i) for i in state) for state in product([1,2,3,4], repeat=n_dice)]
  new_states = ["".join(sorted(draw + "".join(uniques))) for draw in draws]
  # calculate new states' probabilities
  freqs = Counter(new_states).items()
  throw_probs = {k: v/(4**n_dice) for k, v in freqs}
  
  probs.update(throw_probs)
  
  return probs

For example, if we start with "1112", we throw the three "1" dice again, but keep the die with a "2". Therefore the probability mass is distributed over throws containing a "2".

transition_probs("1112")

## {'1111': 0, '1112': 0.015625, '1113': 0, '1114': 0, '1122': 0.046875, '1123': 0.046875, '1124': 0.046875, '1133': 0, '1134': 0, '1144': 0, '1222': 0.046875, '1223': 0.09375, '1224': 0.09375, '1233': 0.046875, '1234': 0.09375, '1244': 0.046875, '1333': 0, '1334': 0, '1344': 0, '1444': 0, '2222': 0.015625, '2223': 0.046875, '2224': 0.046875, '2233': 0.046875, '2234': 0.09375, '2244': 0.046875, '2333': 0.015625, '2334': 0.046875, '2344': 0.046875, '2444': 0.015625, '3333': 0, '3334': 0, '3344': 0, '3444': 0, '4444': 0, 'loss': 0, 'win': 0}

Next, we can generate a Markov matrix containing these transition probabilities for every starting state.

markov_mat = [[v for k, v in transition_probs(state).items()] for state in states]
markov_mat = np.asarray(markov_mat) # turn into numpy array for matrix multiplication
markov_mat.shape

## (37, 37)

Together with an initial throw, we can then simulate a round of the game by multiplying the vector corresponding to the initial throw with the Markov matrix..

first_throw = np.array([v for k, v in throw.items()] + [0, 0])

res = first_throw @ markov_mat
dict(zip(states, res))

## {'1111': 0.000732421875, '1112': 0.0087890625, '1113': 0.0087890625, '1114': 0.0087890625, '1122': 0.01611328125, '1123': 0.0380859375, '1124': 0.0380859375, '1133': 0.01611328125, '1134': 0.0380859375, '1144': 0.01611328125, '1222': 0.0087890625, '1223': 0.0380859375, '1224': 0.0380859375, '1233': 0.0380859375, '1234': 0.087890625, '1244': 0.0380859375, '1333': 0.0087890625, '1334': 0.0380859375, '1344': 0.0380859375, '1444': 0.0087890625, '2222': 0.000732421875, '2223': 0.0087890625, '2224': 0.0087890625, '2233': 0.01611328125, '2234': 0.0380859375, '2244': 0.01611328125, '2333': 0.0087890625, '2334': 0.0380859375, '2344': 0.0380859375, '2444': 0.0087890625, '3333': 0.000732421875, '3334': 0.0087890625, '3344': 0.01611328125, '3444': 0.0087890625, '4444': 0.000732421875, 'loss': 0.15625, 'win': 0.09375}

We can just as easily simulate multiple rounds of the game by taking the matrix to a certain power.
After 20 rounds, it becomes apparent that the probability of winning seems to be 45%.

res = first_throw @ np.linalg.matrix_power(markov_mat, 20)
dict(zip(states, res))

## {'1111': 2.2522812066292564e-06, '1112': 3.60364993060681e-05, '1113': 3.60364993060681e-05, '1114': 3.60364993060681e-05, '1122': 6.756843619887768e-05, '1123': 0.00016216424687730643, '1124': 0.00016216424687730645, '1133': 6.756843619887768e-05, '1134': 0.00016216424687730648, '1144': 6.756843619887768e-05, '1222': 3.603649930606811e-05, '1223': 0.00016216424687730648, '1224': 0.00016216424687730648, '1233': 0.00016216424687730648, '1234': 0.0003783832427137152, '1244': 0.00016216424687730645, '1333': 3.603649930606811e-05, '1334': 0.00016216424687730648, '1344': 0.0001621642468773065, '1444': 3.603649930606811e-05, '2222': 2.2522812066292564e-06, '2223': 3.60364993060681e-05, '2224': 3.60364993060681e-05, '2233': 6.756843619887769e-05, '2234': 0.0001621642468773065, '2244': 6.756843619887769e-05, '2333': 3.6036499306068095e-05, '2334': 0.00016216424687730648, '2344': 0.00016216424687730648, '2444': 3.6036499306068095e-05, '3333': 2.2522812066292564e-06, '3334': 3.60364993060681e-05, '3344': 6.756843619887768e-05, '3444': 3.60364993060681e-05, '4444': 2.2522812066292564e-06, 'loss': 0.5483423210319212, 'win': 0.4484864670291449}

Simulating the game for 1000 rounds confirms this.

res = first_throw @ np.linalg.matrix_power(markov_mat, 1000)
dict(zip(states, res))

## {'1111': 8.17825667702296e-129, '1112': 1.3085210683236735e-127, '1113': 1.3085210683236735e-127, '1114': 1.3085210683236735e-127, '1122': 2.4534770031068862e-127, '1123': 5.888344807456528e-127, '1124': 5.888344807456528e-127, '1133': 2.4534770031068862e-127, '1134': 5.888344807456528e-127, '1144': 2.4534770031068862e-127, '1222': 1.3085210683236735e-127, '1223': 5.888344807456528e-127, '1224': 5.888344807456528e-127, '1233': 5.888344807456528e-127, '1234': 1.3739471217398578e-126, '1244': 5.888344807456528e-127, '1333': 1.3085210683236735e-127, '1334': 5.888344807456528e-127, '1344': 5.888344807456528e-127, '1444': 1.3085210683236735e-127, '2222': 8.17825667702296e-129, '2223': 1.3085210683236735e-127, '2224': 1.3085210683236735e-127, '2233': 2.4534770031068862e-127, '2234': 5.888344807456528e-127, '2244': 2.4534770031068862e-127, '2333': 1.3085210683236735e-127, '2334': 5.888344807456528e-127, '2344': 5.888344807456528e-127, '2444': 1.3085210683236735e-127, '3333': 8.17825667702296e-129, '3334': 1.3085210683236735e-127, '3344': 2.4534770031068862e-127, '3444': 1.3085210683236735e-127, '4444': 8.17825667702296e-129, 'loss': 0.55, 'win': 0.44999999999999984}

What’s interesting is that you have a higher probability of winning (48.3%) if you start with a pair.

dict(zip(states, np.linalg.matrix_power(markov_mat, 1000)[:,-1]))

## {'1111': 0.0, '1112': 0.44999999999999996, '1113': 0.44999999999999996, '1114': 0.44999999999999996, '1122': 0.0, '1123': 0.48333333333333334, '1124': 0.48333333333333334, '1133': 0.0, '1134': 0.48333333333333334, '1144': 0.0, '1222': 0.44999999999999996, '1223': 0.48333333333333334, '1224': 0.48333333333333334, '1233': 0.48333333333333334, '1234': 1.0, '1244': 0.48333333333333334, '1333': 0.44999999999999996, '1334': 0.48333333333333334, '1344': 0.48333333333333334, '1444': 0.44999999999999996, '2222': 0.0, '2223': 0.44999999999999996, '2224': 0.44999999999999996, '2233': 0.0, '2234': 0.48333333333333334, '2244': 0.0, '2333': 0.44999999999999996, '2334': 0.48333333333333334, '2344': 0.48333333333333334, '2444': 0.44999999999999996, '3333': 0.0, '3334': 0.44999999999999996, '3344': 0.0, '3444': 0.44999999999999996, '4444': 0.0, 'loss': 0.0, 'win': 1.0}

And that is how easy and elegant it can be to answer such a question using a Markov chain. Until next time!