Solved – Simulating the game of bingo using a Monte Carlo approach. Is this methodology correct

estimationmonte carloprobabilityrsimulation

A while ago while talking with a friend, he said a friend of his goes to play bingo in order to win money. I told him that this is not possible because his friend is going to lose in the long run, even though they could win some games.

I would like to prove (or disprove!) this using a Monte Carlo approach and R. In order to prove it I decided to simplify the game as follows:

The number of players is fixed (average value per evening)
The cost of each bingo card is fixed (average cost of a bingo card)
Each player buys exactly one bingo card per game (although this is irrelevant since buying more than one bingo card could be assimilated to having more players each with only one bingo card)
The gain/loss at each game per player is a random variable (we can call it X) so that:
$$ X =
\begin{cases}
k-c & \quad p\\
-c & \quad 1-p\\
\end{cases}
$$
where $c > 0$ is the cost of a game/bingo card, $k > 0$ is the prize if the player does a bingo ($k = \mathrm{number\ of\ players} \times c \times 0.53$) and $0 < p < 1$ is the probability of a bingo.

If we know $p$ then we can calculate the expected value of $X$. The problem is now only about estimating $p$.

My first question is: is there any major flaw in this simplification I made that could skew the estimate significantly? The simplification neglects some special prizes that are so rare to win (perhaps 0.01% chance or so) that I felt comfortable leaving them out as a first approximation. Read below for the second question.

The code I used to run the simulation is this one

generate_card <- function(n=15)
{
    card <- sample(1:90,size=15,replace=T)
    return(card)
}

game_df <- function(n_player)
{
    df <- data.frame()
    for(i in 1:n_player)
    {
        df <- rbind(df,generate_card())
    }
    rownames(df) <- as.character(1:n_player)
    return(df)
}


game <- function()
{
    players <- 150
    df <- game_df(players)
    bingo <- F
    numbers_ <- 1:90
    number_drawn_balls <- 0
    pwon <- NA      # Player that won
    while((!bingo) && (length(numbers_) > 0) )
    {
        drawn_number <- sample(numbers_,size=1)
        number_drawn_balls <- number_drawn_balls + 1
        #print(drawn_number)
        numbers_ <- numbers_[numbers_ != drawn_number]

        df[df==drawn_number] <- NA

        bingos <- as.vector(apply(df,1,function(x) sum(is.na(x))))

        for(i in 1:length(bingos))
        {
            if(bingos[i] == 15)
            {
                bingo <- T
                pwon <- i
            }
        }
    }
    return(c(pwon,number_drawn_balls))
}

wins <- NULL
number_of_draws <- NULL
for(k in 1:500)
{
    wins[k] <- game()[1]
    number_of_draws[k] <- game()[2]
}

By running 500 games I calculated the following

hist(number_of_draws,freq=F)
mean(number_of_draws)

hist(wins,freq=F,breaks=1:150)
mean(prop.table(table(wins)))

cost_card <- 3
p <- mean(prop.table(table(wins)))
print(cost_card*150*0.53*p - 3(1-p))

On average 62.5 balls need to be drawn before a player does a bingo.
The probability of a bingo is 0.69% on average, and therefore the expected value $E(X)=-1.30$. The first number (62.5) looks suspiciously high to me although it seems to be consistent with the special prices left out that are usually awarded when a bingo is achieved with number of balls less than say $l > 0$.
But the most suspicious thing to me is the histogram of wins per player (density).

It looks suspicious because I expected a flat histogram where each player should win on average the same number of times. I calculated the probability of a bingo as an average of the probability of a bingo for each player, however, by looking at each single player, some at times might even get a positive expected value and I do not know how to explain/account for this. Is the code right at all? Did I miss something?

Best Answer

Your "histogram of wins" appears to be just a bar chart of the proportion of wins (out of 500 games) for each of the 150 players (the horizontal axis is mislabelled). On average you'd expect each to win $500/150=3\frac{1}{3}$ times, so it's perhaps not surprising that four lucky players have won eight times, & a single unlucky player not at all. A simple way to check is to simulate a few times 500 draws from a multinomial distribution with equal probability across 150 categories (rmultinom(1,500,rep(1/150,150))) & see if your results look similar. A more formal way is to perform Pearson's chi-squared goodness-of-fit test. Note that the more games you simulate, the closer the proportion of wins for each player should get to the expected value.

