rm(list = ls())

# Exact probability
n <- 3
p <- 1/6
prob_exact <- 1 - dbinom(0, n, p)

# Simulation
set.seed(123)
n_sim <- 10000
results <- replicate(n_sim, {
  dice <- sample(1:6, size = 3, replace = TRUE)
  any(dice == 3)
})
prob_sim <- mean(results)

# Comparison
cbind(prob_exact, prob_sim)