# COSC/MATH 201 - Randomness in Practice (Modules 9.2 and 9.5)
# Estimate the chance of a shared birthday by simulation, check it against the exact answer,
# and see how the chance grows with class size.
set.seed(30)
class_size <- 25
birthdays <- sample(1:365, class_size, replace = TRUE)
sort(birthdays)

any(duplicated(birthdays))
birthdays[duplicated(birthdays)]

format(as.Date("2026-01-01") + birthdays[duplicated(birthdays)] - 1, "%B %d")

set.seed(31)
runs <- 2000
shared <- logical(runs)

for (run in 1:runs) {
  birthdays <- sample(1:365, class_size, replace = TRUE)
  shared[run] <- any(duplicated(birthdays))
}

mean(shared)

estimates <- numeric(5)

for (k in 1:5) {
  for (run in 1:runs) {
    birthdays <- sample(1:365, class_size, replace = TRUE)
    shared[run] <- any(duplicated(birthdays))
  }
  estimates[k] <- mean(shared)
}

estimates

exact <- 1 - prod((365 - 0:(class_size - 1)) / 365)
exact

set.seed(32)
sizes <- 2:60
estimate <- numeric(length(sizes))
exact <- numeric(length(sizes))

for (k in seq_along(sizes)) {
  for (run in 1:runs) {
    birthdays <- sample(1:365, sizes[k], replace = TRUE)
    shared[run] <- any(duplicated(birthdays))
  }
  estimate[k] <- mean(shared)
  exact[k] <- 1 - prod((365 - 0:(sizes[k] - 1)) / 365)
}

plot(sizes, estimate, pch = 16, col = "purple",
     ylim = c(0, 1), xlab = "Number of people in the room",
     ylab = "Chance of a shared birthday",
     main = "Shared birthdays: simulation and exact answer")
lines(sizes, exact, col = "darkgreen", lwd = 2)
abline(h = 0.5, lty = 3, col = "gray50")
legend("bottomright", c("Simulation (2,000 rooms each)", "Exact"),
       col = c("purple", "darkgreen"), pch = c(16, NA), lty = c(NA, 1),
       lwd = 2, bty = "n")

sizes[which(exact >= 0.5)[1]]
sizes[which(estimate >= 0.5)[1]]

set.seed(33)
day_weights <- rep(1, 365)
day_weights[182:273] <- 1.5

for (run in 1:runs) {
  birthdays <- sample(1:365, class_size, replace = TRUE, prob = day_weights)
  shared[run] <- any(duplicated(birthdays))
}

mean(shared)