Testing a solution to issue 114 (TruncExpFam)

issue-114.R

# Constants
max_si <- 100000L
tolerances <- c(0, 1e-9)

# Load required packages, create testing function
library(parallel)

test_me <- function(si, quiet = TRUE, fam = "chisq", tolerance = 0) {
  set.seed(si)
  results <- c()
  if (!quiet) cat(" =================== seed ", si, "\n")
  for (lg in c(FALSE, TRUE)) {
    for (lt in c(TRUE, FALSE)) {
      for (i in seq_len(3L)) {
        df <- sample(1:10, 1L)
        pt <- runif(i)
        if (lg) pt <- log(pt)
        a <- min(qtrunc(pt, fam, df, lower.tail = lt, log.p = lg) / 2)
        q_trunc <- qtrunc(pt, fam, df, lower.tail = lt, log.p = lg, a = a)
        q_stats <- qchisq(pt, df, lower.tail = lt, log.p = lg)
        for (ii in seq_along(pt)) {
          results <- c(results, q_trunc[ii] + tolerance >= q_stats[ii])
        }
      }
    }
  }
  return(results)
}

# Testing
for (tol in tolerances) {
  message("Testing tolerance ", tol)
  cl <- makeCluster(detectCores() - 1L)
  clusterEvalQ(cl, library("TruncExpFam"))

  message("Testing on seeds 1 through ", max_si)
  out <- parLapply(cl, seq_len(max_si), test_me, tolerance = tol)
  stopCluster(cl)

  cat(length(out), "seeds tested\n")
  message("Proportion of tests passed on each seed")
  print(summary(data.frame("pct passed" = sapply(out, mean))))
}