Бесконечный цикл в игре с подбрасыванием монет

Вы пересчитываете w в каждой итерации цикла на основе текущего значения heads и tails. Но эти значения всегда будут 1 и 0 (или 0 и 1). Таким образом, w всегда либо -1, либо 1, и никогда не имеет другого значения.

Другая ошибка в вашем коде заключается в том, что вы останавливаетесь только тогда, когда впереди 4 головы. Но по правилам игра должна останавливаться и тогда, когда впереди 4 решки: имеет значение только абсолютная разница.

Логика вашего кода может быть исправлена, но будет работать гораздо более простая логика (обратите внимание, что в следующем коде используются понятные имена переменных, что делает результирующий код более читабельным):

h = function () {
    sides = c('H', 'T')
    diff = 0L
    cost = 0L
    repeat {
        cost = cost + 1L
        flip = sample(sides, 1L)
        if (flip == 'H') diff = diff + 1L
        else diff = diff - 1L
        if (abs(diff) == 4L) return(cost)
    }
}

Вы можете упростить это еще больше, потому что метки сторон монеты на самом деле не имеют значения. Все, что вам нужно, это подбрасывание монеты, которое возвращает один из двух результатов.

Мы можем реализовать это как отдельную функцию. Возвращаемое значение функции не очень важно, пока у нас есть фиксированное соглашение: оно может быть в c('H', 'T'), или c(FALSE, TRUE), или c(0L, 1L) и т. д. Для наших целей было бы удобно возвращать либо -1L, либо 1L, чтобы наша функция h могла напрямую добавлять это значение к diff:

coin_toss = function () {
    sample(c(-1L, 1L), 1L)
}

Но есть и другой способ получения подбрасывания монеты: испытание Бернулли размера 1. И использование испытания Бернулли имеет хорошее свойство: мы можем тривиально расширить нашу функцию, чтобы разрешить несправедливые (предвзятые) подбрасывания монеты. Итак, вот та же функция, но с необязательным параметром bias (по умолчанию монетка подбрасывается честно):

coin_toss = function (bias = 0.5) {
    rbinom(1L, 1L, prob = bias) * 2L - 1L
}

(rbinom(…) возвращает либо 0L, либо 1L. Чтобы преобразовать область значений в c(-1L, 1L), мы умножаем на 2 и вычитаем 1.)

Теперь давайте изменим h, чтобы использовать эту функцию:

h = function (bias = 0.5) {
    cost = 0L
    diff = 0L
    repeat {
        cost = cost + 1L
        diff = diff + coin_toss(bias)
        if (abs(diff) == 4L) return(cost)
    }
}

coin_toss() равно 0 или 1, но, в зависимости от его значения, мы либо

Я хотел бы ответить на ваши вопросы, как в части а), так и в части б). Я буду использовать свои коды, чтобы сэкономить время.

Это крутая игра, в которой программная симуляция может оказаться очень полезной. Суть игры - это «бесконечный цикл», который в конечном итоге заканчивается, когда абсолютная разница количества орлов и решек равна 4. Затем записывается выигрыш. Как упомянул Конрад Рудольф, это игра типа Бернулли. Игра моделируется с помощью кода ниже:

n_games <- 1000 # number of games to play
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)

    n_tails <- length(flip_record) - sum(flip_record) # number of 0s/tails
    n_heads <- sum(flip_record) # number of 1s/heads
    
    if (abs(n_tails - n_heads) == 4) {
      game_payoff <- c(game_payoff, 10 - cost) # record game payoff
      print(paste0("single game payoff: ", 10 - cost)) # print game payoff
      break
    }
  }
}

При большом количестве прогонов, т.е. еще один цикл над этим циклом, мы узнаем, что ожидаемое значение очень близко к -6. Таким образом, игра имеет отрицательное математическое ожидание. Из этого кода следует:

library(ggplot2)
seed <- 122334

# simulation
n_runs <- 100
n_games <- 10000
bias <- 0.5

game_payoff <- c()
expected_value_record <- c()

