R/generate_simulation.R

In classbound: Visualization for Classification Decision Boundaries

#' Simulate data using Gaussian mixture models
#'
#' @description
#' Generates synthetic classification data using Gaussian mixture models constructed by
#' the `MixSim` package. Useful for creating controlled datasets to explore decision
#' boundary behavior when no real dataset is available.
#'
#' @details
#' Class distributions are randomly generated subject to the `MaxOmega` overlap
#' constraint. The simulation is fully reproducible when `seed` is supplied.
#'
#' ## Test data
#'
#' When `test_ratio > 0`, an additional independent dataset is generated using the
#' **same mixture parameters** but a fresh random draw. This is not a split of the
#' training data: the training set has `n` observations and the test set has
#' `round(n * test_ratio)` independently generated observations. The two sets are
#' statistically independent, so the test set is a fair evaluation sample.
#'
#' @param n Number of training observations to generate.
#' @param K Number of classes.
#' @param p Number of numeric features (dimensions).
#' @param MaxOmega Maximum pairwise overlap between mixture components (0 to 1).
#'   Smaller values create better-separated classes.
#' @param class_names Optional character vector of class labels (length `K`).
#'   Defaults to `"Class 1"`, `"Class 2"`, etc.
#' @param seed Optional integer for reproducibility. The global random seed is
#'   restored after the call.
#' @param noise_ratio Numeric in [0, 1). Proportion of `n` to add as uniform background
#'   noise (randomly labeled). Use this to test classifier robustness to contamination.
#' @param test_ratio Numeric in [0, 1). If greater than 0, generates an additional
#'   independent test dataset of size `round(n * test_ratio)` using the same mixture
#'   parameters. This is not a split of the training data.
#'
#' @return If `test_ratio == 0` (default): a data frame with a `Sim` class column and
#'   `p` feature columns (`X1`, `X2`, ...). If `test_ratio > 0`: a list with `$train`
#'   and `$test` data frames.
#'
#' @examples
#' \donttest{
#' # Generate a 3-class, 2-dimensional training dataset
#' train_df <- simulate_mixsim(n = 200, K = 3, p = 2, MaxOmega = 0.05, seed = 42)
#' head(train_df)
#'
#' # Generate training and independent test data
#' sim <- simulate_mixsim(
#'   n = 200, K = 3, p = 2, MaxOmega = 0.05,
#'   seed = 42, test_ratio = 0.3
#' )
#' nrow(sim$train) # 200
#' nrow(sim$test) # 60 (independently generated, not split from train)
#' }
#' @seealso [simu_n()]
#' @export
simulate_mixsim <- function(n, K, p, MaxOmega, class_names = NULL, seed = NULL, noise_ratio = 0, test_ratio = 0) {
  if (!requireNamespace("MixSim", quietly = TRUE)) {
    stop("Package 'MixSim' must be installed to use simulate_mixsim().", call. = FALSE)
  }

  if (!is.null(seed)) {
    if (exists(".Random.seed", envir = globalenv(), inherits = FALSE)) {
      old_seed <- globalenv()$.Random.seed
      on.exit(assign(".Random.seed", old_seed, envir = globalenv()), add = TRUE)
    } else {
      on.exit(rm(".Random.seed", envir = globalenv()), add = TRUE)
    }
    set.seed(seed)
  }

  if (is.null(class_names)) {
    class_names <- paste0("Class ", 1:K)
  }

  sim_params <- MixSim::MixSim(MaxOmega = MaxOmega, K = K, p = p)

  # Helper to generate and format a single dataset
  generate_df <- function(num_obs) {
    raw <- MixSim::simdataset(n = num_obs, Pi = sim_params$Pi, Mu = sim_params$Mu, S = sim_params$S)
    df <- data.frame(raw$X)
    colnames(df) <- paste0("X", seq_len(p))
    df$Sim <- factor(class_names[raw$id], levels = class_names)
    df <- df[, c("Sim", colnames(df)[-ncol(df)])]

    if (noise_ratio > 0) {
      n_noise <- ceiling(num_obs * noise_ratio)
      noise_df <- data.frame(matrix(
        stats::runif(n_noise * p, min = apply(df[, -1, drop = FALSE], 2, min), max = apply(df[, -1, drop = FALSE], 2, max)),
        ncol = p, byrow = TRUE
      ))
      colnames(noise_df) <- colnames(df)[-1]
      noise_df$Sim <- sample(levels(df$Sim), n_noise, replace = TRUE)
      df <- rbind(df, noise_df[, colnames(df)])
    }
    return(df)
  }

  train_df <- generate_df(n)

  if (test_ratio > 0) {
    test_n <- max(1, round(n * test_ratio))
    test_df <- generate_df(test_n)
    return(list(train = train_df, test = test_df))
  }

  return(train_df)
}

