Imagine a game similar to yahtzee, where the goal is to get roll 5 times the same number with 5 dice. Which number it is, doesn't matter, however the rules are a little different:

  1. You get a maximum of 5 tries.
  2. After each try, you must remove the dice which show the number occurring the most. (e.g. if you throw 1,1,2,4,5, you remove both 1's and continue with 3 dice; if you throw 1,1,2,2,5 the rule is to pick the highest number of the most occurring so keep 2,2 and throw 3 in the next turn).

This seems pretty straightforward but imagine the following situation:

  1. Throw 1,2,4,5,6: keep 6 and continue with 4 dice
  2. Throw 2,2,3,5 + 6(kept): keep 2,2 (as it occurs more than 6) and continue with 3 dice
  3. Keep throwing until you've reached either 5 of the same number or you run of of tries

So the value of the dice you're currently keeping out of the game can shift according to your throw. The question is: what is the probability of succeeding in such a game in a maximum of 5 tries

I've written a small script in Python with simulates the game but I can't figure out how to calculate it using exact probabilities. A run of my script with 1 million samples yields a probability of approximately 17%

EDIT: Thanks for all the insight guys. I eventually went for a Markov-Chain solution which takes into account the probabilities of shifting states:

import math as m

def binomial(n,k,p):
    if n>=k:
        n_k = m.factorial(n)/(m.factorial(k)*m.factorial(n-k))
        p_k = p**k
        binom = n_k*p_k*(1-p)**(n-k)
        binom = 0
    return binom

def initialState():
    state_1 = 5*4*3*2/(6**4) #independent of first throw make a different roll each time
    state_3 = 6*binomial(5,3,1/6) #exactly 3 successes in 5
    state_4 = 6*binomial(5,4,1/6) #exactly 4 successes in 5
    state_5 = 6*binomial(5,5,1/6) #exactly 5 successes in 5
    state_2 = 1-state_1-state_3-state_4-state_5
    return [state_1,state_2,state_3,state_4,state_5]

def stateShiftProbability():
    shift = [[0]*5 for _ in range(5)]
    #chance of state 1 staying 1 = none of the 4 thrown dice match the original
    # and we throw no new pairs or more
    shift[0][0] = 5*4*3*2/(6**4)
    #chance of state 1 to state 3 = chance to throw 2 more of the same value or
    # 3 of a new value
    shift[0][2] = binomial(4,2,1/6)+5*binomial(4,3,1/6)
    shift[0][3] = binomial(4,3,1/6)+5*binomial(4,4,1/6)
    shift[0][4] = binomial(4,4,1/6)
    #chance of state 1 to state 2 = 1 of the 4 thrown dice match the original
    # + the chance that we throw a new pair (given that the other 2 dice aren't
    # equal to the orignal dice). Quite complex, so we just take the remaining
    # value
    shift[0][1] = 1-sum(shift[0])
    #from state 2 we follow a similar logic
    shift[1][2] = binomial(3,1,1/6)+5*binomial(3,3,1/6)
    shift[1][3] = binomial(3,2,1/6)
    shift[1][4] = binomial(3,3,1/6)
    shift[1][1] = 1-sum(shift[1])

    shift[2][2] = binomial(2,0,1/6)
    shift[2][3] = binomial(2,1,1/6)
    shift[2][4] = binomial(2,2,1/6)

    shift[3][3] = 5/6
    shift[3][4] = 1/6

    shift[4][4] = 1
    return shift

def getProbabilities(tries):
    initial_state = initialState()
    state_shift = stateShiftProbability()
    current_state = initial_state[:]
    print("After",1,"throw, we are successful", current_state[4]*100,"% of the time")
    for i in range(tries-1):
        new_state = [0]*5
        for j in range(5): #every state j
            # a new state is formed by the probability of a state remaining
            # unchanged + new population shifting towards this state
            for k in range(5):
                new_state[j] += state_shift[k][j]*current_state[k]
        current_state = new_state[:]
        print("After",i+2,"throws, we are successful", current_state[4]*100,"% of the time")

max_number_of_attempts = 5

This results in the exact probabilities I wanted:

After 1 throw, we are successful 0.07716049382716046 % of the time
After 2 throws, we are successful 1.2631458619112939 % of the time
After 3 throws, we are successful 4.6028642525699 % of the time
After 4 throws, we are successful 10.057536103294314 % of the time
After 5 throws, we are successful 17.05067236780467 % of the time
  • 1
    $\begingroup$ rather than a simulation, you could use a computer to trawl through all possible dice combinations, and do a complete calculation. All of your first possible rolls could become objects that would then be populated with the next set of all conditionally possible rolls - it's deterministic, there are no decisions to make $\endgroup$ – Cato Mar 14 '17 at 16:45
  • $\begingroup$ @Cato True, this gives the exact probability. $\endgroup$ – Peter Mar 14 '17 at 16:46
  • $\begingroup$ 'So the value of the dice you're currently keeping out of the game can shift according to your throw' how, do you only keep them out for one round? If not then if you throw a pair, then a second different pair ends the game for you - is that right? $\endgroup$ – Cato Mar 14 '17 at 16:53
  • $\begingroup$ you simulation should be about plus or minus .2% if my statistics memory serves me correctly $\endgroup$ – Cato Mar 14 '17 at 17:01
  • $\begingroup$ @Cato the dice you keep out of the game are theoretically only removed for one round. They just count towards the total for the next throw. So if you keep 2 out and throw a three of a kind of another number, on the next throw you will keep the new triple out and only throw the original pair. The game ONLY ends when you throw 5 in 5 turns or less $\endgroup$ – Jonas Coussement Mar 14 '17 at 17:02

The odds are approximately 17.05067...% so your simulation is pretty accurate. Here are exact odds for the 6 winning rolls:

((1, 1, 1, 1, 1), Fraction(451727069, 19591041024)) ((2, 2, 2, 2, 2), Fraction(1208621087, 58773123072)) ((3, 3, 3, 3, 3), Fraction(1516264807, 58773123072)) ((4, 4, 4, 4, 4), Fraction(1693008047, 58773123072)) ((5, 5, 5, 5, 5), Fraction(17336750593, 528958107648)) ((6, 6, 6, 6, 6), Fraction(20896486973, 528958107648))

And here is the program that I used to calculate it:

import fractions

def die_roll ():
    odds = fractions.Fraction(1, 6)
    for i in range(6):
        yield (i+1, odds);

def roll_dice (n):
    if 0 == n:
        yield(tuple(), 1)
    elif  1 == n:
        for i, odds in die_roll():
            yield ((i,), odds)
        for prev_dice, odds in roll_dice(n-1):
            for i, odds2 in die_roll():
                yield (prev_dice + (i,), odds * odds2)

def roll_dice_grouped (n):
    answer = {}
    for dice_u, odds in roll_dice(n):
        dice = tuple(sorted(dice_u))
        if dice in answer:
            answer[dice] = answer[dice] + odds
            answer[dice] = odds
    return answer

def reroll_dice (dice, odds):
    saw = {1: 0, 2: 0, 3: 0, 4: 0, 5: 0, 6: 0}
    best = {'v': 0, 'c': 0}
    for d in dice:
        saw[d] = saw[d] + 1
        if best['c'] < saw[d]:
            best['c'] = saw[d]
            best['v'] = d
        elif best['c'] == saw[d] and best['c'] < d:
            best['c'] = saw[d]
            best['v'] = d
    held_dice = tuple(best['v'] for _ in range(best['c']))
    for d, odds2 in roll_dice(5 - best['c']):
        yield (d + held_dice, odds * odds2)

turns = [roll_dice_grouped(5)]
for i in range(4):
    this_turn = turns[-1]
    print ("finished turn", i+1)
    next_turn = {}
    for dice in this_turn:
        odds = this_turn[dice]
        for next_dice_u, odds2 in reroll_dice(dice, odds):
            next_dice = tuple(sorted(next_dice_u))
            if next_dice in next_turn:
                next_turn[next_dice] = next_turn[next_dice] + odds2
                next_turn[next_dice] = odds2

for i in range(6):
    dice = tuple(i+1 for _ in range(5))
    print (dice, turns[-1][dice])
  • $\begingroup$ Very interesting solution. I've revised my Markov-chain solution and now I've got the exact same result. You've taken a completely different approach though, nice work! $\endgroup$ – Jonas Coussement Mar 14 '17 at 18:53

It will be a fair amount of work. You need to compute the chance of going from keeping $m$ dice to keeping $n$ dice on one roll. For example, the chance of going from keeping $2$ dice to keeping $3$ dice comes from the chance you get exactly one match of your current set and from the chance all three will match but not match the current set. This is $\frac 1{216}(3\cdot 5 \cdot 5+5)=\frac {80}{216}=\frac {10}{27}$ The $3 \cdot 5 \cdot 5$ comes from the three choices of which die matches and the five choices for each of the other dice. Once you have done this you can compute the chance of keeping each number of dice at each roll. If a decimal result is acceptable you can do it in a spreadsheet, but it will not be exact. For an exact answer you need to keep the fractions. Note that rolling four dice is just like rolling all five, so you can save that much work.

  • $\begingroup$ I was thinking the same thing. In essence I would probably need a markov-chain where the intial states are the (probabilities of the) amount of dice i get the keep after the first roll and then compute the probabilities of the state changes as you say. Then I can check the value of state 5 after 5 runs. Does this sound like a good approach to you? $\endgroup$ – Jonas Coussement Mar 14 '17 at 17:10
  • $\begingroup$ Yes, that is just what I was thinking $\endgroup$ – Ross Millikan Mar 14 '17 at 17:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.