3.1 #8

library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.2.1     ✔ readr     2.2.0
✔ forcats   1.0.1     ✔ stringr   1.6.0
✔ ggplot2   4.0.3     ✔ tibble    3.3.1
✔ lubridate 1.9.5     ✔ tidyr     1.3.2
✔ purrr     1.2.2     
── 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

y <- 1:6
py <- (2*y - 1) / 36  

dd <- data.frame(y = y, probs = py)

(ggplot(data = dd)
  + geom_col(aes(x = y, y = probs), width = 0.1)
  + labs(y = 'P(Y=y)', x = 'y', title = 'pmf of Y (max of two die rolls)')
)

CDF_df <- data.frame(
  y  = c(-1, 0, 1,    2,    3,    4,     5,     6,   7),
  Fy = c(0,  0, 1/36, 4/36, 9/36, 16/36, 25/36, 1,   1)
)

(ggplot(data = CDF_df) 
  + geom_step(aes(x = y, y = Fy))
  + geom_point(aes(x = y, y = Fy), data = filter(CDF_df, y >= 1 & y <= 6))
  + labs(y = 'F(y) = P(Y <= y)', x = 'y', title = 'CDF of Y (max of two die rolls)')
)

N <- 10000
die <- 1:6

set.seed(82804)

many_dice_trials <- replicate(N, sample(die, 2, replace = TRUE), .as = two_rolls) %>%
  mutate(Y = map_dbl(.x = two_rolls, .f = \(dummy) max(dummy)))

many_dice_trials %>% head()

# A tibble: 6 × 3
  .trial two_rolls     Y
   <dbl> <list>    <dbl>
1      1 <int [2]>     3
2      2 <int [2]>     5
3      3 <int [2]>     6
4      4 <int [2]>     6
5      5 <int [2]>     4
6      6 <int [2]>     4

# Simulated mean, variance
many_dice_trials %>%
  summarize(mu_hat = mean(Y), E_Y2 = mean(Y^2), sigma2hat = var(Y))

# A tibble: 1 × 3
  mu_hat  E_Y2 sigma2hat
   <dbl> <dbl>     <dbl>
1   4.47  22.0      2.01

# Simulated median
many_dice_trials %>%
  summarize(median = quantile(Y, 0.5))

# A tibble: 1 × 1
  median
   <dbl>
1      5

# Simulated pmf
(ggplot(data = many_dice_trials)
  + geom_bar(aes(x = Y, y = after_stat(prop)), width = 0.5)
  + labs(y = 'Observed proportion', x = 'y',
         title = 'Simulated pmf of Y (max of two die rolls)')
)

# Simulated CDF
many_trials_with_ecdf <- many_dice_trials %>%
  mutate(empirical_CDF = cume_dist(Y))

many_trials_with_ecdf %>% head()

# A tibble: 6 × 4
  .trial two_rolls     Y empirical_CDF
   <dbl> <list>    <dbl>         <dbl>
1      1 <int [2]>     3         0.253
2      2 <int [2]>     5         0.688
3      3 <int [2]>     6         1    
4      4 <int [2]>     6         1    
5      5 <int [2]>     4         0.448
6      6 <int [2]>     4         0.448

(ggplot(data = many_trials_with_ecdf)
  + geom_step(aes(x = Y, y = empirical_CDF))
  + labs(y = 'Observed proportion', x = 'y',
         title = 'Simulated CDF of Y')
)