Chapter 10.3

Question 8

library(purrrfect)
library(tidyverse)

var_0 = 1
alpha = 0.05

(two_tests <- parameters(~var_true, ~n,
                         seq(1, 10, l = 100), c(5,10,15,20))
  %>% mutate(crit_value_T1 = qchisq(0.95, df = n),
             crit_value_T2 = qnorm(0.95, mean = 0, sd = 1),
             prob_T1_rejects = 1 - pchisq(crit_value_T1 / var_true, df = n),
             prob_T2_rejects = 1 - pnorm(crit_value_T2 / sqrt(var_true), mean = 0, sd = 1))
) %>% head

# A tibble: 6 × 6
  var_true     n crit_value_T1 crit_value_T2 prob_T1_rejects prob_T2_rejects
     <dbl> <dbl>         <dbl>         <dbl>           <dbl>           <dbl>
1     1        5          11.1          1.64          0.0500          0.0500
2     1       10          18.3          1.64          0.0500          0.0500
3     1       15          25.0          1.64          0.0500          0.0500
4     1       20          31.4          1.64          0.0500          0.0500
5     1.09     5          11.1          1.64          0.0711          0.0576
6     1.09    10          18.3          1.64          0.0793          0.0576

(ggplot(data = two_tests)
  + geom_line(aes(x = var_true, y = prob_T1_rejects, col = 'Test 1'))
  + geom_line(aes(x = var_true, y = prob_T2_rejects, col = 'Test 2'))
  + facet_wrap(~n, labeller = label_both)
  + scale_y_continuous(breaks = c(0.05, seq(0.25, 1, by = 0.25)))
  + geom_hline(aes(yintercept = 0.05), linetype = 2)
  + theme_classic()
  + labs(x = expression(sigma[true]^2),
         y = 'P(reject null)', color = '')
)

We see that when the true value of \(\sigma^2\) is 1, both tests have a rejection rate of 0.05, indicating they both have the same size. We see that for every other value of \(\sigma^2\), Test 1 has more power than Test 2, showing that T1 is uniformly more powerful than T2.

Question 9

alpha <- 0.05
theta_0 <- 0.2

# simulate critical values
(beta_crit_values <- parameters(~n, c(5, 10, 15, 20))
  %>% add_trials(10000)
  %>% mutate(YSample = map(n, .f = \(n) rbeta(n, theta_0, 1)))
  %>% mutate(T1 = map_dbl(YSample, \(y) sum(log(y)) / length(y)),
             T2 = map_dbl(YSample, mean)
             )
  %>% summarize(c1_simulated = quantile(T1, 0.95),
                c2_simulated = quantile(T2, 0.95),
                .by = n)
) %>% head()

# A tibble: 4 × 3
      n c1_simulated c2_simulated
  <dbl>        <dbl>        <dbl>
1     5        -2.01        0.365
2    10        -2.71        0.306
3    15        -3.09        0.276
4    20        -3.29        0.265

# join simulated critical values
options(width = 10000)
(beta_decisions <- parameters(~n, ~true_theta,
                              c(5, 10, 15, 20), seq(0.2, 1, l = 20)
                              )
  %>% add_trials(10000)
  %>% inner_join(., beta_crit_values, by = 'n')
  %>% mutate(Ysample = pmap(list(n, true_theta), \(n, t) rbeta(n, t, 1)),
             T1 = map_dbl(Ysample, \(y) sum(log(y)) / length(y)),
             T2 = map_dbl(Ysample, mean))
  %>% mutate(T1_decision = ifelse(T1 > c1_simulated, 1, 0),
             T2_decision = ifelse(T2 > c2_simulated, 1, 0))
  %>% summarize(prob_T1_rejects = mean(T1_decision),
                prob_T2_rejects = mean(T2_decision),
                .by = c(true_theta, n))
) %>% head()

# A tibble: 6 × 4
  true_theta     n prob_T1_rejects prob_T2_rejects
       <dbl> <dbl>           <dbl>           <dbl>
1      0.2       5          0.0555          0.0588
2      0.242     5          0.103           0.0905
3      0.284     5          0.159           0.136 
4      0.326     5          0.238           0.170 
5      0.368     5          0.308           0.220 
6      0.411     5          0.400           0.272 

# visualize
(ggplot(data = beta_decisions)
  + geom_line(aes(x = true_theta, y = prob_T1_rejects, col = 'Test 1'))
  + geom_line(aes(x = true_theta, y = prob_T2_rejects, col = 'Test 2'))
  + facet_wrap(~n, labeller = label_both)
  + scale_y_continuous(breaks = c(0.05, seq(0.25, 1, by = 0.25)))
  + geom_hline(aes(yintercept = 0.05), linetype = 2)
  + theme_classic()
  + labs(x = expression(theta[true]),
         y = 'P(reject null)', color = '')
)

We see that when the true value of \(\theta\) is 0.2, both tests have a rejection rate of 0.05, indicating they both have the same size. We see that for every other value of \(theta\), Test 1 has more power than Test 2, showing that T1 is uniformly more powerful than T2.