Simulation in R for probability puzzle

Question

Let two integers such that 1<a<7 and 2<b<5.We play the following game: We start with a euros and we roll a dice.If the result of the dice is at least b we earn 1 euro and we lose 1 euro otherwise. The game ends if we reach 8 euros (winners) or in 0 (losers). We prefer to be ?

A) Rich: we start with a equals 6 euros , but unlucky because b equals 5. B) Mediocre: we start with a equals 4 and b equals 4 C) Lucky: we start with a equals 2 but b equals 3.

How can I do it in R using simulation in order to decide If I want to be rich ,mediocre or lucky?

My effort


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)

Your code confuses the bet with the die roll. In your description, you always win 1 euro or lose 1 euro, but in your code you always win `bet` euros, where `bet` is a die roll. You also ignore `p`, which you need to use. I would rename your variable `roll <- sample(1:6, size = 1)`, and then rewrite the rest of your function. — Gregor Thomas, Sep 24 '21 at 19:17
Also, `if (stake == 0) return(1) else return(0)` seems very confusing the way you switch things. I would suggest `if(stake == 0) return(0) else return(1)`. — Gregor Thomas, Sep 24 '21 at 19:20
(Side note). This is the gambler's ruin (https://en.wikipedia.org/wiki/Gambler%27s_ruin). I would suggest modeling it with a Markov Chain (https://en.wikipedia.org/wiki/Markov_chain). In addition, the die can be replaced with a coin and still generate the same experiment. The probability of the coin will be sum(p(1:b)) (p is a vector of die probabilities). You will then have an analytical way to study this experiment. See here for details: http://www.columbia.edu/~ks20/stochastic-I/stochastic-I-GRP.pdf — Baraliuh, Sep 24 '21 at 22:49

score 2 · Accepted Answer · answered Sep 24 '21 at 22:37

2

We could define a function f() that uses recursion.

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

answered Sep 24 '21 at 22:37

jay.sf

60,139
8
53
110

is this a R function? I have a problem run it in my computer – Homer Jay Simpson Sep 25 '21 at 09:52
@HungryHomer Just abbreviated function notation, available since R4.1. Replace `\(a, b)` and `\(x)` with `function(a, b)` and `function(x)`. Or even better: update your R. – jay.sf Sep 25 '21 at 10:09

Christopher Belanger · Answer 2 · 2021-09-24T19:58:28.100

Here's an extremely over-engineered solution!

Your loop logic looks fine to me, the main issue is that you were originally rolling a die and just adding that many euros to your stake, and your function wasn't returning anything.

The key line below is euros <- euros + (sample(1:6, 1) >= roll_to_win)*2 - 1 . This does the following:

Gets a random number between 1 and 6,
Checks to see if it's bigger than or equal to roll_to_win, which corresponds to b,
Converts a TRUE result to 1 and a FALSE result to -1 (with trickery! TRUE and FALSE behave like 1 and 0, so if you multiply TRUE by 2 and subtract 1 you get 1, and you get -1 for FALSE.

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))

And here's what it looks like for me:

# 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

No clue if that's right or not though :)

This is correct if roll_to_win must exceed the dice roll in order to win. If you can win by matching the dice roll, then change `sample(1:6, 1) > roll_to_win` to use `>=` — Jordo82, Sep 24 '21 at 19:48
You're right! The initial question said "at least b" and I missed that. I've updated the answer. — Christopher Belanger, Sep 24 '21 at 19:58
It’s correct.I edited also my effort.I calculated the probability of ruin so you calculated the opposite.The results are symmetric the same if you see.Thanks for your help though. — Homer Jay Simpson, Sep 24 '21 at 20:08
Therefore the correct answer is C.I want to to be lucky.Not rich not mediocre. — Homer Jay Simpson, Sep 24 '21 at 20:10

Simulation in R for probability puzzle

2 Answers2