R/utils-bootstrap.R

In eFCM: Exponential Factor Copula Model

bootstrap_fcm <- function(
    R           = R,
    s           = s,
    object      = object,
    theta0      = theta0,
    thres       = thres,
    nu          = nu,
    censorL     = censorL,
    control     = control,
    lower       = lower,
    upper       = upper,
    progress    = progress,
    sample_prop = sample_prop,
    sample_ids  = sample_ids,
    parallel = FALSE,
    ncpus = 4,
    mc.set.seed,
    ...){

  data_mat  <- object$data
  coord_mat <- object$coord
  neigh_idx <- object$neigh[[s]] - 1L

  if(is.null(sample_ids)){
    n = nrow(data_mat)
    sample_size <- round(n * sample_prop)
    sample_ids <- replicate(
      n = R,
      expr = sample.int(n, size = sample_size, replace = TRUE),
      simplify = FALSE
    )
  }

  stat_fun <- function(idx) {
    tryCatch({
      d_sub <- as.matrix(data_mat[idx, , drop = FALSE])
      init <- fit_initial(
        coord = coord_mat,
        data  = d_sub,
        neigh = neigh_idx
      )
      out <- optim(
        par     = log(theta0),
        fn      = model_likelihood,
        data_u  = init$data_u,
        coord   = init$coord,
        thres   = thres,
        nu      = nu,
        censorL = censorL,
        lower   = lower,
        upper   = upper,
        control = control
      )
      return(exp(out$par))
    }, error = function(e) {
      warning(e$message)
      return(rep(NA, length(theta0)))
    })
  }
  requireNamespace("pbmcapply", quietly = TRUE)

  if (parallel) {
    if (.Platform$OS.type == "windows") {
      message("pbmclapply falls back to serial execution on Windows (still shows progress).")
    }
    results_list <- pbmcapply::pbmclapply(
      sample_ids, stat_fun,
      mc.cores       = max(1L, as.integer(ncpus)),
      mc.set.seed    = mc.set.seed,
      mc.preschedule = FALSE
    )
  } else {
    results_list <- pbmcapply::pbmclapply(
      sample_ids, stat_fun,
      mc.cores       = 1L,
      mc.set.seed    = mc.set.seed,
      mc.preschedule = FALSE
    )
  }

  results <- do.call(rbind, results_list)
  return(results)
}