What is strange, however, is that the simulated chance of winning isn't exactly $1/150=0.66\dot6$; if one, & only one, player wins each game. If two people complete their cards when a number's called, is it the first to shout "Bingo!" that wins? Or do they both win? Do they split the prize?

Related Solutions

Solved – How to predict the results of a simple card game

The easiest way is just to simulate the game lots of times. The R code below simulates a single game.

nplayers = 4
#Create an empty data frame to keep track
#of card number, suit and if it's magic
empty.hand = data.frame(number = numeric(52),
  suit = numeric(52),
  magic  = numeric(52))

#A list of players who are in the game
players =list()
for(i in 1:nplayers)
  players[[i]] = empty.hand

#Simulate shuffling the deck
deck = empty.hand
deck$number = rep(1:13, 4)
deck$suit = as.character(rep(c("H", "C", "S", "D"), each=13))
deck$magic = rep(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0), each=4)
deck = deck[sample(1:52, 52),]

#Deal out five cards per person
for(i in 1:length(players)){
  r = (5*i-4):(5*i)
  players[[i]][r,] = deck[r,]
}

#Play the game
i = 5*length(players)+1
current = deck[i,]
while(i < 53){
  for(j in 1:length(players)){
    playersdeck = players[[j]]
    #Need to test for magic and suit also - left as an exercise!
    if(is.element(current$number, playersdeck$number)){
      #Update current card
      current = playersdeck[match(current$number,
        playersdeck$number),]
      #Remove card from players deck
      playersdeck[match(current$number, playersdeck$number),] = c(0,
                   0, 0)
    } else {
      #Add card to players deck
      playersdeck[i,] = deck[i,]
      i = i + 1
    }
    players[[j]] = playersdeck
    #Has someone won or have we run out of card
    if(sum(playersdeck$number) == 0 | i > 52){
      i = 53
      break
    }
  }
}

#How many cards are left for each player
for(i in 1:length(players))
{
  cat(sum(players[[i]]$number !=0), "\n") 
}

Some comments

