10.3

Quarto

Quarto enables you to weave together content and executable code into a finished document. To learn more about Quarto see https://quarto.org.

Running Code

When you click the Render button a document will be generated that includes both content and the output of embedded code. You can embed code like this:

library(tidyverse)

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

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

Warning: package 'tidyr' was built under R version 4.4.3

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

Warning: package 'purrr' was built under R version 4.4.3

Warning: package 'dplyr' was built under R version 4.4.3

Warning: package 'stringr' was built under R version 4.4.3

Warning: package 'forcats' was built under R version 4.4.3

Warning: package 'lubridate' was built under R version 4.4.3

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   4.0.1     ✔ tibble    3.2.1
✔ lubridate 1.9.4     ✔ tidyr     1.3.1
✔ purrr     1.0.4     
── 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

alpha <- 0.05

(two_var_tests <- parameters(~sigma2_true, ~n,
  seq(1,10,length.out = 100),
  c(5,10,15,20))
  %>%
    mutate(
      crit_chi = qchisq(1 - alpha, df = n),
      crit_z = qnorm(1 - alpha),
     prob_test1 = 
       1 - pchisq(crit_chi / sigma2_true, df = n),
     prob_test2 = 
       1 - pnorm(crit_z / sqrt(sigma2_true))
    )
  )

# A tibble: 400 × 6
   sigma2_true     n crit_chi crit_z prob_test1 prob_test2
         <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
 7        1.09    15     25.0   1.64     0.0860     0.0576
 8        1.09    20     31.4   1.64     0.0919     0.0576
 9        1.18     5     11.1   1.64     0.0953     0.0651
10        1.18    10     18.3   1.64     0.115      0.0651
# ℹ 390 more rows

ggplot(two_var_tests)+
  geom_line(aes(x = sigma2_true, y = prob_test1, col = "Test 1"))+
  geom_line(aes(x = sigma2_true, y = prob_test2, col = "Test 2"))+
  facet_wrap(~n, labeller = label_both)+
  geom_hline(yintercept = 0.05, linetype = 2)+
  theme_classic()+
  labs(
    x = expression(sigma^2),
    y = "P(reject H0)",
    color = ""
  )

You can add options to executable code like this

library(tidyverse)
library(purrrfect)

alpha = 0.05
theta0 = 0.2

(crit_vals <- parameters(~n, c(5,10,15,20)) %>%
  add_trials(10000)%>%
  mutate(
    x = map(n, ~ runif(.x, 0, theta0)),
    T1 = map_dbl(x, ~ sum(log(.x))),
    xbar = map_dbl(x, mean)
  )%>%
  group_by(n) %>%
  summarise(
    crit_phi1 = quantile(T1, 1 - alpha),
    crit_phi2 = quantile(xbar, 1 -alpha),
    .groups = "drop"
  )
)

# A tibble: 4 × 3
      n crit_phi1 crit_phi2
  <dbl>     <dbl>     <dbl>
1     5     -10.0     0.143
2    10     -21.5     0.130
3    15     -33.5     0.124
4    20     -45.5     0.121

(params <- parameters(
  ~theta_true, ~n, 
  seq(0.2,1,length.out = 20),
  c(5,10,15,20)) %>%
  
  add_trials(10000)%>%
  mutate(
    X = map2(n, theta_true, ~ runif(.x,0,.y)),
    T1 = map_dbl(X, ~ sum(log(.x))),
    xbar = map_dbl(X, mean)
  )%>%
  left_join(crit_vals, by = "n")%>%
  mutate(
    reject_phi1 = T1 > crit_phi1,
    reject_phi2 = xbar > crit_phi2
  )%>%
  group_by(n, theta_true) %>%
  summarise(
    prob_phi1 = mean(reject_phi1),
    prob_phi2 = mean(reject_phi2),
    .groups = "drop"
  )
)

# A tibble: 80 × 4
       n theta_true prob_phi1 prob_phi2
   <dbl>      <dbl>     <dbl>     <dbl>
 1     5      0.2      0.0467    0.0457
 2     5      0.242    0.168     0.241 
 3     5      0.284    0.315     0.492 
 4     5      0.326    0.452     0.682 
 5     5      0.368    0.569     0.806 
 6     5      0.411    0.659     0.882 
 7     5      0.453    0.724     0.921 
 8     5      0.495    0.771     0.944 
 9     5      0.537    0.821     0.963 
10     5      0.579    0.852     0.974 
# ℹ 70 more rows

 params %>%
  pivot_longer(
    cols = starts_with("prob"),
    names_to = "test",
    values_to = "prob"
  ) %>%
  mutate(
    test = recode(test,
      prob_phi1 = "phi1: sum log X (UMP)",
      prob_phi2 = "phi2: mean"
    )
  ) %>%
  ggplot(aes(x = theta_true, y = prob, color = test, linetype = test)) +
  geom_line(size = 1) +
  facet_wrap(~n, labeller = label_both) +
  geom_hline(yintercept = 0.05, linetype = 2) +
  theme_classic() +
  labs(
    x = expression(theta[true]),
    y = "P(reject H0)",
    color = "",
    linetype = ""
  )

Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.

The echo: false option disables the printing of code (only output is displayed).