Amoy Game Copyrights @ Zhen (Peter) Ye, May 2024 ###################################################################################################################################### # Here is a game I call the Amoy game, from my hometown Xiamen. It's been played for hundred of years during the Mid-Autumn Festival. # The game can be used for coding combinatorics and involves the following three steps. # (1) Calculate the probability of success with the number of 4s (once, twice, thrice, fourth, fifth and sixth times) appearing in dice throw. # (2) Simulate 1 million throw and tabulate the result for comparison. # (3) Graphing the tables and exporting the graph to R markdown. # There are prizes to be won with the game for the people who calculate the result correctly using R. #################################################################################################### # Prizes are you got a certificate and you may leave the class early. # Post-doc: six 4 in six-sided dice # PhD: five 4 in six-sided dice # Master: three 4 in six-sided dice # Bachelor: two 4 in six-sided dice # High school: one 4 in six-sided dice ################################################################################################# # Here is one start-up example to calculate the probability of 4s appearing once in dice throw k <- 1 # number of "4s" desired n <- 6 # total number of three-sided dice p <- 1/6 # probability of rolling a "4" on one die # Calculate the probability of rolling exactly one "4s" in six three-sided dice prob_1_4s_six_sided_dice <- dbinom(k, n, p) # Print the probability prob_1_4s_six_sided_dice # Alternatively, we use choose # Define variables n <- 6 k <- 1 p <- 1/6 # Using choose to calculate probabilities prob_1_4s <- choose(n, k) * p^k * (1-p)^(n-k) prob_1_4s ################################################################################################# # # # For your class of students, you get them to use dbinom and choose to find out the probability.# # # ################################################################################################# # Now use shorter codes in block to create a comparison table # Calculate probabilities using dbinom prob_6_4s_six_sided_dice <- dbinom(6, 6, 1/6) prob_5_4s_six_sided_dice <- dbinom(5, 6, 1/6) prob_4_4s_six_sided_dice <- dbinom(4, 6, 1/6) prob_3_4s_six_sided_dice <- dbinom(3, 6, 1/6) prob_2_4s_six_sided_dice <- dbinom(2, 6, 1/6) prob_1_4s_six_sided_dice <- dbinom(1, 6, 1/6) # Calculate probabilities using choose and manual formula prob_6_4s <- choose(6, 6) * (1/6)^6 * (5/6)^0 prob_5_4s <- choose(6, 5) * (1/6)^5 * (5/6)^1 prob_4_4s <- choose(6, 4) * (1/6)^4 * (5/6)^2 prob_3_4s <- choose(6, 3) * (1/6)^3 * (5/6)^3 prob_2_4s <- choose(6, 2) * (1/6)^2 * (5/6)^4 prob_1_4s <- choose(6, 1) * (1/6)^1 * (5/6)^5 # Create a table to compare the results comparison_table <- matrix(c( prob_6_4s_six_sided_dice, prob_5_4s_six_sided_dice, prob_4_4s_six_sided_dice, prob_3_4s_six_sided_dice, prob_2_4s_six_sided_dice, prob_1_4s_six_sided_dice, prob_6_4s, prob_5_4s, prob_4_4s, prob_3_4s, prob_2_4s, prob_1_4s ), nrow = 2, byrow = TRUE) # Assign column names colnames(comparison_table) <- c("6_4s", "5_4s", "4_4s", "3_4s", "2_4s", "1_4s") # Assign row names rownames(comparison_table) <- c("dbinom", "choose") # Print the table print(comparison_table) #################################################################################### # R code for simulation and tabulating the result # #################################################################################### # Simulate dice roll # Create a data frame # Tabulate the result # Compare the result of the simulation with calculated result #################################################################################### # The code is shortened to the below. # Simulation Function simulate_rolls <- function(n_dice, n_sims, successes_needed) { successes <- 0 for (i in 1:n_sims) { rolls <- sample(1:6, n_dice, replace = TRUE) number_of_fours <- sum(rolls == 4) if (number_of_fours == successes_needed) { successes <- successes + 1 } } return(successes) } #################################################################################### # Run Simulations and Tabulate Results n_sims <- 1000000 n_dice <- 6 successes_needed <- 1:6 # Create an empty data frame to store results results_df <- data.frame( 'Number of 4s' = successes_needed, 'Probability' = numeric(length(successes_needed)) ) # Perform simulations for each number of "4s" needed for (i in 1:length(successes_needed)) { result <- simulate_rolls(n_dice, n_sims, successes_needed[i]) probability_from_simulation <- result / n_sims results_df[i, 'Number of Successes'] <- result results_df[i, 'Probability'] <- probability_from_simulation } # Print the results in a tabulated format print(results_df) ########################################################################################### #Graphing the result dataframe library(ggplot2) # Use ggplot2 to plot the results names(results_df) library(ggplot2) # Updated ggplot code to display full probability values on labels ggplot(data = results_df, aes(x = factor(Number.of.4s), y = Probability)) + geom_bar(stat = "identity", fill = "skyblue") + # Lighter blue for bars geom_label(aes(label = as.character(Probability)), # Display full value of Probability position = position_stack(vjust = 0.5), fill = "orange", # Background color of the label remains black color = "steelblue", # Text color remains white size = 3.5) + # Text size remains the same labs(title = "Probability of Rolling Certain Numbers of '4s'", x = "Number of '4s'", y = "Probability", caption = "Data based on 1 million simulations of rolling 6 dice.") + theme_minimal() ##################################################################################### ##### # # # Copyrights @ Zhen (Peter) Ye, May 2024 #####