8.2 Problems 5 and 6

library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.6
✔ forcats   1.0.1     ✔ stringr   1.6.0
✔ ggplot2   4.0.1     ✔ tibble    3.3.1
✔ lubridate 1.9.4     ✔ tidyr     1.3.2
✔ purrr     1.2.1     
── 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

library(purrr)

B <- 10000
theta0 <- 2
n_grid <- c(5, 10, 20, 40, 80, 160, 320)
z <- qnorm(0.975)

one_run <- function(n) {
  y <- rbeta(n, shape1 = theta0, shape2 = 1)
  M <- max(y)
  ybar <- mean(y)
  s <- sd(y)

  a <- 0.025^(1/n)
  b <- 0.975^(1/n)
  exact_L <- log(b) / log(M)
  exact_U <- log(a) / log(M)

  cover_exact <- (exact_L <= theta0 && theta0 <= exact_U)

  mu_L <- ybar - z * s / sqrt(n)
  mu_U <- ybar + z * s / sqrt(n)

  eps <- 1e-8
  mu_L <- pmin(pmax(mu_L, eps), 1 - eps)
  mu_U <- pmin(pmax(mu_U, eps), 1 - eps)

  asy_L <- mu_L / (1 - mu_L)
  asy_U <- mu_U / (1 - mu_U)

  cover_asy <- (asy_L <= theta0 && theta0 <= asy_U)

  tibble(cover_exact = cover_exact, cover_asy = cover_asy)
}

coverage_df <- map_dfr(n_grid, function(n) {
  sims <- replicate(B, one_run(n), simplify = FALSE) %>% bind_rows()
  tibble(
    n = n,
    exact = mean(sims$cover_exact),
    asymptotic = mean(sims$cover_asy)
  )
})

coverage_long <- coverage_df %>%
  pivot_longer(cols = c(exact, asymptotic), names_to = "method", values_to = "coverage")

coverage_long %>%
  ggplot(aes(x = n, y = coverage, group = method)) +
  geom_line() +
  geom_point() +
  geom_hline(yintercept = 0.95, linetype = 2) +
  scale_x_log10(breaks = n_grid) +
  labs(title = "Coverage vs n (true theta = 2)", y = "coverage", x = "n") +
  theme_classic()

library(tidyverse)
library(purrr)

B <- 10000
n_grid <- c(5, 20, 50, 100)

pops <- tribble(
  ~pop,           ~rfun,                                     ~mu_true, ~sig2_true,
  "UNIF(0,8)",     function(n) runif(n, 0, 8),                4,        16/3,
  "EXP(scale=4)",  function(n) rexp(n, rate = 1/4),           4,        16,
  "GEO(p=0.2)",    function(n) rgeom(n, prob = 0.2),          4,        20
)

cover_once <- function(y, mu_true, sig2_true) {
  n <- length(y)
  ybar <- mean(y)
  s <- sd(y)
  s2 <- var(y)

  tcrit <- qt(0.975, df = n - 1)
  mu_L <- ybar - tcrit * s / sqrt(n)
  mu_U <- ybar + tcrit * s / sqrt(n)
  cover_mu <- (mu_L <= mu_true && mu_true <= mu_U)

  chi_L <- qchisq(0.025, df = n - 1)
  chi_U <- qchisq(0.975, df = n - 1)

  sig2_L <- (n - 1) * s2 / chi_U
  sig2_U <- (n - 1) * s2 / chi_L
  cover_sig2 <- (sig2_L <= sig2_true && sig2_true <= sig2_U)

  c(cover_mu = cover_mu, cover_sig2 = cover_sig2)
}

results <- pops %>%
  crossing(n = n_grid) %>%
  mutate(sim = pmap(list(rfun, mu_true, sig2_true, n), function(rfun, mu0, v0, n) {
    mat <- replicate(B, {
      y <- rfun(n)
      cover_once(y, mu0, v0)
    })
    tibble(
      cover_mu = mean(mat["cover_mu", ]),
      cover_sig2 = mean(mat["cover_sig2", ])
    )
  })) %>%
  unnest(sim)

mu_table <- results %>%
  select(pop, n, mu_true, cover_mu) %>%
  pivot_wider(names_from = pop, values_from = cover_mu) %>%
  arrange(n)

sig2_table <- results %>%
  select(pop, n, sig2_true, cover_sig2) %>%
  pivot_wider(names_from = pop, values_from = cover_sig2) %>%
  arrange(n)

mu_table

# A tibble: 4 × 5
      n mu_true `EXP(scale=4)` `GEO(p=0.2)` `UNIF(0,8)`
  <dbl>   <dbl>          <dbl>        <dbl>       <dbl>
1     5       4          0.882        0.883       0.940
2    20       4          0.918        0.918       0.947
3    50       4          0.931        0.929       0.950
4   100       4          0.940        0.943       0.952

sig2_table

# A tibble: 12 × 5
       n sig2_true `EXP(scale=4)` `GEO(p=0.2)` `UNIF(0,8)`
   <dbl>     <dbl>          <dbl>        <dbl>       <dbl>
 1     5     16             0.823       NA          NA    
 2     5     20            NA            0.823      NA    
 3     5      5.33         NA           NA           0.986
 4    20     16             0.735       NA          NA    
 5    20     20            NA            0.726      NA    
 6    20      5.33         NA           NA           0.996
 7    50     16             0.702       NA          NA    
 8    50     20            NA            0.703      NA    
 9    50      5.33         NA           NA           0.998
10   100     16             0.689       NA          NA    
11   100     20            NA            0.686      NA    
12   100      5.33         NA           NA           0.998