Discussion_8

11. Assume that you are playing craps with dice that are loaded in the following
way: faces two, three, four, and five all come up with the same probability
(1/6) + r. Faces one and six come up with probability (1/6) − 2r, with 0 <
r < .02. Write a computer program to find the probability of winning at craps
with these dice, and using your program find which values of r make craps a
favorable game for the player with these dice.

There are two functions: simulate_craps that simulates one game of craps using the loaded dice, and prob_win_craps that calculates the probability of winning at craps with the loaded dice for a given value of r.

To find which values of r make craps a favorable game for the player, we can loop over a range of r values and calculate the probability of winning for each value.

We are getting negative probabilities error because the probability distribution for the loaded dice includes negative probabilities for faces one and six when r > 0.01.

To avoid this error, we can check for negative probabilities and set them to zero before simulating the game of craps.

# Define function to simulate one game of craps
simulate_craps <- function(p_one, p_six, p_others) {
  # Roll the dice
  roll1 <- sample(1:6, 1, prob = c(p_one, p_others, p_others, p_others, p_six, p_others))
  roll2 <- sample(1:6, 1, prob = c(p_one, p_others, p_others, p_others, p_six, p_others))
  total <- roll1 + roll2
  
  # Determine outcome of the game
  if (total %in% c(7, 11)) {
    outcome <- "win"
  } else if (total %in% c(2, 3, 12)) {
    outcome <- "lose"
  } else {
    point <- total
    repeat {
      roll1 <- sample(1:6, 1, prob = c(p_one, p_others, p_others, p_others, p_six, p_others))
      roll2 <- sample(1:6, 1, prob = c(p_one, p_others, p_others, p_others, p_six, p_others))
      total <- roll1 + roll2
      if (total == point) {
        outcome <- "win"
        break
      } else if (total == 7) {
        outcome <- "lose"
        break
      }
    }
  }
  
  # Return outcome of the game
  return(outcome)
}

# Define function to calculate probability of winning at craps with loaded dice
prob_win_craps <- function(r) {
  # Define probabilities of rolling each face of the loaded dice
  p_one <- max(1/6 - 2*r, 0)
  p_six <- max(1/6 - 2*r, 0)
  p_others <- (1/6) + r
  
  # Simulate many games of craps and count the number of wins
  num_wins <- 0
  num_games <- 100000
  for (i in 1:num_games) {
    outcome <- simulate_craps(p_one, p_six, p_others)
    if (outcome == "win") {
      num_wins <- num_wins + 1
    }
  }
  
  # Calculate probability of winning
  prob_win <- num_wins / num_games
  
  # Return probability of winning
  return(prob_win)
}

# Loop over values of r from 0 to 0.02 in increments of 0.001
for (r in seq(0, 0.02, 0.001)) {
  # Calculate probability of winning with loaded dice for current value of r
  prob_win <- prob_win_craps(r)
  
  # Print value of r and probability of winning
  cat("r =", r, ", probability of winning =", prob_win, "\n")
}

## r = 0 , probability of winning = 0.49423 
## r = 0.001 , probability of winning = 0.49516 
## r = 0.002 , probability of winning = 0.49277 
## r = 0.003 , probability of winning = 0.49271 
## r = 0.004 , probability of winning = 0.49525 
## r = 0.005 , probability of winning = 0.49389 
## r = 0.006 , probability of winning = 0.49402 
## r = 0.007 , probability of winning = 0.49437 
## r = 0.008 , probability of winning = 0.49482 
## r = 0.009 , probability of winning = 0.49536 
## r = 0.01 , probability of winning = 0.49829 
## r = 0.011 , probability of winning = 0.49646 
## r = 0.012 , probability of winning = 0.49404 
## r = 0.013 , probability of winning = 0.49835 
## r = 0.014 , probability of winning = 0.49871 
## r = 0.015 , probability of winning = 0.49888 
## r = 0.016 , probability of winning = 0.50129 
## r = 0.017 , probability of winning = 0.49917 
## r = 0.018 , probability of winning = 0.50223 
## r = 0.019 , probability of winning = 0.49983 
## r = 0.02 , probability of winning = 0.50125

With the max() function added to the definition of p_one and p_six to set any negative probabilities to zero. The program also includes a loop over values of r from 0 to 0.02 in increments of 0.001.