R: monte carlo binomial gambling simulation

I am trying to figure out how to simulate a biased gambling problem using monte carlo simulations.

The problem is: Simulate two players tossing a coin; A and B. Player A has a 0.55 chance to win. Player B has a 0.45 chance. Each player starts with $3. If one player wins they take $1 from the other. The game ends when either one player has all the money OR 25 iterations have been played. I then want to plot the relative frequencies of players winning, then run this many time in order to get the estimates for player A and player B winning all the money.

What I am stuck on is getting the montecarlo simulations going and the calculating the probability of one player accumulating all the other player's money.

So far I can generate the data frame for one game and plot it.

Game <- c('Bethany', 'Algernon')   #outcomes in the game

#initialise an empty df
Games_data <- data.frame(Game = numeric(),
                        winner = character(),
                        Bethany_bank = numeric(),
                        Algernon_bank = numeric(),
                        Bethany_Freq = numeric(),
                        Algernon_Freq =  numeric()
)

#intialise variables
count <- 26
i <- 1
temp_Bethany_bank <- 3
temp_Algernon_bank <- 3

#populate the data frame until 25 games or someone wins
while(i < count) {
  temp_game <- i
  temp_winner <- sample(Game, prob =c(0.55, 0.45), size = 1)
  
  if(temp_winner == 'Bethany') {
    temp_Bethany_bank <- temp_Bethany_bank + 1
    temp_Algernon_bank <- temp_Algernon_bank - 1
  } else {
      temp_Bethany_bank <- temp_Bethany_bank - 1
      temp_Algernon_bank <- temp_Algernon_bank + 1}
  
  temp_Bethany_freq = 0.0
  temp_Algernon_freq = 0.0
  
  temp <- data.frame(Game = temp_game,
                     winner = temp_winner,
                     Bethany_bank = temp_Bethany_bank,
                     Algernon_bank = temp_Algernon_bank,
                     Bethany_Freq = temp_Bethany_freq,
                     Algernon_Freq = temp_Algernon_freq
                     )
  
  Games_data <- rbind(Games_data, temp)
  
  Games_data$Bethany_Freq <- cumsum(Games_data$winner == 'Bethany') / 1:nrow(Games_data)
  Games_data$Algernon_Freq <- cumsum(Games_data$winner == 'Algernon') / 1:nrow(Games_data)
  
  if(Games_data$Bethany_bank[i] <= 0 || Games_data$Algernon_bank[i] <= 0) {break} else {i <- i + 1}
}

#show the dataframe and the plot:
Games_data

ggplot(data = Games_data) +
  geom_point(aes(x = Game, y =  Bethany_Freq), color = 'coral', alpha = 0.8) + #Bethany's wins
  geom_point(aes(x = Game, y =  Algernon_Freq), color = 'steelblue', alpha = 0.8) + #Bethany's wins
  geom_line(aes(x = Game, y =  Bethany_Freq), color = 'coral') +
  geom_line(aes(x = Game, y =  Algernon_Freq), color = 'steelblue') +
  theme_classic() + 
  labs(title = "Relative frequency plots for Bethany vs Algernon over 25 games") 

How can I run this many times, like 100 or 1000, store the outputs in an object, plot all the trials and then calculate the probability of one player getting all the money?

Cheers in advance!

Wrap your code in a function, call it using replicate, then plot them all using, say, marrangeGrob from gridExtra.

mc_sim <- function() {
... your code

# Additional code to account for a "no result" game.
  if(i==count)
    Games_data <- rbind(Games_data, 
                        data.frame(Game = i,
                                   winner = "no result",
                                   Bethany_bank = temp_Bethany_bank,
                                   Algernon_bank = temp_Algernon_bank,
                                   Bethany_Freq = temp_Bethany_freq,
                                   Algernon_Freq = temp_Algernon_freq))

    Games_data
}

Run 9 simulations and save the plots to a list.

sim <- replicate(9, mc_sim(), simplify=FALSE)

sim_plots <- lapply(sim, \(x) ggplot(data = x) +
  geom_point(aes(x = Game, y =  Bethany_Freq), color = 'coral', alpha = 0.8) + #Bethany's wins
  geom_point(aes(x = Game, y =  Algernon_Freq), color = 'steelblue', alpha = 0.8) + #Bethany's wins
  geom_line(aes(x = Game, y =  Bethany_Freq), color = 'coral') +
  geom_line(aes(x = Game, y =  Algernon_Freq), color = 'steelblue') +
  theme_classic() 
)

Plot them

library(gridExtra)
marrangeGrob(sim_plots, nrow=3, ncol=3,
             top="Relative frequency plots for Bethany vs Algernon over 25 games")

The probabilities of each player winning can be calcualted using the relative frequencies.

prop.table(table(sapply(sim, \(x) x$winner[nrow(x)])))

 Algernon   Bethany 
0.4444444 0.5555556 

With more simulations, the probabilities converge to about 35% and 62% with 3% being "no result".

sim <- replicate(1000, mc_sim(), simplify=FALSE)
prop.table(table(sapply(sim, \(x) x$winner[nrow(x)])))
 #Algernon   Bethany no result 
 #   0.346     0.623     0.031