如何实现蒙特卡洛模拟偏态赌博并计算获胜概率？

步骤1：将单局模拟封装为函数 #

把单局模拟的逻辑打包成可重复调用的函数，既方便批量执行，也能同时返回每局的详细过程数据和最终胜负结果。

library(ggplot2)
library(dplyr)

# 单局模拟函数
simulate_single_game <- function() {
  Game <- c('Bethany', 'Algernon')
  
  # 初始化本局数据框
  Games_data <- data.frame(
    Game = numeric(),
    winner = character(),
    Bethany_bank = numeric(),
    Algernon_bank = numeric(),
    Bethany_Freq = numeric(),
    Algernon_Freq = numeric()
  )
  
  count <- 26
  i <- 1
  temp_Bethany_bank <- 3
  temp_Algernon_bank <- 3
  
  # 模拟对局直到终止条件触发
  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 <- data.frame(
      Game = temp_game,
      winner = temp_winner,
      Bethany_bank = temp_Bethany_bank,
      Algernon_bank = temp_Algernon_bank,
      Bethany_Freq = 0.0,
      Algernon_Freq = 0.0
    )
    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)
    
    # 检查终止条件：一方破产或完成25局
    if(temp_Bethany_bank <= 0 || temp_Algernon_bank <= 0) {
      break
    } else {
      i <- i + 1
    }
  }
  
  # 判定本局最终结果
  final_result <- case_when(
    temp_Bethany_bank <= 0 ~ "Algernon_win",
    temp_Algernon_bank <= 0 ~ "Bethany_win",
    TRUE ~ "Draw"
  )
  
  # 返回本局过程数据和最终结果
  list(game_data = Games_data, final_result = final_result)
}

步骤2：执行批量蒙特卡洛模拟 #

设定模拟次数（比如1000次），通过循环批量调用单局函数，收集所有模拟的过程数据和最终结果。

# 设置模拟次数
num_simulations <- 1000

# 执行批量模拟并添加模拟ID标记
simulation_results <- lapply(1:num_simulations, function(sim_id) {
  result <- simulate_single_game()
  result$game_data$simulation_id <- sim_id
  result
})

# 合并所有模拟的过程数据为一个大数据框
all_game_data <- do.call(rbind, lapply(simulation_results, function(x) x$game_data))

# 整理所有模拟的最终结果
final_results <- data.frame(
  simulation_id = 1:num_simulations,
  result = sapply(simulation_results, function(x) x$final_result)
)

步骤3：绘制所有试验的相对频率图 #

提供两种可视化方式，根据需求选择：

方式1：所有模拟叠加图（直观展示整体趋势） #

把所有模拟的频率变化线画在同一张图上，用低透明度区分单条线，能清晰看到整体的频率走向：

ggplot(all_game_data) +
  geom_line(aes(x = Game, y = Bethany_Freq, group = simulation_id), color = 'coral', alpha = 0.1) +
  geom_line(aes(x = Game, y = Algernon_Freq, group = simulation_id), color = 'steelblue', alpha = 0.1) +
  theme_classic() +
  labs(
    title = "所有模拟的相对频率变化",
    x = "对局次数",
    y = "获胜相对频率"
  )

方式2：分面展示单局情况（适合小次数模拟） #

如果模拟次数较少（比如100次），可以用分面逐个展示每局的频率变化：

ggplot(all_game_data) +
  geom_line(aes(x = Game, y = Bethany_Freq), color = 'coral') +
  geom_line(aes(x = Game, y = Algernon_Freq), color = 'steelblue') +
  facet_wrap(~simulation_id, ncol = 10) +
  theme_classic() +
  theme(axis.text = element_text(size = 6)) +
  labs(
    title = "每局模拟的相对频率变化",
    x = "对局次数",
    y = "获胜相对频率"
  )

步骤4：计算双方赢走全部资金的概率 #

通过统计所有模拟的最终结果，计算各类情况的发生概率：

# 统计各类结果的次数
result_counts <- table(final_results$result)

# 计算概率并保留4位小数
probabilities <- prop.table(result_counts)

cat("Bethany赢走全部资金的概率：", round(probabilities["Bethany_win"], 4), "\n")
cat("Algernon赢走全部资金的概率：", round(probabilities["Algernon_win"], 4), "\n")
cat("25局后平局的概率：", round(probabilities["Draw"], 4), "\n")

内容的提问来源于stack exchange，提问作者Lavacave