9.1

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(purrrfect)

Attaching package: 'purrrfect'

The following objects are masked from 'package:base':

    replicate, tabulate

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

N <- 10000

(sim_results <- parameters(~mu, ~n, 
            c(-2,-1,0,1,2),
            c(5,10,20,40,80))%>%
    add_trials(N) %>%
    mutate(true_theta = mu^2) %>%
    mutate(Ysample = pmap(list(mu,n), .f = \(mu,n) rnorm(n, mean = mu,sd = 1)))%>%
    mutate(
      Ybar1 = map_dbl(Ysample, mean),
      Ybar2 = map_dbl(Ysample, \(y) mean(y^2)),
      estimator1 = Ybar1^2 - (1/n),
      estimator2 = Ybar2 - 1
    ) %>%
    group_by(mu,n) %>%
    summarise(mean_est1 = mean(estimator1),
              mean_est2 = mean(estimator2),
              var_est1 = var(estimator1),
              var_est2 = var(estimator2),
              .groups = "drop")
  
)

# A tibble: 25 × 6
      mu     n mean_est1 mean_est2 var_est1 var_est2
   <dbl> <dbl>     <dbl>     <dbl>    <dbl>    <dbl>
 1    -2     5     3.97      3.97    3.28     3.61  
 2    -2    10     3.99      3.99    1.59     1.77  
 3    -2    20     4.01      4.01    0.826    0.927 
 4    -2    40     4.00      4.00    0.400    0.449 
 5    -2    80     4.00      4.00    0.207    0.231 
 6    -1     5     1.01      0.993   0.894    1.20  
 7    -1    10     0.999     0.994   0.410    0.593 
 8    -1    20     1.01      1.01    0.206    0.298 
 9    -1    40     1.00      1.00    0.104    0.152 
10    -1    80     0.997     0.999   0.0489   0.0734
# ℹ 15 more rows

ggplot(sim_results, aes(x = mu))+
  geom_point(aes(y = mean_est1, color = "Estimator 1"))+
  geom_point(aes(y = mean_est2, color = "Estimator 2"))+
  geom_line(aes(y = mu^2, color = "True Theta"))+
  facet_wrap(~n)+
  labs(title = "Unbiasedness Check",
       y = "Estimate of Theta")+
  theme_minimal()

ggplot(sim_results, aes(x = n)) + 
  geom_line(aes(y = var_est1, color = "Estimator 1"))+
  geom_line(aes(y = var_est2, color = "Estimator 2"))+
  facet_wrap(~mu) + 
  labs(title = "Variance comparison",
       y = "Variance")+
  theme_minimal()

You can add options to executable code like this

[1] 4

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