for (j in seq_len(n_runs)) {
  
  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)
      # print(flip_record)
      
      n_tails <- length(flip_record) - sum(flip_record) # number of 0s/tails
      n_heads <- sum(flip_record) # number of 1s/heads
      
      if (abs(n_tails - n_heads) == 4) {
        game_payoff <- c(game_payoff, 10 - cost) # record game payoff
        print(paste0("single game payoff: ", 10 - cost))
        break
      }
    }
  }
  expected_value_record <- c(expected_value_record, mean(game_payoff))
  game_payoff <- c()
}

# plot expected value
expected_value_record <- cbind.data.frame("run" = seq_len(length(expected_value_record)), expected_value_record)

ggplot(data = expected_value_record) +
  geom_line(aes(x = run, y = expected_value_record)) +
  scale_x_continuous(breaks = c(seq(1, max(expected_value_record$run), by = 3), max(expected_value_record$run))) +
  labs(
    title = "Coin flip experiment: expected value in each run. ", 
    caption = paste0("Number of runs: ", n_runs, ". ", "Number of games in each run: ", n_games, "."), 
    x = "Run", 
    y = "Expected value") +
  geom_hline(yintercept = mean(expected_value_record$expected_value_record), size = 1.4, color = "red") +
  annotate(
    geom = "text",
    x = 0.85 * n_runs,
    y = max(expected_value_record$expected_value_record),
    label = paste0("Mean across runs: ", mean(expected_value_record$expected_value_record)),
    color = "red") +
  theme(plot.title = element_text(hjust = 0.5), plot.caption = element_text(hjust = 0.5))

Графика:

Давайте теперь посмотрим на часть б) вопроса с другой симуляцией. Цикл был завернут в функцию, которую с помощью sapply мы прогоняем последовательность вероятностей:

library(ggplot2)
seed <- 122334

# simulation function
coin_game <- function(n_runs, n_games, bias = 0.5){
  game_payoff <- c()
  expected_value_record <- c()
  
  for (j in seq_len(n_runs)) {
    
    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)
        # print(flip_record)
        
        n_tails <- length(flip_record) - sum(flip_record) # number of 0s/tails
        n_heads <- sum(flip_record) # number of 1s/heads
        
        if (abs(n_tails - n_heads) == 4) {
          game_payoff <- c(game_payoff, 10 - cost) # record game payoff
          break
        }
      }
    }
    expected_value_record <- c(expected_value_record, mean(game_payoff))
    game_payoff <- c()
  }
  return(expected_value_record)
}

# run coin_game() on a vector of probabilities - introduce bias to find fair game conditions
n_runs = 1
n_games = 1000
expected_value_record <- sapply(seq(0.01, 0.99, by = 0.01), coin_game, n_runs = n_runs, n_games = n_games)

# plot expected value
expected_value_record <- cbind.data.frame("run" = seq_len(length(expected_value_record)), "bias" = c(seq(0.01, 0.99, by = 0.01)), expected_value_record)

ggplot(data = expected_value_record) +
  geom_line(aes(x = bias, y = expected_value_record)) +
  scale_x_continuous(breaks = c(seq(min(expected_value_record$bias), max(expected_value_record$bias), by = 0.1), max(expected_value_record$bias))) +
  scale_y_continuous(breaks = round(c(0, seq(min(expected_value_record$expected_value_record), max(expected_value_record$expected_value_record), length.out = 10)), digits = 4)) +
  labs(
    title = "Coin flip experiment: expected value for each probability level", 
    caption = paste0("Number of runs per probability level: ", n_runs, ". ", "Number of games in each run: ", n_games, "."), 
    x = "Probability of success in Bernoulli trial", 
    y = "Expected value") +
  geom_hline(yintercept = 0, size = 1.4, color = "red") +
  geom_text(aes(x = 0.1, y = 0, label = "Fair game", hjust = 1, vjust = -1), size = 4, color = "red") +
  theme(plot.title = element_text(hjust = 0.5), plot.caption = element_text(hjust = 0.5))

Графика:

Изучение фрейма данных expect_value_record показывает, что игра честна, когда значения вероятности находятся в диапазонах: 0,32–0,33 или 0,68–0,69.

Последний код легко изменить, чтобы выжать из него более надежные числа.