Background and Motivation
The game of Yahtzee is a model for decision making under uncertainty: specifically, for how to make optimal decisions when confronted with a sequence of them over time. The techniques used to analyze this game apply directly, for example, to business decisions that recur, such as whether and how much to make periodic investments in a business prospect. Thus, some time and effort spent in understanding this simple game can provide insights into developing strategies in complex, real-world situations.
Decision Analysis
Abstractly, a "prospect" is a set $\Omega$ of possible "states" of a game, business, or other situation. Associated with each state $\omega\in\Omega$ is a payoff $X(\omega)$ (a number which may be positive or negative). In the simplest case of a single decision, a set of actions or "options" $\mathcal{D}(\omega)$ is available to the decision maker; they may depend on the state $\omega$. Each option $d\in\mathcal{D}(\omega)$ leads to a probability distribution over the states $F(d)$ (which is either known or can be estimated accurately). Thus, associated with each option $d$ is the expected value of the payoff,
$$v(d) = \mathbb{E}_{F(d)}(X).$$
The best options are those with the highest expected values. Let $d^{*}(\omega)$ be any one of the best options for any state $\omega$. The function
$$X^\prime(\omega) = v(d^{*}(\omega))$$
that assigns to each state the value of choosing a best option is the payoff available under optimal decision making.
(A generalization of this theory accounts for risk by replacing expected payoff by expected utility.)
This figure diagrams a prospect consisting of two states named "11" and "001", shown in the right hand column. Associated with them are payoffs $1$ and $0$, respectively, indicating the payoff $X$. The possible starting states, shown in the left hand column, are the same two states. The set of all possible options is shown in the middle stack of orange nodes. Not all options are available from each of the beginning states: those that are, are shown with arrows. For instance, option "2" is not available from state "001" but it is available from state "11".
The probability distribution $F(d)$ associated with each option $d$ is shown by arrows pointing to the states: the chance of each state is given by a number on each arrow. (Arrows with no numbers have probability $1$.) The figure uses heavier arrows to show greater chances. The expectation of each distribution has been worked out and indicated next to each option. For instance, the value of option "1" is $3/4$ because its probability distribution assigns a chance of $3/4$ to reach state "11", whose payoff is $1$, and a chance of $1/4$ to reach "001", whose payoff is $0$. The expectation therefore is $3/4\times 1 + 1/4\times 0 = 3/4$, as shown.
The best options available from the beginning states are shown with heavy arrows. For instance, the four options available from state "001" are "0", "1", "01", and "001", with expected payoffs of $3/4$, $3/4$, $1/2$, and $0$, respectively. The first two are maximal, so the best options are "0" and "1".
In this fashion such an influence diagram contains all the information needed to work out the best decision to make in any beginning state, as well as its expected value. This one shows that $X^\prime$ assigns the value $1$ to "11" and $3/4$ to "001".
Now consider a prospect that is a sequence of decisions which will take place at discrete times indexed by $i=1, 2, \ldots, n$. The decisions and their circumstances may vary over time, so we must index everything (including the states and the options available from them). Optimal decisions can be made by "back-propagation" through this system, beginning at the end. At some time $i$, $i=2, 3, \ldots, n$, assume the best payoff $X_i$ is known. Then $X_{i-1} = X_i^\prime$ is the best payoff that will be available at time $i-1$. Iterating from $X_n = X$ down to $X_1$ produces a sequence of payoffs (random variables) along with a sequence of optimal decisions to make in order to achieve those best payoffs. Finally, the value of this prospect is the expectation of $X_1$ over a probability distribution governing the starting state at step $i=1$.
Computation
The calculations are amazingly simple, because they consist of a sequence of (a) finding an expectation for each option for each state and (b) selecting a best option for each state, thereby updating $X$ to $X^\prime$. However, for large problems they can get threateningly complicated: if there are many options for each of many states, this can be a lot of work. (Monte-Carlo methods come to the fore in such cases: when we cannot do the full calculation, we can sample from the possibilities and at least find approximately best options.)
In the present case, there is so much commonality among the options (regardless of the state) and so much similarity of the situation from one roll to the next, that with some care in designing computational data structures the entire calculation can be speedily done. Instead of maintaining a set of options for each state, the union of all the options $\mathcal{D}=\cup_{\omega\in\Omega} \mathcal{D}(\omega)$ can be computed once and for all. We then only have to record which elements of $\mathcal {D}$ are available for each state. This suggests writing a universal updating function decide which will compute $X^\prime$ from $X$. Its arguments will be the states $\Omega$ (which in Yahtzee are the same at each roll and so do not need to be indexed), the payoff vector $X$, the options $\mathcal{D}$, and a logical matrix available indicating the options available at each state.
As an example, here is an R implementation. It assumes that states and options can be used to index into the arrays as appropriate. When more than one option is best, it chooses randomly from the best ones (using sample). If we assume the existence of a function argmax that identifies the highest value in an array indexed by a set of options and returns that value and one of the options associated with it, the code needs all of two lines:
values <- options %*% payoff # Expected payoff for each option
x <- lapply(states, argmax) # Best option and its value, for each state
The first line (a fast matrix multiplication) computes expectations for all options in $\mathcal{D}$, then the second line applies argmax to all the states in $\Omega$ to obtain $X^\prime$ along with the associated best options. You can find these lines, as well as an implementation of argmax, in the decide function listed at the end of this post.
Solution
The rest is a matter of computing the states and options for Yahtzee, along with the function $F(d)$ giving the probability distribution of states for every option. This can be a real pain because the calculations are not simple. Three ideas make the computation relatively simple:
The states, options, and probabilities can be found simply by rolling five dice in sequence and keeping track of all outcomes and probabilities as you go along. No combinatorial mathematics is needed at all.
The options available from a state are precisely those where the state has a nonzero probability of resulting from the option, as given in (1).
The computation can be made highly efficient by adopting a suitable notation for the states.
In this case, the actual values showing on the dice don't matter: all we care about is how many unique values there are, how many values show twice, and so on. A good data structure for this information is a variable-length vector whose $i^\text{th}$ entry counts the number of values showing on exactly $i$ dice. We can terminate the vector after the last non-zero entry. To illustrate, the seven possible states in the game of Yahtzee are
$$(5),(31),(12),(201),(011),(1001),(00001).$$
(For convenience I have omitted the commas between digits. From now on, I will omit the parentheses, too.) The first state, $5$, represents five distinct values, such as $1, 3, 4, 5, 6$ showing on the dice. The state $011$ is a full house, such as $4,4,1,1,1$ showing on the dice. The state $00001$ is five of a kind: the yahtzee.
Options can be written using the same notation if we indicate what is showing on any dice that are not re-rolled. For example, suppose the initial roll shows $1,3,3,5,6$. Its code is $(31)$ because there are three unique values $1,5,6$ and one doubled value $3,3$. If we re-roll $1,5,6$, then the dice that remain can be indicated by the code $(01)$. There are 19 options in total, including re-rolling all dice and the seven states corresponding to re-rolling no dice. Incidentally, it is convenient to break the rule that all vectors end with a nonzero entry and write $(0)$ for the state where there are no dice: that makes all vectors non-empty.
Computing the probability distribution $F(d)$ associated with an option $d$ is the hardest part of the problem. This notation makes it easy to figure out what happens when a single die is rolled. The result is either (a) to increment one of the nonzero entries in the option vector by one, when the roll coincides with one of the faces already showing, or (b) to increment the first entry in the vector by one, when the roll shows a distinct new face. In fact, (b) can be considered a special case of (a) if we include an (invisible) zeroth entry in each option vector equal to the number of faces not yet appearing. For example, starting at the option $01$ in which the player has picked up three dice and left two of the same face showing, we examine the augmented vector $(5,0,1)$ (the $5$ is the number of faces not appearing among the two dice). It can result in the vector $(4,1,1)$ in exactly five different ways and in the vector $(5,0,0,1)$ in one way (the new roll matches the faces that were showing).
Generally, starting with the option $d = (d_1,d_2,\ldots,d_k)$ we may compute the following:
Construct the augmented vector $(d_0,d_1,d_2,\ldots,d_k)$ where $d_0$ is the number of faces $f=6$ minus the number of nonzero $d_i$.
For each position $i$ where $d_i \ne 0$, move to the vector $(d_0,d_1,\ldots,d_{i-1},d_i-1,d_{i+1}+1,d_{i+2},\ldots, d_k)$ with probability $d_i/f$. (For $i=k$, view $d_{i+1} = d_{k+1} = 0$.)
By performing this (very simple calculation) repeatedly until all five dice have been rolled, we will have determined all possible states reachable from $d$ along with the chances of those states, $F(d)$. For details, see the outcomes function in the code.
Answers
The chance of achieving a yahtzee is the expected value of the yahtzee indicator $X$ where $X(\text{yahtzee})=1$ and $X(\omega)=0$ otherwise.
The command yahtzee() (using the R code at the end of this post) produces extensive output showing all relevant data structures and intermediate results. Primary among these is the sequence of payoffs $X_i$:
$payoffs
5 31 12 201 011 1001 00001
After 1 0.0126314586 0.02906379 0.02906379 0.09336420 0.09336420 0.3055556 1
After 2 0.0007716049 0.00462963 0.00462963 0.02777778 0.02777778 0.1666667 1
After 3 0.0000000000 0.00000000 0.00000000 0.00000000 0.00000000 0.0000000 1
(This agrees with the table at http://en.wikipedia.org/wiki/Yahtzee#Any_Yahtzee. That table collapses "31" and "12" into a common category and "201" and "011" into another common category.)
The line at the bottom is the yahtzee indicator $X_3$: the state "00001" is five-of-a-kind, the yahtzee. The penultimate line is $X_2 = X_3^\prime$: these are the expected payoffs when the best options are chosen in each state. The first line is $X_1 = X_2^\prime$. Given the outcome of the initial roll of five dice, it tells us the expected payoff (with best play) for each possibility, one per column. For instance, if there is one pair showing after the first roll (state "31"), then the chance of achieving a Yahtzee after the third roll is $0.029\ldots$ with best play.
What is that best play? This is contained in the moves data structure of the output:
$moves
5 31 12 201 011 1001 00001
After 1 "0" "01" "01" "001" "001" "0001" "00001"
After 2 "0" "01" "01" "001" "001" "0001" "00001"
Note that these give a best move from each state, not necessarily all best moves, and therefore this table can change from one run to the next. (There's another best move to "1" from state "5", which is intuitively obvious: rerolling all five dice is the same as re-rolling any four of them.) Nevertheless, it shows that the strategy does not vary from one roll to another and that it amounts to the following:
If all five values are different, re-roll all the dice (or, equivalently, all dice but one).
Otherwise, retain only one of the faces with highest multiplicity and reroll all the others. (In particular, the best option from "12", which is two pairs, is to retain just one of those pairs and re-roll three dice rather than re-rolling the one odd die.)
Finally, the expected value can be found by consulting the probability distribution $F(d)$ associated with option "0" (no faces showing):
$options
5 31 12 201 011 1001 00001
0 0.09259259 0.4629630 0.2314815 0.15432099 0.03858025 0.019290123 0.0007716049
When multiplied by the first line of the payoffs matrix ("After 1") and added up, this gives the chance of achieving a yahtzee with best play. It is equal to $2,783,176 / 6^{10} \approx 0.04602864.$
Incidentally, it is not always the case that a static strategy is the best. If you want to achieve a full house "011", the best options actually change between the first and second rolls of the dice: try it and see!
Working Code
The following implementation is in R. I chose R for its popularity on this site (and its free availability), but it's a problematic language for working with any kind of dynamic data structure. It probably has caused many bugs to lurk which I haven't yet ferreted out--but it does work correctly on smaller problems that can be hand-checked and performs well when stepped through in debugging mode on larger problems.
This implementation handles arbitrary numbers of faces (all equiprobable), numbers of dice, numbers of rolls per move, and target states. (It is easily modified to solve a problem where the target is an arbitrary probability distribution of states, at some cost in flexibility: after all, the code itself computes what those states are, so it's not simple to design a way to input a probability distribution over them.) As a test, the command
yahtzee("11", 2, 3, 3)
finds optimal play to achieve one pair out of three fair coins, given three rolls. Its output includes
$payoffs
11 001
After 1 1 0.9375
After 2 1 0.7500
After 3 1 0.0000
$moves
11 001
After 1 "2" "1"
After 2 "2" "0"
The last two lines of payoffs and the last line of moves summarize the calculations in the figure above. Examining the rest of the output and comparing it to the figure ought to clarify the data structures and help understand the code.
#
# Perform one step of a decision analysis.
#
decide <- function(states, payoff, options, available) {
#
# Find the best option (and its value) from state `state` given a payoff
# vector `payoff` indexed by all states.
#
# Rows of `options` are probability distributions over the states
# Columns of `available` are logical indicators of which options are
# available from each state.
#
argmax <- function(state) {
valid <- available[, state]==TRUE
expectations <- values[valid]
names(expectations) <- rownames(available)[valid]
#i <- which.max(expectations) # Returns definite results
best <- which(expectations == max(expectations))
if (length(best) > 1) best <- sample(best, 1)
return (list(option=names(expectations)[best], value=expectations[best]))
}
#
# Return a list of best (option, value) pairs indexed by the states.
#
values <- options %*% payoff # Expected payoff for each option
x <- lapply(states, argmax) # Best option and its value, for each state
names(x) <- states
return (x)
}
#==============================================================================#
#
# Generate all states and options by rolling dice.
#
yahtzee <- function(target, n.faces=6, n.dice=5, n.rolls=3) {
#
# R-specific helper functions.
# These manage some data types and data structures.
#
state.start <- "0"
identity.matrix <- function(inames) {
x <- outer(inames, inames, '==')
rownames(x) <- colnames(x) <- inames
return (x)
}
df.as.matrix <- function(x) {
rows <- unique(x$start) #$
cols <- unique(x$end) #$
a <- matrix(0, length(rows), length(cols), dimnames=list(rows, cols))
for (i in 1:nrow(x)) a[x$start[[i]], x$end[[i]]] <- x$prob[[i]] #$
return (a)
}
state.name <- function(s) {
n <- max(which(s != 0));
ifelse(n > 0, paste(s <- s[1:n], collapse=""), state.start)
}
name.state <- function(x) as.numeric(unlist(strsplit(x, "")))
#
# Function to compute the transition probabilities among all possible options.
#
outcomes <- function(s, f=n.faces) {
x <- name.state(s)
transition <- function(i) {
r <- (x - (1:n == i) + (1:n == i+1))
p <- x[i]
return(list(start=s, end=state.name(r[-1]), prob=p/f))
}
x <- c(f - sum(x), x, 0)
n <- length(x)
result <- as.data.frame(t(sapply(which(x > 0), transition)))
return(result)
}
#----------------------------------------------------------------------------#
#
# Verify that `target` is legitimate.
#
if (!missing(target)) {
s <- name.state(target)
if (sum(s != 0) > n.faces || sum(s * 1:(length(s))) != n.dice)
stop(paste0("Target `", target, "` is inconsistent with ",
n.faces, " faces and ", n.dice, " dice."))
}
#
# Roll dice to find all possible states and the transition probabilities
# between them.
#
states <- state.start
x <- list()
empty <- subset(outcomes(states, n.faces), subset=FALSE)
for (die in 1:n.dice) {
x[[die]] <- empty
for (s in states) {
x[[die]] <- rbind(x[[die]], outcomes(s, n.faces))
}
states <- unlist(unique(x[[die]]$end)) #$
}
#
# Convert the transition graph `x` into an options matrix.
#
x <- lapply(x[1:n.dice], df.as.matrix)
options <- x[[1]]
for (i in 2:n.dice) options <- rbind(options %*% x[[i]], x[[i]])
options <- rbind(options, identity.matrix(colnames(x[[n.dice]])))
available <- ifelse(options != 0, 1, 0)
#
# Allocate arrays for the results.
#
payoffs <- matrix(NA, n.rolls, length(states))
colnames(payoffs) <- states
rownames(payoffs) <- sapply(1:n.rolls, function(i) paste("After", i))
moves <- payoffs[-n.rolls, ]
#
# Compute the final payoff vector.
#
if (missing(target)) target <- paste(c(rep(0, n.dice-1), 1), collapse="")
payoffs[n.rolls, ] <- 0
payoffs[n.rolls, target] <- 1
#
# Perform the decision analysis, roll by roll.
#
if (n.rolls > 1) {
for (i in (n.rolls-1):1) {
x <- decide(states, payoffs[i+1, ], options, available)
payoffs[i, ] <- unlist(lapply(x, function(s) s$value))
moves[i, ] <- unlist(lapply(x, function(s) s$option))
}
}
#
# Compute the expectation based on the probabilitie of the starting state
# (that is, of the outcome of the first roll).
#
expectation <- as.vector(options["0", , drop=FALSE] %*% payoffs[1, ])
return (list(expectation=expectation, payoffs=payoffs, moves=moves,
options=options, available=available,
target=target, parameters=c(faces=n.faces, dice=n.dice, rolls=n.rolls)))
}
