sim_method.R

library(gamlss.dist)
library(semTools)

johnsonSU <- function(num_n, num_x, mean, std, skew, kurt, corr) {
  corr_matrix <- matrix(corr, nrow = num_x, ncol = num_x)
  diag(corr_matrix) <- 1.0
  
  # Check is Normal distribution skew = 0, kurt = 0
  is_normal <- abs(skew) < 1e-4 && abs(kurt) < 1e-4
  
  # Multivariate Standard Normal
  mu_z <- rep(0, num_x)
  Z_data <- MASS::mvrnorm(n = num_n, mu = mu_z, Sigma = corr_matrix)
  
  sim_data <- matrix(0, nrow = num_n, ncol = num_x)
  epsilon <- 1e-6 
  
  if (!is_normal) {
    mapped_nu <- skew * 1.5 
    mapped_tau <- max(0.5, 3.0 / (abs(kurt) + 1.0)) 
  }
  
  # NORTA Transform Inverse CDF
  for (j in 1:num_x) {
    # Transform Z to U (Uniform distribution)
    U <- pnorm(Z_data[, j]) 
    U <- pmax(pmin(U, 1 - epsilon), epsilon) 
    
    if (is_normal) {
      raw_data <- stats::qnorm(U)
    } else {
      raw_data <- gamlss.dist::qJSU(U, mu = 0, sigma = 1, nu = mapped_nu, tau = mapped_tau)
    }
    sim_data[, j] <- ((raw_data - base::mean(raw_data)) / stats::sd(raw_data)) * std + mean
  }
  
  df_sim <- as.data.frame(sim_data)
  colnames(df_sim) <- paste0("x", 1:num_x)
  
  return(df_sim)
}

mvrnonnormVM <- function(num_n, num_x, mean, std, skew, kurt, corr) {
  # mvrnonnorm: Generate Non-normal Data using Vale and Maurelli (1983) method
  # Create matrix for multivariate data set (Multiple x same parameter)
  mean_mat <- rep(mean, num_x)
  std_mat  <- rep(std,  num_x)
  skew_mat <- rep(skew, num_x)
  kurt_mat <- rep(kurt, num_x)
  
  corr_matrix <- matrix(corr, nrow = num_x, ncol = num_x)
  diag(corr_matrix) <- 1.0
  cov_matrix <- outer(std_mat, std_mat) * corr_matrix
  
  # Start process to simulate data
  # set.seed(1) 
  simulated_data <- semTools::mvrnonnorm(
    n        = num_n, 
    mu       = mean_mat, 
    Sigma    = cov_matrix, 
    skewness = skew_mat, 
    kurtosis = kurt_mat
  )
  df_sim <- as.data.frame(simulated_data)
  # colnames(df_sim) <- gsub("V", "x", colnames(df_sim))
  
  return(df_sim)
}