Practice set 4.7

Introduction

Along with practicing parametric grid simulations, I also want you to practice assessing and analyzing your simulation results (a la Question 2). This type of analysis above and beyond simply obtaining simulation results will be required on future assignments (including, possibly, the coding portion of your midterm exam).

Submission instructions

Submit a link to your rendered html for this assignment.

Question 1

Recall the following problem. Let \(X\) represents the number of cars driving past a parking ramp in any given hour, and suppose that \(X\sim POI(\lambda)\). Given \(X=x\) cars driving past a ramp, let \(Y\) represent the number that decide to park in the ramp, with \(Y|X=x\sim BIN(x,p)\). We’ve shown analytically that unconditionally, \(Y\sim POI(\lambda p)\). Verify this result via a simulation study over a grid of \(\lambda \in \{10, 20, 30\}\) and \(p\in \{0.2, 0.4, 0.6\}\) using the framework presented in your notes, with 10,000 simulated outcomes per \((\lambda, p)\) combination. Create a faceted plot of the overlaid analytic and empirial CDFs as well as the p-p plot.

I’ll get you started:

library(tidyverse)
library(purrrfect)

N <- 10000
set.seed(227)
(thinned_pois <- parameters(~ lambda, ~ p,
                            c(10, 20, 30), c(0.2, 0.4, 0.6)
)
  %>% add_trials(N)
  %>% mutate(X = map_int(lambda, \(l) rpois(1, lambda = l)))
  %>% mutate(Y = pmap_int(list(X, p), \(var1, var2) rbinom(1, size = var1, prob = var2)))
  %>% mutate(Fhat_Y = cume_dist(Y), .by = c(lambda, p))
  %>% mutate(F_Y = ppois(Y, lambda = lambda*p))
) %>% head

(ggplot(aes(x = Y), data = thinned_pois)
  + 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(lambda ~ p, labeller = label_both, scales = 'free_x')
)

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

# A tibble: 6 × 7
  lambda     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

Question 2

In this problem we will study the relationship between \(\alpha\) and \(Var(Y)\) for fixed \(\mu\) for the beta distribution.

Recall that if \(Y\sim BETA(\alpha,\beta)\) and with \(\mu \equiv E(Y)\), \(\beta = \alpha \cdot \frac{1-\mu}{\mu}\).

Use this fact to simulate 10,000 realizations of \(Y\sim BETA(\alpha,\beta)\) for each combination of \(\alpha \in \{2,4,8,16\}\) and \(\mu\in\{0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7\}\).

N <- 10000
set.seed(227)
(beta_sim <- parameters(~ alpha, ~ mu,
                        c(2, 4, 8, 16), seq(0.1, 0.7, by = 0.1)
)
  %>% add_trials(N)
  %>% mutate(beta = alpha*(1-mu)/mu)
  %>% mutate(Y = pmap_dbl(list(alpha, beta), \(a, b) rbeta(1, shape1 = a, shape2 = b)))
  %>% mutate(f_Y = dbeta(Y, shape1 = alpha, shape2 = beta))
)

# A tibble: 280,000 × 6
   alpha    mu .trial  beta      Y   f_Y
   <dbl> <dbl>  <dbl> <dbl>  <dbl> <dbl>
 1     2   0.1      1    18 0.152   3.17
 2     2   0.1      2    18 0.0508  7.16
 3     2   0.1      3    18 0.0773  6.73
 4     2   0.1      4    18 0.168   2.52
 5     2   0.1      5    18 0.154   3.08
 6     2   0.1      6    18 0.0764  6.77
 7     2   0.1      7    18 0.0418  6.92
 8     2   0.1      8    18 0.0485  7.12
 9     2   0.1      9    18 0.189   1.84
10     2   0.1     10    18 0.170   2.46
# ℹ 279,990 more rows

A)

Create a faceted plot of the empirical densities overlaid with the analytic densities for each \((\alpha, \mu)\) combination.

(ggplot(aes(x = Y), data = beta_sim)
  + geom_histogram(aes(y = after_stat(density)),
                   fill = 'goldenrod',
                   binwidth = 0.02, center = 0.01)
  + geom_line(aes(y = f_Y), col = 'cornflowerblue', linewidth = 1)
  + labs(y = 'Beta density', x = 'y')
  + theme_classic(base_size = 10)
  + facet_grid(alpha ~ mu, labeller = label_both, scales = 'free_y')
)

B)

Use summarize() to find the simulated variance of \(Y\) for each \((\alpha,\mu)\) combination.

(beta_var <- beta_sim
  %>% summarize(var_Y = var(Y), .by = c(alpha, mu))
)

# A tibble: 28 × 3
   alpha    mu   var_Y
   <dbl> <dbl>   <dbl>
 1     2   0.1 0.00424
 2     2   0.2 0.0142 
 3     2   0.3 0.0278 
 4     2   0.4 0.0405 
 5     2   0.5 0.0505 
 6     2   0.6 0.0562 
 7     2   0.7 0.0536 
 8     4   0.1 0.00223
 9     4   0.2 0.00778
10     4   0.3 0.0148 
# ℹ 18 more rows

C)

Create a plot of the simulated variances as a function of \(\alpha\). Use lines as your geometry, and color-code each line by \(\mu\).

(ggplot(aes(x = alpha, y = var_Y, color = factor(mu)), data = beta_var)
  + geom_line()
  + geom_point()
  + scale_x_continuous(breaks = c(2, 4, 8, 16))
  + labs(x = expression(alpha), y = 'Simulated Var(Y)', color = expression(mu))
  + theme_classic(base_size = 12)
)

D)

Comment on how the combination of \(\alpha\) and \(\mu\) impact variance.

For a fixed μ, the variance of Y goes down as α increases, with the biggest drop happening between α = 2 and α = 4. After that, the changes get smaller. Since α + β = α/μ when μ is fixed, a larger α makes the distribution more concentrated around its mean. For a fixed α, the variance generally increases as μ goes from 0.1 to 0.6, but the curves for μ = 0.6 and 0.7 are almost the same and cross each other. This happens because the numerator of the variance formula reaches its highest point around μ = 2/3. The differences in variance between the μ values are also bigger when α = 2 and get smaller as α increases, meaning μ has a bigger effect on variance when α is small.