R 中的概率谜题模拟

Simulation in R for probability puzzle

让满足 1

A) Rich:我们从a开始等于6欧元,但不走运,因为b等于5。 B) 平庸:我们从 a 等于 4 和 b 等于 4 开始 C) 幸运:我们从 a 等于 2 但 b 等于 3.

我如何在 R 中使用模拟来决定我是想变得富有、平庸还是幸运?

我的努力


gamble <- function(a,n,p) {
  stake <- a
  while (stake > 0 & stake < n) {
    bet <- sample(c(-1,1),1,prob=c(1-p,p))
    stake <- stake + bet
  }
  if (stake == 0) return(1) else return(0)
}   

a <- 6
n <- 8
p <- 1/3
trials <- 100000
simlist <- replicate(trials, gamble(a, n, p))
mean(simlist) # Estimate of probability that gambler is ruined

a <- 4
n <- 8
p <- 1/2
trials <- 100000
simlist <- replicate(trials, gamble(a, n, p))
mean(simlist) 


a <- 2
n <- 8
p <- 2/3
trials <- 100000
simlist <- replicate(trials, gamble(a, n, p))
mean(simlist) 


这是一个设计过度的解决方案!

我觉得你的循环逻辑没问题,主要问题是你最初是在掷骰子,只是把那么多欧元加到你的赌注上,而你的函数没有返回任何东西。

下面的关键行是euros <- euros + (sample(1:6, 1) >= roll_to_win)*2 - 1 。这会执行以下操作:

  • 获取1到6之间的随机数,
  • 检查是否大于等于roll_to_win,对应b
  • TRUE 结果转换为 1,将 FALSE 结果转换为 -1(有技巧!TRUEFALSE 表现得像 1 和 0,所以如果你TRUE 乘以 2 减 1 得到 1,FALSE.
  • 得到 -1
library(tidyverse)

# set up our gambling function. We've hard-coded that it's 8 euros to win.
gamble <- function(start_euros,
                   roll_to_win){
  
  euros <- start_euros
  win_euros <- 8
  
  while(euros > 0 & euros < win_euros){
    # for debugging, print euros to console
    #message(euros)
    
    # roll a 6-sided die, see if the result is
    # bigger than the value we need to win, and
    # add or subtract 1 from euros based on result
    euros <- euros + (sample(1:6, 1) >= roll_to_win)*2 - 1    
    
  }
  
  # return 0 for loss, 1 for win
  return (euros / win_euros)
}

# set up the number of reps we'll do
reps <- 10000

# build a nested tibble with our values for a and b and generate trials
results <- tribble(~a,~ b, 
                   6, 5,
                   4, 4,
                   2, 3) %>%
  mutate(trials = purrr::map2(a,b, function(x, y) replicate(reps, gamble(x,y))))

# get observed frequency by taking mean 
results %>%
  mutate(prob = purrr::map_dbl(trials, mean))

这对我来说是这样的:

# A tibble: 3 x 4
      a     b trials          prob
  <dbl> <dbl> <list>         <dbl>
1     6     5 <dbl [10,000]> 0.240
2     4     4 <dbl [10,000]> 0.498
3     2     3 <dbl [10,000]> 0.751

不知道这是否正确:)

我们可以定义一个使用递归的函数f()

f <- \(a, b) {
  s <- sample(6, 1)
  l <- length(a)
  a <- c(a, a[l] + ifelse(s >= b, 1, -1))
  if (a[l + 1] %in% c(0, 8)) return(a)
  else return(f(a, b))
}

R <- 1e5
set.seed(42)
r0 <- replicate(R, Map(f, c(6, 4, 2), c(5, 4, 3)))
rowMeans(matrix(rapply(r0, \(x) x[length(x)] == 8), 3))
# [1] 0.24668 0.49885 0.75435