6.1 code problem 9

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.2     ✔ tibble    3.3.0
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.1.0     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

lambda_grid <- c(0.5, 1, 1.5, 2)
n_rep <- 10000

sim_df <- map_dfr(
  lambda_grid,
  ~ tibble(
    lambda = .x,
    X = rexp(n_rep, rate = .x),
    Y = rexp(n_rep, rate = .x)
  ) %>%
    mutate(U = X + Y)
)
theory_df <- crossing(
  lambda = lambda_grid,
  u = seq(0, 15, length.out = 400)
) %>%
  mutate(density = dgamma(u, shape = 2, rate = lambda))

ggplot() +
  geom_histogram(
    data = sim_df,
    aes(x = U, y = after_stat(density)),
    bins = 40
  ) +
  geom_line(
    data = theory_df,
    aes(x = u, y = density),
    linewidth = 1
  ) +
  facet_wrap(~ lambda, scales = "free") +
  labs(
    x = "U = X + Y",
    y = "Density",
    title = "Simulated sums of two Exp(λ) vs Gamma(2, λ) density"
  ) +
  theme_classic()