#' Simulate multivariate normal data for multiple classes
#'
#' @description
#' Generates synthetic classification data by drawing independent multivariate normal
#' samples for each class. Useful when you want explicit control over each class's mean
#' and covariance structure.
#'
#' @details
#' ## Test data
#'
#' When `test_ratio > 0`, an additional independent test dataset is generated by
#' drawing fresh samples of size `round(ns * test_ratio)` for each class using the
#' same `means` and `covs`. This is not a split of the training data. The training
#' set has `sum(ns)` observations; the test set has `sum(round(ns * test_ratio))`
#' independently generated observations.
#'
#' @param means A list of numeric vectors, one per class, specifying the class means.
#'   Each vector must have length equal to the number of features.
#' @param covs A list of covariance matrices, one per class.
#' @param ns A numeric vector of sample sizes, one per class.
#' @param class_names Optional character vector of class labels (length equal to
#'   `length(ns)`). Defaults to `"Class 1"`, `"Class 2"`, etc.
#' @param seed Optional integer for reproducibility. The global random seed is
#'   restored after the call.
#' @param noise_ratio Numeric in [0, 1). Proportion of `sum(ns)` to add as uniform
#'   background noise (randomly labeled).
#' @param test_ratio Numeric in [0, 1). If greater than 0, generates an additional
#'   independent test dataset. This is not a split of the training data.
#'
#' @return If `test_ratio == 0` (default): a data frame with a `Sim` class column
#'   and feature columns (`X1`, `X2`, ...). If `test_ratio > 0`: a list with
#'   `$train` and `$test` data frames.
#'
#' @examples
#' \donttest{
#' means <- list(c(0, 0), c(3, 3), c(0, 5))
#' covs <- list(diag(2), diag(2), diag(2))
#' ns <- c(60, 60, 60)
#'
#' train_df <- simu_n(means, covs, ns, seed = 1)
#' head(train_df)
#'
#' # With independent test data
#' sim <- simu_n(means, covs, ns, seed = 1, test_ratio = 0.3)
#' nrow(sim$train) # 180
#' nrow(sim$test) # 54 (independently generated, not split from train)
#' }
#' @seealso [simulate_mixsim()]
#' @export
simu_n <- function(means, covs, ns, class_names = NULL, seed = NULL, noise_ratio = 0, test_ratio = 0) {
  num_classes <- length(ns)

  if (!requireNamespace("MASS", quietly = TRUE)) {
    stop("Package 'MASS' must be installed to use simu_n().", call. = FALSE)
  }

  if (!is.null(seed)) {
    if (exists(".Random.seed", envir = globalenv(), inherits = FALSE)) {
      old_seed <- globalenv()$.Random.seed
      on.exit(assign(".Random.seed", old_seed, envir = globalenv()), add = TRUE)
    } else {
      on.exit(rm(".Random.seed", envir = globalenv()), add = TRUE)
    }
    set.seed(seed)
  }

  if (is.null(class_names)) {
    class_names <- paste0("Class ", 1:num_classes)
  }

  generate_df <- function(sizes) {
    sim_data <- lapply(1:num_classes, function(i) {
      if (sizes[i] == 0) {
        return(NULL)
      }
      sim <- MASS::mvrnorm(n = sizes[i], mu = means[[i]], Sigma = covs[[i]])
      if (sizes[i] == 1) sim <- matrix(sim, nrow = 1)
      df <- data.frame(Sim = class_names[i], sim)
      colnames(df)[-1] <- paste0("X", seq_len(ncol(sim)))
      df
    })

    df <- do.call(rbind, sim_data)
    if (is.null(df)) {
      return(NULL)
    }

    df$Sim <- factor(df$Sim, levels = class_names)

    if (noise_ratio > 0) {
      n_noise <- ceiling(sum(sizes) * noise_ratio)
      p <- ncol(df) - 1
      noise_df <- data.frame(matrix(
        stats::runif(n_noise * p, min = apply(df[, -1, drop = FALSE], 2, min), max = apply(df[, -1, drop = FALSE], 2, max)),
        ncol = p, byrow = TRUE
      ))
      colnames(noise_df) <- colnames(df)[-1]
      noise_df$Sim <- sample(levels(df$Sim), n_noise, replace = TRUE)
      df <- rbind(df, noise_df[, colnames(df)])
    }

    return(df)
  }

  train_df <- generate_df(ns)

  if (test_ratio > 0) {
    test_ns <- pmax(1, round(ns * test_ratio))
    test_df <- generate_df(test_ns)
    return(list(train = train_df, test = test_df))
  }

  return(train_df)
}