Probability of Yahtzee straight with strict re-roll rules

Question

Assuming you have 5 fair dice and up to 3 re-roll attempts (during which you can re-roll any dice, as per standard Yahtzee rules) what is the probability of succesfully rolling a "Small Straight" (aka 1-2-3-4, 2-3-4-5, or 3-4-5-6)? Also, what is the probability of rolling a "Large Straight" (aka 1-2-3-4-5 or 2-3-4-5-6)?

In either case, assume the following rules for keeping and re-rolling dice:

• Keep exactly one of each 2, 3, 4, and 5 result
• Only keep a 1 or 6 result if it is already part of a Small Straight or Large Straight

Unfortunately I'm quite unfamiliar with even basic probability problems let alone something with this level of complexity, and online research hasn't helped much. I'm much more comfortable backing into probabilities by simulating roll outcomes like this in Python, though I've hit a wall with those efforts as well for straights in particular, even with these defined rules. Would deeply appreciate some help.

Thanks!

Answer 1

(I have edited the answer)

Large Straight

Let me see if I got the strategy rules right. Does the following simulation code (in Python) model the situation correctly?

import random

keepers = set([2,3,4,5])

def simuLS(doPrint=False):
    kept = set()
    for roundI in range(3):
        ds = [random.randint(1, 6) for _ in range(5-len(kept))]
        for d in ds:
            if d in keepers: kept.add(d)
        if 1 in ds and all(x in kept for x in [2,3,4]):
            kept.add(1)
        if 6 in ds and all(x in kept for x in [3,4,5]):
            kept.add(6)
        gotLS = len(kept)==5
        if doPrint:
            print ("Rolled", ds)
            print ("Kept are ", kept)
            print ("Got LS: ", gotLS)
        if gotLS: return True
    return False

#Let's see how it goes
#simuLS(doPrint=True)

simuN = 10000
p = sum(1.0 for _ in range(simuN) if simuLS())/simuN
print (float(p))

The output is around 0.24.

Are you familiar with Markov chains? We can construct a Markov chain with states all the combinations of the keeper-numbers $\{1,2,3,4,5,6\}$ at most of length 4 and having 1 or 6 only if 2,3,4 and 3,4,5 are present respectively. We also have one absorbing state 'LS' that means we have gotten a large straight. We can calculate the transition probabilities by brute force by considering all possible rolls of the not kept dice (there are 5 - length of the state of them) and seeing what additional keepers we got or did we in fact get a LS. I wrote the following code in SageMath, since it handles matrices and rational numbers without any efforts.

import itertools

keepers = set([2,3,4,5])

states = [tuple(S) for S in Subsets(range(1,7)) if len(S)<5 and
          (1 not in S or all(x in S for x in [2,3,4]))
          and (6 not in S or all(x in S for x in [3,4,5]))]
startState = tuple() #the empty set
targetState = 'LS' #have gotten large straight
states.append(targetState)
stateToInd = {s: i for i,s in enumerate(states)}
startI = stateToInd[startState]
targetI = stateToInd[targetState]
mat = matrix(QQ, len(states))

for i1,s1 in enumerate(states[:-1]): #don't do for the LS-state
    d = 5-len(s1) #how many dice are (re-)thrown
    s1Set = set(s1)
    p = 1/6**d #probability for each outcome
    for xs in itertools.product(*[range(1, 7) for _ in range(d)]):
        s2Set = s1Set.union(set(xs)&keepers)
        if 1 in xs and all(x in s2Set for x in [2,3,4]): s2Set.add(1)
        if 6 in xs and all(x in s2Set for x in [3,4,5]): s2Set.add(6)
        
        if len(s2Set)==5: #got LS
            mat.add_to_entry(i1, targetI, p)
        else:
            s2 = tuple(sorted(s2Set))
            mat.add_to_entry(i1, stateToInd[s2], p)
mat.add_to_entry(targetI, targetI, 1) #absorbing state

#look at the transition matrix
print (len(states))
print (states)
show(mat)

#print ([sum(row) for row in mat]) #sanity check

sol = (mat**3)[startI][targetI]
print (sol, "=", float(sol))

You can run Sage code online: Sage Cell Server

The output is $\frac{295073737}{1224440064} = 0.24098667274578822$.

I'm not sure the grouping together of states, that I mentioned in my answer before editing, can be done anymore, since now the transitions to the states containing 1 or 6 depend on what we already have. At least not for all same-size-groups. Anyhow, we have a $19\times 19$-matrix, so it is doable.

Small Straight

We can use the previous Markov chain also here:

startRow = (mat**3)[startI]
endsForSmall = [stateToInd[s] for s in states if len(s)>=4] + [targetI]
smallSol = sum(startRow[j] for j in endsForSmall)
print (smallSol, "=", float(smallSol))

Answer: $\frac{2966131613}{4897760256} = 0.605609800799527$

For simulation:

import random

keepers = set([2,3,4,5])

def simuSS(doPrint=False):
    targets = [(1,2,3,4),
               (2,3,4,5),
               (3,4,5,6)]
    kept = set()
    for roundI in range(3):
        ds = [random.randint(1, 6) for _ in range(5-len(kept))]
        for d in ds:
            if d in keepers: kept.add(d)
        if 1 in ds and all(x in kept for x in [2,3,4]):
            kept.add(1)
        if 6 in ds and all(x in kept for x in [3,4,5]):
            kept.add(6)
        gotSS = any(all(x in kept for x in y) for y in targets)
        if doPrint:
            print ("Rolled", ds)
            print ("Kept are ", kept)
            print ("Got SS: ", gotSS)
        if gotSS: return True
    return False

#Let's see how it goes
#simuSS(doPrint=True)

simuN = 10000
p = sum(1.0 for _ in range(simuN) if simuSS())/simuN
print (float(p))
#--> 0.60