PS 7_3

Question 4

library(tidyverse)

Warning: package 'ggplot2' was built under R version 4.5.2

── 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   4.0.1     ✔ 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

(sim <- parameters(~i, ~n, ~lambda,
                   c(3,5,9), c(10,15,20), c(0.05, 0.1)
                   )
  %>% add_trials(10000)
  %>% mutate(y_sample = pmap(list(i, n, lambda),
                             .f = \(i, n, l) rexp(n, rate = l)))
  %>% mutate(y_sorted = map(y_sample, sort))
  %>% mutate( yi = pmap_dbl(list(i, y_sorted), \(idx, y) pluck(y, idx)),
              yi_prev = pmap_dbl(list(i, y_sorted), \(idx, y) pluck(y, idx - 1)),
              T_stat = yi - yi_prev,
              f_t = (lambda * (n - i + 1)) * exp(-(lambda * (n - i + 1)) * T_stat))
)

# A tibble: 180,000 × 10
       i     n lambda .trial y_sample   y_sorted      yi yi_prev T_stat    f_t
   <dbl> <dbl>  <dbl>  <dbl> <list>     <list>     <dbl>   <dbl>  <dbl>  <dbl>
 1     3    10   0.05      1 <dbl [10]> <dbl [10]>  3.69   2.61  1.08   0.260 
 2     3    10   0.05      2 <dbl [10]> <dbl [10]>  6.68   6.61  0.0700 0.389 
 3     3    10   0.05      3 <dbl [10]> <dbl [10]>  8.20   4.77  3.43   0.101 
 4     3    10   0.05      4 <dbl [10]> <dbl [10]>  5.74   4.22  1.52   0.217 
 5     3    10   0.05      5 <dbl [10]> <dbl [10]>  3.76   1.89  1.87   0.189 
 6     3    10   0.05      6 <dbl [10]> <dbl [10]>  2.44   0.668 1.77   0.197 
 7     3    10   0.05      7 <dbl [10]> <dbl [10]>  4.94   1.30  3.64   0.0931
 8     3    10   0.05      8 <dbl [10]> <dbl [10]>  3.49   3.31  0.179  0.372 
 9     3    10   0.05      9 <dbl [10]> <dbl [10]>  7.06   4.22  2.84   0.128 
10     3    10   0.05     10 <dbl [10]> <dbl [10]>  8.02   0.529 7.49   0.0200
# ℹ 179,990 more rows

# Overlay
sim %>%
  ggplot(aes(x = T_stat)) +
  geom_histogram(aes(y = after_stat(density)), fill = "goldenrod", bins = 40) +
  geom_line(aes(y = f_t), color = "cornflowerblue", linewidth = 1) +
  facet_grid(n ~ i + lambda, scales = "free", labeller = label_both) +
  coord_cartesian(xlim = c(0, 20)) +
  theme_classic()