Problem 1

#Problem 1
(thinned_binom<- parameters(~alpha, ~p,
                             c(10, 20, 30), c(0.2, 0.4, 0.6)
                            )
  %>% add_trials(N)
  %>% mutate(X = map_int(alpha, \(l) rpois(1, l)))
  %>% mutate(Y = pmap_int(list(X, p), \(var1, var2) rbinom(1, size= var1, prob= var2)))
  %>% mutate(Fhat_Y = cume_dist(Y), .by = c(alpha, p))
  %>% mutate(F_Y = ppois(Y, alpha *p))
  )%>% head
## # A tibble: 6 × 7
##   alpha     p .trial     X     Y Fhat_Y   F_Y
##   <dbl> <dbl>  <dbl> <int> <int>  <dbl> <dbl>
## 1    10   0.2      1    11     3  0.858 0.857
## 2    10   0.2      2     7     1  0.413 0.406
## 3    10   0.2      3     8     0  0.138 0.135
## 4    10   0.2      4     7     3  0.858 0.857
## 5    10   0.2      5     9     1  0.413 0.406
## 6    10   0.2      6     6     1  0.413 0.406
(ggplot(aes(x=Y),data=thinned_binom)
+geom_step(aes(y=Fhat_Y, color = 'Simulated CDF'))
+geom_step(aes(y=F_Y, color = 'Analytic CDF'))
+labs(y=expression(P(Y<= y)), x= 'y', color=' ')
+theme_classic(base_size = 12)
+facet_grid(alpha~p, labeller = label_both)
)

(ggplot(aes(x=Fhat_Y, y = F_Y), data=thinned_binom)
  +geom_point()
  +labs(y=expression(P(Y <= y)), x= expression(hat(P)(Y <= y)))
  +geom_abline(intercept=0, slope=1)
  +theme_classic(base_size = 12)
  +facet_grid(alpha~p, labeller = label_both)
)

Problem 2a

(beta_sim <- parameters(~ alpha, ~ mu,
                       c(2, 4, 8, 16),
                       c(0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7))
  %>%add_trials(N) 
  %>%mutate(beta = alpha * (1 - mu) / mu) 
  %>%mutate(Y = pmap_dbl(list(alpha, beta), \(a, b) rbeta(1, a, b)))
)
## # A tibble: 280,000 × 5
##    alpha    mu .trial  beta      Y
##    <dbl> <dbl>  <dbl> <dbl>  <dbl>
##  1     2   0.1      1    18 0.126 
##  2     2   0.1      2    18 0.114 
##  3     2   0.1      3    18 0.176 
##  4     2   0.1      4    18 0.0743
##  5     2   0.1      5    18 0.106 
##  6     2   0.1      6    18 0.0469
##  7     2   0.1      7    18 0.191 
##  8     2   0.1      8    18 0.0111
##  9     2   0.1      9    18 0.101 
## 10     2   0.1     10    18 0.130 
## # ℹ 279,990 more rows
(beta_curves <- beta_sim
  %>%distinct(alpha, mu, beta)
  %>%mutate(x = list(seq(0, 1, length.out = 501)))
  %>%unnest(x) 
  %>%mutate(f_Y = dbeta(x, alpha, beta))
)
## # A tibble: 14,028 × 5
##    alpha    mu  beta     x   f_Y
##    <dbl> <dbl> <dbl> <dbl> <dbl>
##  1     2   0.1    18 0     0    
##  2     2   0.1    18 0.002 0.661
##  3     2   0.1    18 0.004 1.28 
##  4     2   0.1    18 0.006 1.85 
##  5     2   0.1    18 0.008 2.39 
##  6     2   0.1    18 0.01  2.88 
##  7     2   0.1    18 0.012 3.34 
##  8     2   0.1    18 0.014 3.77 
##  9     2   0.1    18 0.016 4.16 
## 10     2   0.1    18 0.018 4.52 
## # ℹ 14,018 more rows
(ggplot() +
  geom_histogram(data = beta_sim,
                 aes(x = Y, y = after_stat(density)),
                 binwidth = 0.02) +
  geom_line(data = beta_curves, aes(x = x, y = f_Y)) +
  facet_grid(alpha ~ mu, labeller = label_both)
)

Problem 2b

(var_tbl <- beta_sim
    %>% group_by(alpha, mu)
    %>% summarize(var_hat = var(Y), .groups = "drop"))
## # A tibble: 28 × 3
##    alpha    mu var_hat
##    <dbl> <dbl>   <dbl>
##  1     2   0.1 0.00429
##  2     2   0.2 0.0143 
##  3     2   0.3 0.0273 
##  4     2   0.4 0.0398 
##  5     2   0.5 0.0501 
##  6     2   0.6 0.0552 
##  7     2   0.7 0.0531 
##  8     4   0.1 0.00218
##  9     4   0.2 0.00762
## 10     4   0.3 0.0145 
## # ℹ 18 more rows

Problem 2c

(ggplot(var_tbl, aes(x = alpha, y = var_hat, color = factor(mu), group = mu)) 
    +geom_line() 
    +geom_point() 
    +labs(x = expression(alpha), y = "Variance", color = expression(mu)) 
    +theme_classic(base_size = 12))