## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 6.5,
  fig.height = 3.8,
  fig.align = "center",
  dpi = 110
)
options(digits = 4)
set.seed(1)
has_ggplot2 <- requireNamespace("ggplot2", quietly = TRUE)
# vignette code runs with at most 2 workers (the limit R CMD check enforces)
sim_cores <- 2L

## ----load---------------------------------------------------------------------
library(mcstatsim)

## ----first-sim----------------------------------------------------------------
sim_mean <- function(n, mu) {
  x <- rnorm(n, mean = mu)
  data.frame(estimate = mean(x))
}

## ----first-grid---------------------------------------------------------------
grid <- expand.grid(n = c(10, 50), mu = c(0, 1))
grid

## ----first-run----------------------------------------------------------------
res <- runsim(
  n = 1000,
  grid_params = grid,
  sim_func = sim_mean,
  seed = 101,
  num_cores = sim_cores,
  show_progress = FALSE
)
head(res)

## ----first-summary------------------------------------------------------------
summarise_sim(
  res,
  by = c("n", "mu"),
  measures = list(bias = ~ calc_bias(estimate, mu[1]))
)

## ----contract-----------------------------------------------------------------
sim_ci <- function(n, mu) {
  x <- rnorm(n, mean = mu)
  est <- mean(x)
  se <- sd(x) / sqrt(n)
  crit <- qt(0.975, df = n - 1)
  data.frame(
    est   = est,
    lo    = est - crit * se,
    hi    = est + crit * se,
    truth = mu           # constant within a scenario, handy for summarise_sim()
  )
}

## ----repro--------------------------------------------------------------------
a <- runsim(200, grid, sim_mean, seed = 7, num_cores = 1, show_progress = FALSE)
b <- runsim(200, grid, sim_mean, seed = 7, num_cores = 2, show_progress = FALSE)
identical(a$estimate, b$estimate)

## ----repro-attr---------------------------------------------------------------
attr(a, "seed")

## ----on-error-----------------------------------------------------------------
flaky <- function(n, mu) {
  x <- rnorm(n, mu)
  if (mean(x) < -0.5) stop("solver failed to converge")
  data.frame(est = mean(x))
}

res_f <- runsim(300, expand.grid(n = 15, mu = c(-1, 0, 1)), flaky,
                seed = 4, num_cores = sim_cores, show_progress = FALSE,
                on_error = "omit")
nrow(res_f)

## ----failure-summary----------------------------------------------------------
failure_summary(res_f)

## ----checkpoint---------------------------------------------------------------
ckpt <- tempfile(fileext = ".rds")

# first call: imagine this is interrupted partway through
r1 <- runsim(400, grid, sim_mean, seed = 9, save_path = ckpt,
             checkpoint_every = 50, num_cores = sim_cores, show_progress = FALSE)

# a later call with the same settings resumes instead of restarting
r2 <- runsim(400, grid, sim_mean, seed = 9, save_path = ckpt,
             checkpoint_every = 50, num_cores = sim_cores, show_progress = FALSE)

identical(r1$estimate, r2$estimate)
file.remove(ckpt)

## ----summarise----------------------------------------------------------------
res_ci <- runsim(1000, expand.grid(n = c(15, 60), mu = 0), sim_ci,
                 seed = 202, num_cores = sim_cores, show_progress = FALSE)

perf <- summarise_sim(
  res_ci,
  by = "n",
  measures = list(
    bias     = ~ calc_bias(est, truth[1]),
    coverage = ~ calc_coverage(lo, hi, truth[1]),
    width    = ~ calc_width(lo, hi)
  )
)
perf

## ----digits-------------------------------------------------------------------
summarise_sim(res_ci, by = "n",
              measures = list(bias = ~ calc_bias(est, truth[1])),
              digits = NULL)

## ----calc-direct--------------------------------------------------------------
set.seed(3)
theta_hat <- rnorm(500, mean = 2.1, sd = 0.4)
calc_bias(theta_hat, true_param = 2)
calc_rmse(theta_hat, true_param = 2)

## ----check-estimates----------------------------------------------------------
divergent <- function(n, mu) {
  x <- rnorm(n, mu)
  data.frame(est = mean(x), ratio = mu / (mu - 1))   # Inf when mu == 1
}
res_d <- runsim(50, expand.grid(n = 20, mu = c(0, 1, 2)), divergent,
                seed = 5, num_cores = sim_cores, show_progress = FALSE)
check_estimates(res_d, by = "mu")

## ----mcse-target--------------------------------------------------------------
mcse_target(mcse_current = perf$coverage_mcse, n_current = perf$n_sim,
            target = 0.005)

## ----plot-data, include = FALSE-----------------------------------------------
grid2 <- expand.grid(n = c(15, 60), dist = c("normal", "lognormal"),
                     stringsAsFactors = FALSE)

sim_two <- function(n, dist) {
  x <- if (dist == "normal") rnorm(n) else rlnorm(n) - exp(0.5)
  est <- mean(x)
  se  <- sd(x) / sqrt(n)
  crit <- qt(0.975, df = n - 1)
  data.frame(est = est, lo = est - crit * se, hi = est + crit * se, truth = 0)
}

res2 <- runsim(1500, grid2, sim_two, seed = 303, num_cores = sim_cores,
               show_progress = FALSE)
perf2 <- summarise_sim(
  res2, by = c("n", "dist"),
  measures = list(coverage = ~ calc_coverage(lo, hi, truth[1]))
)

## ----plot-convergence, eval = has_ggplot2, fig.height = 4.4-------------------
plot_convergence(res2, "est", by = c("n", "dist"))

## ----plot-zip, eval = has_ggplot2, fig.height = 4-----------------------------
plot_zip(res2[res2$n == 15, ], "lo", "hi", "truth", by = "dist")

## ----plot-performance, eval = has_ggplot2-------------------------------------
plot_performance(perf2, "coverage", x = "n", colour = "dist", ref = 0.95)

## ----case-sim-----------------------------------------------------------------
sim_location <- function(n, eps, estimator) {
  contaminated <- rbinom(n, 1, eps) == 1
  x <- rnorm(n)
  x[contaminated] <- rnorm(sum(contaminated), sd = 8)

  est <- switch(estimator,
    mean    = mean(x),
    trimmed = mean(x, trim = 0.2),
    median  = median(x)
  )
  data.frame(estimate = est, truth = 0)
}

design <- expand.grid(
  n = 40,
  eps = c(0, 0.05, 0.15),
  estimator = c("mean", "trimmed", "median"),
  stringsAsFactors = FALSE
)

## ----case-run-----------------------------------------------------------------
sim_res <- runsim(2000, design, sim_location, seed = 2024,
                  num_cores = sim_cores, show_progress = FALSE)

## ----case-summary-------------------------------------------------------------
perf_loc <- summarise_sim(
  sim_res,
  by = c("estimator", "eps"),
  measures = list(
    bias = ~ calc_bias(estimate, truth[1]),
    rmse = ~ calc_rmse(estimate, truth[1])
  )
)
perf_loc[order(perf_loc$eps, perf_loc$rmse), ]

## ----case-plot, eval = has_ggplot2--------------------------------------------
plot_performance(perf_loc, "rmse", x = "eps", colour = "estimator")