Simulation exercise: solutions

This document includes only one possible solution. Yours may vary. For the simulations, your specific quantities will also vary.

Setup (given)

library(tidyverse)

set.seed(365)

counties <- left_join(
    read_csv("https://princeton-ddss.github.io/r-bootcamp/files/county_data.csv"),
    read_csv("https://princeton-ddss.github.io/r-bootcamp/files/county_population.csv"),
    join_by(GEOID == fips)
)

Part 1: the birthday problem

One room

# Without `replace = TRUE`, no two people *can* share a birthday.
room <- sample(1:365, size = 23, replace = TRUE)

# `anyDuplicated()` returns the index of the first duplicate, or `0` if there isn't one, so comparing to zero returns T/F
anyDuplicated(room) > 0

[1] FALSE

## this also works; returns T/F:
#  any(duplicated(room))

Wrap it in a function

shared_birthday <- \(k) {
    birthdays <- sample(1:365, size = k, replace = TRUE)
    anyDuplicated(birthdays) > 0
}

replicate(5, shared_birthday(23))

[1] FALSE FALSE  TRUE  TRUE FALSE

Many rooms

rooms <- replicate(10000, shared_birthday(23))

mean(rooms)

[1] 0.5071

How many replications did you need?

reps <- c(10, 100, 1000, 10000, 100000)

ests <- sapply(
    reps,
    \(r) replicate(r, shared_birthday(23)) |> mean()
)

tibble(
    replications = reps,
    estimate = ests
)

# A tibble: 5 × 2
  replications estimate
         <dbl>    <dbl>
1           10    0.6  
2          100    0.46 
3         1000    0.518
4        10000    0.505
5       100000    0.508

Sweep across room sizes

ks <- 2:60
n <- 1000

prob <- sapply(
    ks,
    \(k) replicate(n, shared_birthday(k)) |> mean()
)

birthday_curve <- tibble(
    k = ks,
    prob = prob
)

head(birthday_curve, 3)

# A tibble: 3 × 2
      k  prob
  <int> <dbl>
1     2 0    
2     3 0.008
3     4 0.013

Plot it

ggplot(birthday_curve, aes(x = k, y = prob)) +
    geom_line(linewidth = 1) +
    geom_hline(yintercept = 0.5, linetype = "dotted") +
    geom_vline(xintercept = 23, linetype = "dotted") +
    labs(
        x = "People in the room",
        y = "Pr(at least one shared birthday)",
        title = "The birthday problem, simulated"
    ) +
    theme_minimal(base_size = 15)

min(birthday_curve$k[birthday_curve$prob > 0.5])

[1] 23

Check against the exact answer

exact_prob <- function(k) 1 - prod((365 - 0:(k - 1)) / 365)

birthday_curve <- birthday_curve |>
    mutate(exact = sapply(k, exact_prob))

exact_prob(23)

[1] 0.5072972

birthday_curve |>
    pivot_longer(
        cols = c(prob, exact),
        names_to = "source",
        values_to = "p"
    ) |>
    mutate(source = if_else(
        source == "prob",
        "Simulated",
        "Exact"
    )) |>
    ggplot(aes(x = k, y = p, color = source)) + ## switching to + now
    geom_line(linewidth = 1) +
    labs(
        x = "People in the room",
        y = "Pr(at least one shared birthday)",
        color = NULL
    ) +
    theme_minimal(base_size = 15) +
    theme(legend.position = "bottom")

The two curves are hard to tell apart. The largest gap anywhere on the sweep is given by:

round(max(abs(birthday_curve$prob - birthday_curve$exact)), 3)

Part 2: sampling and the bootstrap

The population

pop_all <- counties$pop[!is.na(counties$pop)]

ggplot(tibble(pop = pop_all), aes(x = pop)) +
    geom_histogram(bins = 60) +
    scale_x_continuous(
        labels = scales::label_number(scale_cut = scales::cut_short_scale())
    ) +
    labs(x = "County population", y = "Count") +
    theme_minimal(base_size = 15)

## note that the data are extremely skewed

mean(pop_all)

[1] 104575.2

sd(pop_all)

[1] 335250.6

Sampling distributions

sizes <- c(2, 5, 30, 100, 1000)

samp <- map_dfr(sizes, \(n) {
    tibble(
        n = n,
        xbar = replicate(10000, sample(pop_all, size = n, replace = T) |> mean())
    )
})

samp |>
    mutate(n = factor(n, levels = sizes, labels = paste("n =", sizes))) |>
    ggplot(aes(x = xbar)) +
    geom_histogram(bins = 50) +
    # The x-axes are **free**, so read the axis labels rather than the apparent widths.
    facet_wrap(~n, scales = "free", nrow = 1) +
    scale_x_continuous(
        labels = scales::label_number(scale_cut = scales::cut_short_scale())
    ) +
    labs(x = "Sample mean of county population", y = "Count") +
    theme_minimal(base_size = 13) +
    theme(axis.text.x = element_text(angle = 45, hjust = 1))

The width of those distributions

samp |>
    summarize(sim_sd = sd(xbar), .by = n) |>
    mutate(
        theory = sd(pop_all) / sqrt(n),
        ratio = sim_sd / theory
    )

# A tibble: 5 × 4
      n  sim_sd  theory ratio
  <dbl>   <dbl>   <dbl> <dbl>
1     2 246003. 237058. 1.04 
2     5 147210. 149929. 0.982
3    30  60133.  61208. 0.982
4   100  33912.  33525. 1.01 
5  1000  10638.  10602. 1.00 

Every simulated spread lands within about 5% of \(\sigma/\sqrt{n}\), but n >= 1000 gets us close.

Note the \(\sqrt{\ }\): to halve your uncertainty you need four times the data.

The middle 95%

quantile(samp$xbar[samp$n == 100], c(0.025, 0.975))

     2.5%     97.5% 
 57467.46 192179.57 

95% of samples of 100 counties give a mean between those two numbers. A sample landing outside happens one time in twenty.

The bootstrap

my_counties <- counties |>
    filter(!is.na(income)) |>
    slice_sample(n = 50)

mean(my_counties$income)

[1] 65099.86

boot50 <- replicate(
    10000,
    sample(my_counties$income, size = 50, replace = T) |>
        mean()
)

ggplot(tibble(m = boot50), aes(x = m)) +
    geom_histogram(bins = 40) +
    labs(x = "Bootstrap mean income", y = "Count") +
    theme_minimal(base_size = 15)

sd(boot50)

[1] 2116.174

quantile(boot50, c(0.025, 0.975))

    2.5%    97.5% 
61125.75 69343.31