Untitled

#Libraries:

# Load the MASS package for the fitdistr function
library(MASS)

http://cs.gettysburg.edu/projects/pig/pigCompare.html

— This is a link to things that could help

Bank

This page is to practice creating the bank game, and then trying new strategies to win? Should you bank low? Are you really going to win by being patient and banking really low? Should you just follow everyone else? What about how often do big winners win? What is a consistent strategy? What is the standard deviation of a round? How often do we roll 7’s?

Distribution of Dice?

First we need a virtual die to be rolled:

# Function to simulate rolling multiple dice with a specified number of sides
roll_dice <- function(num_dice = 1, num_sides = 6) {
  results <- sample(1:num_sides, num_dice, replace = TRUE)
  return(results)
}

# Test the function
num_dice <- 10
num_sides <- 6
rolled_values <- roll_dice(num_dice, num_sides)
print(rolled_values)

 [1] 6 1 1 6 5 2 2 2 3 3

Now that we have a virtual die, let’s graph 1000 dice rolls. Can you already picture the distribution?

#set values and run 1000 dice
num_dice <- 1000
num_sides <- 6
rolled_values <- roll_dice(num_dice, num_sides)


# Plot a histogram of the rolled values
hist(rolled_values, breaks = seq(0.5, max(rolled_values) + 0.5, by = 1), 
     main = "1000 Dice Rolls", xlab = "Dice Roll Result", ylab = "Frequency")

Uniform distribution transformed?! Distribution of Two Dice Visualized

That was fun but now let’s look at the sum of two dice. Do you already know what this distribution is?

num_dice <- 1000
num_sides <- 6
rolled_values <- roll_dice(num_dice, num_sides)
rolled_second <- roll_dice(num_dice, num_sides)




roll_sum <- rolled_values + rolled_second

table(roll_sum)

roll_sum
  2   3   4   5   6   7   8   9  10  11  12 
 24  69  95 110 135 183 123 104  81  48  28 

# Plot a histogram of the summed values
hist(roll_sum, breaks = seq(1, max(roll_sum) + 0.5, by = 1), 
     main = "Histogram of Sum of Two Dice Rolls", xlab = "Sum of Two Dice Rolls", ylab = "Frequency")

Okay, This graph is interesting to me… because even though it takes a full second for the computer to compute this, it still seems to still have some skiwompus— ness. Isn’t it weird that even with 100 samples, we are still struggling to see the perfect distribution? I guess we just need more samples? We could improve the way our code works, but later we need to know if it was doubles or not, and I dont really want to rewrite it just for that lol… just deal with a slow program for now… sorry.

Did you lose the game?

Now it’s time to write a program to play the game! We want to have the game running, and then to be able to add players… tbh, I didn’t really think about how this would work in R. Does R even have classes? Would using classes make the most sense?

# Initialize an empty data frame to store the results
results <- data.frame(round = integer(), total = integer(), doubles = logical())

# Function to simulate rolling two dice and adding their sum to the total
roll_dice_and_sum <- function(num_sides = 6) {
  dice1 <- sample(1:num_sides, 1, replace = TRUE)
  dice2 <- sample(1:num_sides, 1, replace = TRUE)
  total <- dice1 + dice2
  doubles <- dice1 == dice2
  return(list(total = total, doubles = doubles))
}

run_game <- function(results = NULL, doubles = FALSE, safe_3 = FALSE, num_rounds = 10, num_sided_dice = 6){
  # Initialize an empty data frame to store the results
  if(is.null(results))
  {
    results <- data.frame(round = integer(), num_rolls = integer(), round_total = integer(), stringsAsFactors = FALSE)
  }
  

  # Initialize an empty data frame to store the results

  for (round_number in 1:num_rounds) {
  # Initialize variables for the current round
  round_total <- 0
  num_rolls <- 0
  
  # Simulate dice rolls until a 7 is rolled
  while (TRUE) {
    # Simulate a single dice roll
    dice_result <- roll_dice_and_sum(num_sided_dice)
    round_value <- dice_result$total
    
    #check for doubles and Update round_total
    if (dice_result$doubles && doubles)
    {
      round_total <- round_total * 2
    }
    else{
    
      round_total <- round_total + round_value
    }
    
    # Increment the number of rolls
    num_rolls <- num_rolls + 1
    
    if (num_rolls <= 3 && dice_result$total == 7 && safe_3) 
    {
    #update the value of this roll 
    dice_result$value = 50 
  }
  
    
    # Check if a 7 was rolled, ending the round
    if (round_value == 7) {
      break
    }
  }
  
  # Store summary information for the round
  results <- rbind(results, data.frame(round = round_number, num_rolls = num_rolls, round_total = round_total))
  
  # Increment the round number
  round_number <- round_number + 1
  }
  return(results)
}


hello_world <- run_game()

There’s the function… now let’s make some data!

game_results <- run_game(num_rounds = 1000)

Let’s run some statistics on our first test :)

# Compute summary statistics for the total value of each round
round_summary <- summary(game_results$round_total)
print(round_summary)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   7.00   15.00   28.00   40.44   55.00  328.00 

# Fit an exponential distribution to the total value of each round
fit <- fitdistr(game_results$round_total, "exponential")

# Print the estimated parameters of the fitted exponential distribution
print(fit$estimate)

      rate 
0.02472922 

# Plot the histogram of the total value of each round
hist(game_results$round_total, breaks = 20, freq = FALSE, main = "Histogram of Round Total", xlab = "Total Value", ylab = "Density")

# Add the fitted exponential curve to the plot
curve(dexp(x, rate = fit$estimate["rate"]), add = TRUE, col = "red", lwd = 2)

It’s interesting that it is ALMOST an exponential… but the first values feel much too high? What is up with that? Are the doubles really throwing the spread so terribly? Let’s move on to looking at the score with doubles

# Compute mean, standard deviation, and quartiles for the total value of each round
mean_total <- mean(game_results$round_total)
sd_total <- sd(game_results$round_total)
quantiles_total <- quantile(game_results$round_total, probs = c(0.25, 0.5, 0.75))

print("Mean:")

[1] "Mean:"

print(mean_total)

[1] 40.438

print("Standard Deviation:")

[1] "Standard Deviation:"

print(sd_total)

[1] 37.65033

print("Quartiles:")

[1] "Quartiles:"

print(quantiles_total)

25% 50% 75% 
 15  28  55 

With doubles:

double_game <- run_game(doubles = TRUE, safe_3 = TRUE, num_rounds = 1000)


# Generate x values representing possible values of the geometric distribution
x_values <- 1:max(double_game$num_rolls)

# Plot the histogram of the total value of each round
hist(double_game$num_rolls, freq = FALSE, main = "Histogram of Round Total",
     xlab = "Total Value", ylab = "Density", breaks = 15)

# Fit a geometric distribution to the total value of each round
fit_geom <- fitdistr(double_game$num_rolls, "geometric")

# Calculate the PMF of the geometric distribution
pmf_values <- dgeom(x_values, prob = fit_geom$estimate["prob"])

# Plot the PMF of the geometric distribution
points(x_values, pmf_values, type = "h", col = "blue", lwd = 2)