You will need to add a couple of lines for magic cards and suits, but data structure is already there. I presume you didn't want a complete solution? ;)
To estimate the average game length, just place the above code in a function and call lots of times.
Rather than dynamically increasing a vector when a player gets a card, I find it easier just to create a sparse data frame that is more than sufficient. In this case, each player has a data frame with 52 rows, which they will never fill (unless it's a 1 player game).
There is a small element of strategy with this game. What should you do if you can play more than one card. For example, if 7H comes up, and you have in your hand 7S, 8H and the JC. All three of these cards are "playable".

Probability – Expected Rolls to Roll Every Number on a Dice an Odd Number of Times Explained

You can think about your problem as a Markov chain, i.e., a set of states with certain transition probabilities between states. You start in one state (all cards face up) and end up in an absorbing state (all cards face down). Your question is about the expected number of steps until you reach that absorbing state, either for a single chain, or for the expected minimum number of steps across $n$ independent Markov chains running simultaneously.

And there are actually two slightly different ways to look at this. The first one, as whuber comments, is to consider the six cards as six different bits $\{0,1\}$ and consider the state as a six-vector in $\{0,1\}^6$, i.e., the six-dimensional discrete hypercube. We start out at the vertex $(0,0,0,0,0,0)$, and the absorbing state is $(1,1,1,1,1,1)$. A step can take us to an adjacent vertex, in which exactly one bit is flipped with respect to the original state. That is, transitions take us from one vertex to any neighboring one with Hamming distance exactly one, and each such neighbor has an equal probability of being the next state.

There is some literature on random walks and Markov chains on discrete cubes with Hamming distances, but nothing I could locate at short notice. We have a very nice thread on Random walk on the edges of a cube, which might be interesting.

The second way to look at this is to use the fact that all cards are interchangeable (assuming a fair die). Then we can use just seven different states, corresponding to the number of cards face down. We start in the state $i=0$, and the absorbing state is $i=6$. The transition probabilities depend on the state we are in:

From $i=0$ (all cards face up), we will flip one card down and end up with one card face down with certainty: we have the transition probability $p_{01}=1$ (and $p_{0j}=0$ for $j\neq 1$).
From $i=1$, we can reach $j=0$ with probability $p_{10}=\frac{1}{6}$ and $j=2$ with probability $p_{12}=\frac{5}{6}$.

Overall, we get the following transition matrix:

$$ T=\begin{pmatrix} 0 & \frac{6}{6} & 0 & 0 & 0 & 0 & 0 \\ \frac{1}{6} & 0 & \frac{5}{6} & 0 & 0 & 0 & 0 \\ 0 & \frac{2}{6} & 0 & \frac{4}{6} & 0 & 0 & 0 \\ 0 & 0 & \frac{3}{6} & 0 & \frac{3}{6} & 0 & 0 \\ 0 & 0 & 0 & \frac{4}{6} & 0 & \frac{2}{6} & 0 \\ 0 & 0 & 0 & 0 & \frac{5}{6} & 0 & \frac{1}{6} \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 \end{pmatrix} $$

We start with certainty in the state $i=0$. We can encode the probabilities for each state at a certain point with a vector $v\in[0,1]^7$, and our starting state corresponds to $v_0=(1,0,0,0,0,0,0)$.

Here is a fundamental fact about Markov chains (which is easy to see and to prove via induction): the probabilities for the state after $k$ transitions are given by $v_k=(T')^kv_0$. (That is $T$ transposed. You can also work with row vectors $v$, then you don't need to transpose, but "$v_0T^k$" takes a little getting used to.)

Thus, the probability that we have ended up in the absorbing state $i=6$ after $k$ steps is precisely the last entry in that vector, or $v_k[6]=((T')^kv_0)[6]$. Of course, we could already have been in the absorbing state after $k-1$ steps. So the probability that our Markov chain ends up in the absorbing state for the first time after $k$ steps is

$$ p_k := ((T')^kv_0)[6]-((T')^{k-1}v_0)[6]. $$

We can numerically calculate $p_k$ for a large enough number of $k\leq K$ such that $\sum_{k=0}^Kp_k\approx 1$, and there may even be a closed form solution. Then, given $p_k$, we can calculate the expectation as

$$ \sum_{k=0}^\infty kp_k \approx \sum_{k=0}^K kp_k. $$

Next, assume we have $n$ players, and we want to know after how many steps the game will end, i.e., when the first player has all their cards face down. We can easily calculate the probability $q_k^n$ that at least one player has all cards face down after $k$ or fewer steps by noting that

$$ \begin{align*} q_k^n &= P(\text{at least one player has all cards face down after $k$ or fewer steps}) \\ &= 1-P(\text{all $n$ players need at least $k+1$ steps}) \\ &= 1-P(\text{ONE player needs at least $k+1$ steps})^n \\ &= 1-\bigg(\sum_{j=k+1}^\infty p_j\bigg)^n \\ &= 1-\bigg(1-\sum_{j=0}^k p_j\bigg)^n. \end{align*} $$

From this, we can derive the probability $p^n_k$ that a game of $n$ players ends after exactly $k$ steps:

$$ p^n_k = q_k^n-q_{k-1}^n = \bigg(1-\sum_{j=0}^{k-1} p_j\bigg)^n-\bigg(1-\sum_{j=0}^k p_j\bigg)^n. $$

And from this, in turn, we can again calculate the expected length of a game with $n$ players:

$$ \sum_{k=0}^\infty kp^n_k \approx \sum_{k=0}^K kp^n_k. $$

As I wrote above, there may be a closed form solution for the $p_k$, but for now, we can numerically evaluate them using R. I'm using $K=10,000$, so that $\sum_{k=0}^K p_k=1$ up to machine accuracy.

max_steps <- 10000
state_probabilities <- matrix(NA,nrow=max_steps+1,ncol=7,dimnames=list(0:max_steps,6:0))
state_probabilities[1,] <- c(1,0,0,0,0,0,0)
transition_matrix <- rbind(
    c(0,6,0,0,0,0,0),
    c(1,0,5,0,0,0,0),
    c(0,2,0,4,0,0,0),
    c(0,0,3,0,3,0,0),
    c(0,0,0,4,0,2,0),
    c(0,0,0,0,5,0,1),
    c(0,0,0,0,0,0,6))/6

for ( kk in 1:max_steps ) {
    state_probabilities[kk+1,] <- t(transition_matrix)%*%state_probabilities[kk,]
}

probs <- diff(state_probabilities[,7])
sum(probs)  # yields 1
sum(probs*seq_along(probs)) # yields 83.2

plot(probs[1:400],type="h",xlab="Number of steps",ylab="Probability",las=1)

Next, this is how we get the probabilities $p^4_k$ for $n=4$ players:

n_players <- 4

probs_minimum <- sapply(1:max_steps,
    function(kk)(1-sum(probs[1:(kk-1)]))^n_players-(1-sum(probs[1:kk]))^n_players)
head(probs_minimum)
plot(probs_minimum[1:400],type="h",xlab="Number of steps",ylab="Probability",
    las=1,main=paste(n_players,"players"))

Of course, four persons finish more quickly than a single person. For $n=4$, we get an expected value of

sum(probs_minimum*seq_along(probs_minimum))
[1] 25.44876

Finally, I like to confirm calculations like this using simulation.

n_sims <- 1e5
steps_minimum <- rep(NA,n_sims)
pb <- winProgressBar(max=n_sims)
for ( ii in 1:n_sims ) {
    setWinProgressBar(pb,ii,paste(ii,"of",n_sims))
    set.seed(ii)    # for reproducibility
    states <- matrix(FALSE,nrow=6,ncol=n_players)
    n_steps <- 0
    while ( TRUE ) {
        n_steps <- n_steps+1
        for ( jj in 1:n_players ) {
            roll <- sample(1:6,1)
            states[roll,jj] <- !states[roll,jj]
        }
        if ( any ( colSums(states) == 6 ) ) {
            steps_minimum[ii] <- n_steps
            break
        }
    }
}
close(pb)

The distribution of the numbers of steps needed in our $10^5$ simulated games matches the calculated $p^4_k$ rather well:

result <- structure(rep(0,length(probs_minimum)),.Names=seq_along(probs_minimum))
result[names(table(steps_minimum))] <- as.vector(table(steps_minimum))/n_sims
cbind(result,probs_minimum)[1:30,]
    result probs_minimum
1  0.00000    0.00000000
2  0.00000    0.00000000
3  0.00000    0.00000000
4  0.00000    0.00000000
5  0.00000    0.00000000
6  0.06063    0.06031414
7  0.00000    0.00000000
8  0.08072    0.07919228
9  0.00000    0.00000000
10 0.08037    0.08026479
11 0.00000    0.00000000
12 0.07382    0.07543464
13 0.00000    0.00000000
14 0.06826    0.06905406
15 0.00000    0.00000000
16 0.06409    0.06260212
17 0.00000    0.00000000
18 0.05668    0.05654555
19 0.00000    0.00000000
20 0.05180    0.05100393
21 0.00000    0.00000000
22 0.04570    0.04598101
23 0.00000    0.00000000
24 0.04078    0.04144437
25 0.00000    0.00000000
26 0.03749    0.03735245
27 0.00000    0.00000000
28 0.03241    0.03366354
29 0.00000    0.00000000
30 0.03026    0.03033861

Finally, the mean of the steps needed in the simulated games also matches the calculated expectation quite well:

mean(steps_minimum)
[1] 25.43862

Best Answer

Related Solutions

Solved – How to predict the results of a simple card game

Probability – Expected Rolls to Roll Every Number on a Dice an Odd Number of Times Explained

Related Question