ps10_3

Question 8

library(tidyverse)

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

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

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

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.2.0
✔ 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

power_data <- parameters(~sigma2_true, ~n, seq(1, 10, l = 100), c(5, 10, 15, 20)) %>%
  mutate(
    prob_phi1 = 1 - pchisq(qchisq(0.95, df = n) / sigma2_true, df = n),
    prob_phi2 = 1 - pnorm(qnorm(0.95) / sqrt(sigma2_true))
  )

ggplot(power_data, aes(x = sigma2_true)) +
  geom_line(aes(y = prob_phi1, color = "phi_1 (UMP)")) +
  geom_line(aes(y = prob_phi2, color = "phi_2 (Mean)")) +
  facet_wrap(~n, labeller = label_both) +
  geom_hline(yintercept = 0.05, linetype = 2) +
  scale_y_continuous(breaks = c(0.05, seq(0.25, 1, by = 0.25))) +
  theme_classic() +
  labs(x = expression(sigma[true]^2), y = "P(reject null)", color = "")

• Both keep the 5% error rate, but the UMP test’s power dominates and scales with sample size while the mean test stays pretty flat.

Question 9

library(tidyverse)
library(purrrfect)

theta0 <- 0.2

(many_nulls <- parameters(~n, c(5, 10, 15, 20)) %>% 
  add_trials(10000) %>% 
  mutate(
    y_sample = map(n, .f = \(n) rbeta(n, shape1 = theta0, shape2 = 1)),
    T1 = map_dbl(y_sample, .f = \(y) mean(log(y))),
    T2 = map_dbl(y_sample, .f = \(y) mean(y))
  )
) %>% head()

# A tibble: 6 × 5
      n .trial y_sample     T1      T2
  <dbl>  <dbl> <list>    <dbl>   <dbl>
1     5      1 <dbl [5]> -4.91 0.0401 
2     5      2 <dbl [5]> -9.00 0.00398
3     5      3 <dbl [5]> -4.59 0.107  
4     5      4 <dbl [5]> -2.15 0.238  
5     5      5 <dbl [5]> -1.94 0.277  
6     5      6 <dbl [5]> -6.49 0.0824 

(crit_vals <- many_nulls %>% 
  group_by(n) %>% 
  summarize(
    c_T = quantile(T1, 0.95),
    c_Y = quantile(T2, 0.95),
    .groups = "drop"
  )
)

# A tibble: 4 × 3
      n   c_T   c_Y
  <dbl> <dbl> <dbl>
1     5 -1.95 0.379
2    10 -2.68 0.312
3    15 -3.08 0.282
4    20 -3.30 0.263

# Power on grid of theta(a)
(power_sim <- parameters(~n, ~theta, c(5, 10, 15, 20), seq(0.2, 1, l = 20)) %>% 
  left_join(crit_vals, by = "n") %>% 
  add_trials(10000) %>% 
  mutate(
    y_sample = map2(n, theta, .f = \(n, theta) rbeta(n, shape1 = theta, shape2 = 1)),
    T1 = map_dbl(y_sample, .f = \(y) mean(log(y))),
    T2 = map_dbl(y_sample, .f = \(y) mean(y)),
    phi1_rejects = T1 > c_T,
    phi2_rejects = T2 > c_Y
  ) %>% 
  group_by(n, theta) %>% 
  summarize(
    prob_phi1 = mean(phi1_rejects),
    prob_phi2 = mean(phi2_rejects),
    .groups = "drop"
  )
) %>% head()

# A tibble: 6 × 4
      n theta prob_phi1 prob_phi2
  <dbl> <dbl>     <dbl>     <dbl>
1     5 0.2      0.0485    0.0483
2     5 0.242    0.0925    0.0766
3     5 0.284    0.150     0.110 
4     5 0.326    0.222     0.155 
5     5 0.368    0.286     0.196 
6     5 0.411    0.370     0.244 

# Power curves
ggplot(power_sim, aes(x = theta)) +
  geom_line(aes(y = prob_phi1, color = "phi_1 (UMP)")) +
  geom_line(aes(y = prob_phi2, color = "phi_2 (Mean)")) +
  facet_wrap(~n, labeller = label_both) +
  geom_hline(yintercept = 0.05, linetype = 2) +
  scale_y_continuous(breaks = c(0.05, seq(0.25, 1, by = 0.25))) +
  theme_classic() +
  labs(x = expression(theta), y = "P(reject null)", color = "")

• We see that phi_1 is in-fact UMP than phi_2 regardless of sample size and theta value.