R/simulate_p.R

In BriceonWiley/IntegratedLikelihood.R:

simulate_p <- function(
  N, p = 0.1, phi = 0.1, theta = 0.1, k = NULL, alpha = 0.05, f = 0.1,
  len = 1e4, check = NULL, summary = FALSE, closed_form = FALSE
) {
  # browser()

  # messaging
  print_n <- str_pad(as.character(N), 3, "left", "0")
  print_p <- str_pad(as.character(p), 4, "right", "0")

  if (is.null(k)) message(paste0("N is ", N, ", and p is ", print_p, "."))
  else {
    print_k <- str_pad(as.character(k), 4, "right", "0") %>%
      paste(collapse = ", ")
    message(paste0("N is ", N, ", p is ", print_p, " and k is an element of {", print_k, "}."))
  }


  # producing and checking samples
  samples <- produce_samples(N, p, phi, theta, len, f)
  if(!is.null(check)) samples <- check_zeros(samples, check)


  # producing intervals
  mWalds <- pmap_dfr(samples, ~ mWald_p(..1, ..2, ..3, ..4, ..5, ..6, alpha))
  intLiks <-
    pmap_dfr(
      samples,
      ~ intLik_p(..1, ..2, ..3, ..4, ..5, ..6, alpha, closed_form)
      )
  counts <-
    c(
      "N" = N, "M" = N * 0.9, "n" = N * 0.1, "B" = len,
      "NaN" = sum(is.nan(mWalds$mW_Lower))
    )
  proportions <- c("p" = p, "phi" = phi, "theta" = theta)
  intervals <- bind_cols(mWalds, intLiks)

  if(!is_null(k)) {
    adj_intLiks <-
      expand_grid(k, samples) %>%
      pmap_dfr(~ adj_intLik_p(..2, ..3, ..4, ..5, ..6, ..7, ..1)) %>%
      {
        x <- .;
        map(k, ~ select(x, ends_with(paste0("k=", .x))))
      } %>%
      pmap_dfc(
        .l = list(
          seq(1, length(k) * len, by = len),
          seq(len, length(k) * len, by = len),
          .
        ),
        .f = ~ slice(..3, ..1:..2)
      )

    proportions <- c(proportions, "k" = k)
    intervals <- bind_cols(intervals, adj_intLiks)
  }

  results <- list(
    "counts" = counts,
    "proportions" = proportions,
    "samples" = samples,
    "intervals" = intervals
  )


  # producing final product
  if (summary) {
    coverages <- coverage_p(results, ignore_NaN = TRUE)
    widths <- width_p(results, ignore.NaN = TRUE)
    summarized <-
      bind_rows(coverages, widths) %>%
      select(-N, -p) %>%
      t() %>%
      `colnames<-`(c("Coverage Probability", "Average Width"))
    message(
      paste0(
        "Removed ", results$counts['NaN'],
        " bad intervals in summary calculations."
      )
    )
    list(
      "counts" = results$counts,
      "proportions" = results$proportions,
      "summary" = round(summarized, 4)
    )
  } else results
}