# ============================================================================ # 04-04 Probability Simulator # Chapters 05-01, 05-02, 05-03 # # Run from the project folder: Rscript r/analysis.R # ============================================================================ set.seed(2026) rolls <- read.csv("data/dice_rolls.csv") n <- nrow(rolls) # ── 1. EXACT DISTRIBUTION OF THE SUM OF TWO DICE ──────────────────────────── cat("\n===== 1. THEORETICAL vs EMPIRICAL =====\n") S <- expand.grid(d1 = 1:6, d2 = 1:6) S$sum <- S$d1 + S$d2 theory <- as.numeric(table(S$sum)) / 36 emp <- as.numeric(table(factor(rolls$sum, levels = 2:12))) / n comp <- data.frame( sum = 2:12, ways = as.integer(table(S$sum)), theoretical = round(theory, 4), empirical = round(emp, 4), difference = round(emp - theory, 4) ) print(comp, row.names = FALSE) cat(sprintf("\nCHECK sum of theoretical probabilities = %.4f\n", sum(theory))) cat(sprintf("Largest gap between theory and %d simulated rolls: %.4f\n", n, max(abs(emp - theory)))) # ── 2. PROBABILITY RULES ON THE DICE ──────────────────────────────────────── cat("\n===== 2. PROBABILITY RULES =====\n") cat(sprintf("P(sum = 7) theory %.4f empirical %.4f\n", mean(S$sum == 7), mean(rolls$sum == 7))) cat(sprintf("P(sum >= 10) theory %.4f empirical %.4f\n", mean(S$sum >= 10), mean(rolls$sum >= 10))) cat(sprintf("P(doubles) theory %.4f empirical %.4f\n", mean(S$d1 == S$d2), mean(rolls$die1 == rolls$die2))) cat(sprintf("P(at least one 6) theory %.4f empirical %.4f\n", mean(S$d1 == 6 | S$d2 == 6), mean(rolls$die1 == 6 | rolls$die2 == 6))) cat(sprintf(" via the complement: 1 - (5/6)^2 = %.4f\n", 1 - (5/6)^2)) cat(sprintf("P(sum = 7 OR doubles) = P(7)+P(dbl)-P(both) = %.4f\n", mean(S$sum == 7) + mean(S$d1 == S$d2) - mean(S$sum == 7 & S$d1 == S$d2))) cat(sprintf("P(sum >= 10 | doubles) = %.4f\n", mean(S$sum >= 10 & S$d1 == S$d2) / mean(S$d1 == S$d2))) # ── 3. LAW OF LARGE NUMBERS ───────────────────────────────────────────────── cat("\n===== 3. LAW OF LARGE NUMBERS =====\n") running <- cumsum(rolls$sum == 7) / seq_len(n) for (k in c(10, 50, 100, 500, 1000)) { cat(sprintf(" after %5d rolls: %.4f (error %.4f)\n", k, running[k], abs(running[k] - 6/36))) } # ── 4. CONTINGENCY TABLE ──────────────────────────────────────────────────── cat("\n===== 4. CONTINGENCY TABLE =====\n") tab <- matrix(c(45, 55, 20, 15, 30, 35), nrow = 2, byrow = TRUE, dimnames = list(exercise = c("Exercises", "Does not"), health = c("Excellent", "Good", "Poor"))) print(addmargins(tab)) N <- sum(tab) cat(sprintf("\nmarginal P(Exercises) = %.4f\n", sum(tab[1, ]) / N)) cat(sprintf("marginal P(Excellent) = %.4f\n", sum(tab[, 1]) / N)) cat(sprintf("joint P(Exercises AND Excellent) = %.4f\n", tab[1, 1] / N)) cat(sprintf("conditional P(Excellent | Exercises) = %.4f <- divide by the ROW total\n", tab[1, 1] / sum(tab[1, ]))) cat(sprintf("conditional P(Exercises | Excellent) = %.4f <- divide by the COLUMN total\n", tab[1, 1] / sum(tab[, 1]))) cat(sprintf("union P(Exercises OR Excellent) = %.4f\n", (sum(tab[1, ]) + sum(tab[, 1]) - tab[1, 1]) / N)) expected <- outer(rowSums(tab), colSums(tab)) / N cat("\nExpected counts IF independent:\n"); print(round(expected, 2)) cat(sprintf("Observed 45 vs expected %.1f -> the variables are DEPENDENT\n", expected[1, 1])) # ── 5. BAYES' THEOREM ─────────────────────────────────────────────────────── cat("\n===== 5. BAYES' THEOREM =====\n") bayes <- function(prior, sens, fpr) { evidence <- sens * prior + fpr * (1 - prior) c(evidence = evidence, posterior = sens * prior / evidence) } sens <- 0.95; fpr <- 0.10 cat(sprintf("Test: %.0f%% sensitive, false-positive rate %.0f%%\n\n", 100 * sens, 100 * fpr)) cat(" prevalence P(positive) P(disease | positive)\n") for (p in c(0.001, 0.005, 0.01, 0.05, 0.10, 0.30, 0.50)) { b <- bayes(p, sens, fpr) cat(sprintf(" %6.3f %6.4f %6.4f\n", p, b[1], b[2])) } cat("\nNatural frequencies at 1% prevalence, per 10,000 people:\n") p <- 0.01 cat(sprintf(" with disease %6.0f of whom %5.0f test positive\n", 10000 * p, 10000 * p * sens)) cat(sprintf(" without disease %6.0f of whom %5.0f test positive\n", 10000 * (1 - p), 10000 * (1 - p) * fpr)) cat(sprintf(" => P(disease | positive) = %.0f / %.0f = %.4f\n", 10000 * p * sens, 10000 * p * sens + 10000 * (1 - p) * fpr, bayes(p, sens, fpr)[2])) # ── 6. COUNTING ───────────────────────────────────────────────────────────── cat("\n===== 6. COUNTING =====\n") cat(sprintf("5! = %d\n", factorial(5))) cat(sprintf("10P3 (order matters) = %d\n", factorial(10) / factorial(7))) cat(sprintf("10C3 (order does not) = %d\n", choose(10, 3))) cat(sprintf("49C6 lottery combos = %s\n", format(choose(49, 6), big.mark = ","))) cat(sprintf("P(jackpot) = %.3e\n", 1 / choose(49, 6))) cat(sprintf("P(match exactly 5) = %.3e (1 in %s)\n", choose(6,5)*choose(43,1)/choose(49,6), format(round(choose(49,6)/(choose(6,5)*choose(43,1))), big.mark = ","))) cat(sprintf("MISSISSIPPI arrangements = %s\n", format(factorial(11)/prod(factorial(c(1,4,4,2))), big.mark = ","))) bday <- function(k) 1 - prod((365 - 0:(k - 1)) / 365) cat("\nBirthday problem:\n") for (k in c(10, 23, 30, 50, 70)) cat(sprintf(" k = %2d -> %.4f\n", k, bday(k))) # ── 7. MONTY HALL ─────────────────────────────────────────────────────────── cat("\n===== 7. MONTY HALL (100,000 simulated games) =====\n") trials <- 1e5 car <- sample(1:3, trials, replace = TRUE) pick <- sample(1:3, trials, replace = TRUE) stay_wins <- mean(pick == car) switch_wins <- mean(pick != car) # switching wins exactly when the first pick was wrong cat(sprintf(" stay: %.4f (theory 1/3 = %.4f)\n", stay_wins, 1/3)) cat(sprintf(" switch: %.4f (theory 2/3 = %.4f)\n", switch_wins, 2/3)) cat(" Switching doubles your chance -- the host's choice is not independent\n") cat(" of where the car is, which is what makes the intuition fail.\n") # ── 8. PLOTS ──────────────────────────────────────────────────────────────── png("probability_plots.png", width = 1100, height = 850) par(mfrow = c(2, 2), mar = c(4.5, 4.5, 3, 1)) barplot(rbind(theory, emp), beside = TRUE, names.arg = 2:12, col = c("#5B2A86", "#0FA3A3"), border = NA, main = "Two dice: theory vs 1,000 rolls", xlab = "Sum", ylab = "Probability") legend("topright", c("theoretical", "empirical"), fill = c("#5B2A86", "#0FA3A3"), bty = "n") plot(running, type = "l", col = "#5B2A86", ylim = c(0, 0.35), main = "Law of Large Numbers", xlab = "Rolls", ylab = "Running P(sum = 7)") abline(h = 6/36, col = "#0FA3A3", lwd = 2, lty = 2) prev <- seq(0.001, 0.5, length.out = 300) post <- sapply(prev, function(p) bayes(p, sens, fpr)[2]) plot(prev, post, type = "l", col = "#5B2A86", lwd = 2, main = "Bayes: posterior vs base rate", xlab = "Prevalence", ylab = "P(disease | positive)") abline(v = 0.01, lty = 2, col = "grey50") ks <- 1:60 plot(ks, sapply(ks, bday), type = "l", col = "#0B7A7A", lwd = 2, main = "Birthday problem", xlab = "People in the room", ylab = "P(at least one shared birthday)") abline(h = 0.5, lty = 2, col = "grey50"); abline(v = 23, lty = 2, col = "grey50") par(mfrow = c(1, 1)) dev.off() cat("\nWrote probability_plots.png\n")