PS 8.1

9. Consider the previous question. Use simulation studies to compare the three estimators θ ̂_1, θ ̂_2, and θ ̂_3 of the median of an exponential when β=5. Use 10,000 replications. Create a summary table structured like the one below, with all results rounded to at most 3 decimal places. Full credit only given if your table is clearly readable in your rendered document, making life easier for both you and anyone you share the table with to interpret it. Discuss the results of your simulation study, comparing the properties of these estimators.

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

library(purrrfect)

Attaching package: 'purrrfect'

The following objects are masked from 'package:base':

    replicate, tabulate

theta <- 5 * log(2)

sim_data <-
  parameters(~n,
             c(3, 11, 51)) %>%
  add_trials(10000) %>%
  mutate(
    Ysample = map(n, ~ rexp(.x, rate = 1/5)),
    ybar    = map_dbl(Ysample, mean),
    ymed    = map_dbl(Ysample, median)
  ) %>%
  mutate(
    theta_hat1 = ybar,
    theta_hat2 = log(2) * ybar,
    theta_hat3 = ymed
  )

summary_table <-
  sim_data %>%
  pivot_longer(theta_hat1:theta_hat3,
               names_to = "estimator",
               values_to = "estimate") %>%
  summarize(
    bias     = mean(estimate - theta),
    variance = var(estimate),
    mse      = bias^2 + variance,
    .by = c(n, estimator)
  ) %>%
  mutate(across(c(bias, variance, mse), round, 3))

Warning: There was 1 warning in `mutate()`.
ℹ In argument: `across(c(bias, variance, mse), round, 3)`.
Caused by warning:
! The `...` argument of `across()` is deprecated as of dplyr 1.1.0.
Supply arguments directly to `.fns` through an anonymous function instead.

  # Previously
  across(a:b, mean, na.rm = TRUE)

  # Now
  across(a:b, \(x) mean(x, na.rm = TRUE))

library(knitr)
kable(summary_table)

n	estimator	bias	variance	mse
3	theta_hat1	1.521	8.260	10.573
3	theta_hat2	-0.009	3.968	3.969
3	theta_hat3	0.715	9.346	9.857
11	theta_hat1	1.506	2.278	4.547
11	theta_hat2	-0.019	1.095	1.095
11	theta_hat3	0.200	2.293	2.333
51	theta_hat1	1.534	0.501	2.856
51	theta_hat2	0.000	0.241	0.241
51	theta_hat3	0.055	0.503	0.506