You and four statistician colleagues find a $100 bill on the floor of your department’s faculty lounge. None of you have change, so you agree to play a game of chance to divide the money probabilistically. The five of you sit around a table. The game is played in turns. Each turn, one of three things can happen, each with an equal probability: The bill can move one position to the left, one position to the right, or the game ends and the person with the bill in front of him or her wins the game. You have tenure and seniority, so the bill starts in front of you. What are the chances you win the money?
There’s definitely a starting advantage for me (i.e., the professor with tenure and seniority). This problem can be modeled as a Markov process with absorbing states. At each turn, one of three things can happen. If either the $100 bill passes left or right, these are transitions to non-absorbing states, because the game does not end, and another turn follows. However, if a player wins, this is an absorbing state, as this ends the game. So, there are five possible absorbing states, one for the possibility of each player winning. And there are ten possible non-absorbing states, one for each combination of five professors passing either left or right. What we want to know is: What probability is associated with the absorbing state for the starting player, in this case, that of the senior professor?
We start by writing out one initial state matrix and two transition matrices: one matrix for transitions from non-absorbing states to non-absorbing states, and another matrix for transitions from non-absorbing states to absorbing states. For convenience, let’s number our professors from one to five, and let’s place the senior professor at position number three in the middle. So, because the bill starts with the senior professor at position three, the initial state may be written:
For transitions from non-absorbing states to non-absorbing states, let’s remember that the professors are sitting around a table, so the professors at positions one and five can pass the bill to one another. We can write out the first transition matrix :
So for example, in the first column of , the in the second and fifth rows indicates that there is a one-third probability – if the bill is with the first professor – the bill being passed to either the second or fifth professor. Note that columns sum only to , because for each turn there is also a one-third chance that the process enters into an absorbing (i.e., winning) state. Our second transition matrix – representing transitions to absorbing states – may then be written:
From here, we can obtain the fundamental matrix , where, is the inverse of , and is the identity matrix of the same dimensions as . We multiply by , because we are only interested in the outcome for the initial state, where the bill begins at position three with the senior professor.
One way to interpret this is the expected number of turns the bill is any any position. Unsurprisingly, the elements of this vector sum to three, and in expectation, the bill will spend the greatest number of turns in the third position. This is the same result we would get if we solved for the expected number of trials until first success, i.e., the expected number of turns until any professor won the bill, which is simply , where is the probability of success on each trial.
With the same transition matrices, we can also obtain the stable matrix , given the initial state , which is also the solution to the riddle:
The stable matrix gives us the probability of absorption in each position; note that these sum to one. Thus, the senior professor in the third position will win the $100 bill with probability 0.455; the professors to right and left of the senior professor both win the bill with probability 0.182; and the two professors on the other side of the table win the bill with probability 0.091.
Implementation in R
R is a great environment for matrix math, as long as one does not forget to enclose matrix operators with the “%” symbol. These matrices are just a little too large to comfortably work with on paper. Here we define our transition matrices and the identity matrix with the same five-by-five dimensions.
Simulation is a quick and easy way to check these results. Here, we simply use a number of nested loops to play this game 100,000 times. It has been said that usage of loops like this is “not R-like,” because R is vectorized, and there are usually more efficient approaches. I agree for the most part; however, for simple tasks like this, which require minimal computational power, I find nested loops more “human readable” and intuitive. Throwing in a few print() commands allows for easy debugging and allows one to “see” the simulation (or game, in this case) played out.
# for reproducability
set.seed(405)# to store results
out<-matrix(nrow=100000,ncol=2)# start sim
for(iin1:100000){cp<-3# reset game
winner<-0player<-3# start with senior prof in position 3
round<-0#start iteration of game
while(winner==0){#play until there is a winner
turn<-sample(1:3,1)round<-round+1if(turn==2){winner<-1}#end game
elseif(turn==1){#pass left
if(cp!=1){cp<-cp-1}elseif(cp==1){cp<-5}# b/c profs in circle
}elseif(turn==3){#pass right
if(cp!=5){cp<-cp+1}elseif(cp==5){cp<-1}# b/c profs in circle
}}# save game results
out[i,1]<-cpout[i,2]<-round}out<-data.frame(out)colnames(out)<-c("winner","num.rounds")# print results
print(round(table(out$winner)/100000,3))
##
## 1 2 3 4 5
## 0.091 0.181 0.454 0.183 0.091
print(mean(out$num.rounds))
## [1] 2.99247
In this simulation, the senior professor in the third position wins the $100 with probability 0.454 – about the same as we got with our transition matrices above. The 0.001 deviation is due partly to chance, and also because base R’s round() function rounds values of 0.5 to even numbers, rather than rounding up. Likewise, the average number of turns-per-game in the simulation was 2.99 turns, very close to the expected value of 3.00.
Plots of Simulation Results
library(ggplot2)ggplot(out,aes(x=factor(winner)))+geom_bar(aes(y=(..count..)/sum(..count..)),fill="#45B29D")+ylab("Probability of winning $100 bill\n")+xlab("\nProfessor position\n(Senior professor at position three)")+ggtitle("Simulated probability of winning,\ngiven 100,000 iterations\n")+theme(panel.background=element_rect(fill="#ffffff"))
ggplot(out,aes(x=num.rounds))+geom_bar(aes(y=(..count..)/sum(..count..)),fill="#45B29D")+ylab("Fraction of games won in number of turns\n")+xlab("\nNumber of turns")+ggtitle("Simulated number of turns until win,\ngiven 100,000 iterations\n")+theme(panel.background=element_rect(fill="#ffffff"))
RnotR PROBABILITY · EXPECTATION · TRANSITION PROBABILITIES riddler