Probability of throwing the same number on 5 dice after 5 throws while removing the most frequently occuring number

Question

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)
    else:
        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
    print(shift)
    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
getProbabilities(max_number_of_attempts)

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

Answer 1

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)
    else:
        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
        else:
            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
            else:
                next_turn[next_dice] = odds2
    turns.append(next_turn)

for i in range(6):
    dice = tuple(i+1 for _ in range(5))
    print (dice, turns[-1][dice])

Answer 2

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.