Create a function in R for simulating a coin flip game

Question

In a coin flip game, you flip a fair coin until the difference between the number of heads and number of tails is 3. You are paid $8 at the end, but you have to pay $1 for each flip of the coins. Use N =100000 simulations and find the expected amount you could win.

I was able to use the following code for 1 game but it breaks for N=100,000 as it is too long. How can I turn this into a function that will give me expected payoff for 100,000 games?

#For 1 game
# number of games to play
n_games <- 1 
bias <- 0.5
game_payoff <- c()

for (i in seq_len(n_games)) {
  
  cost <- 0
  flip_record <- c()
  payoff <- c()
  
  repeat{
    cost <- cost + 1
    flip <- rbinom(1, 1, prob = bias)
    flip_record <- c(flip_record, flip)
    # number of 0s/tails
    n_tails <- length(flip_record) - sum(flip_record) 
    # number of 1s/heads
    n_heads <- sum(flip_record) 
    
    if (abs(n_tails - n_heads) == 3) {
      # record game payoff
      game_payoff <- c(game_payoff, 8 - cost)
      # print game payoff
      print(paste0("single game payoff: ", 8 - cost)) 
      break
    }
  }
}

Answer 1

Here is a solution.
I have added an argument bankruptcy in order to avoid very long games.
In the example run below, only 10,000 games (n_games) are played.

game <- function(prob = 0.5, pay = 8L, bankruptcy = 1e6) {
  cost <- 0L
  difference <- 0L
  flip_record <- integer(bankruptcy)
  # loop until the difference between the number of heads and 
  # number of tails is 3 or the coin was flipped
  # more times than what the player can afford to flip 
  while(abs(difference) < 3 && cost < bankruptcy) {
    flip <- sample(c(-1, 1), 1, prob = c(prob, 1 - prob))
    cost <- cost + 1L
    difference <- difference + flip
    flip_record[cost] <- flip
  }
  win <- (abs(difference) == 3L)*pay
  flip_record <- flip_record[seq.int(cost)]
  list(win = win, cost = cost, flips = flip_record)
}

set.seed(2022)
n_games <- 1e4
prob <- 0.5

# the games are independent, this could use parallelization
system.time(
  games_played <- replicate(n_games, game(), simplify = FALSE)
)
#>    user  system elapsed 
#>   22.62   16.36   39.00

# extract the relevant values from the returned list
Wins <- sapply(games_played, `[[`, 'win')
Costs <- sapply(games_played, `[[`, 'cost')

# example differences sequences (game 1 and game 10)
Flips <- sapply(games_played, `[[`, 'flips')
cumsum(Flips[[1]])
#>  [1] -1 -2 -1 -2 -1 -2 -1  0  1  0  1  2  3

cumsum(Flips[[10]])
#> [1]  1  0 -1  0  1  2  3

# it pays more times than not
table(Wins - Costs > 0)
#> 
#> FALSE  TRUE 
#>  4258  5742

# the proportion of times the player turns a profit
mean(Wins - Costs > 0)
#> [1] 0.5742

# but when the players loose they loose more
# than the profit they get
mean(Wins - Costs) * n_games
#> [1] -9898

# the same number, different way of computing it
# and with formatted output
formatC(sum(Wins - Costs), big.mark = ",")
#> [1] "-9,898"

# the distribution of wins/losses is left skewed
# with a long tail => it doesn't pay to play
profit <- table(Wins - Costs)
barplot(profit)

^{Created on 2022-10-06 with reprex v2.0.2}

Answer 2

It will be more performant to vectorize over the games:

N <- 1e5
d <- rep(1L, N) # current flip difference for all N games; always 1 after the first flip
n <- 1L # the number of flips so far
flip <- c(-1L, 1L)
cost <- 0L # initialize the total cost of playing N games
while (length(d)) {
  d <- d + sample(flip, length(d), TRUE) # each flip adds or subtracts 1 from the difference
  idx <- which(abs(d) < 3L) # which games have not ended?
  cost <- cost + (n <- n + 1L)*(length(d) - length(idx)) # add the cost from the ended games
  d <- d[idx] # remove the ended games from the simulation
}

# expected payoff
8 - cost/N
#> [1] -0.99824