R/biloop_transform.R

In dcorBSS: Distance-Correlation Based Methods for Blind Source Separation and Dependence Analysis

#' Biloop Transformation
#'
#' Apply the biloop transformation to a numeric vector or matrix.
#'
#' The transformation maps each univariate variable to a two-dimensional
#' nonlinear embedding \eqn{x \mapsto (u(x), v(x))} based on two tangential circles. 
#' The data may first be robustly standardized (advised).
#' For multivariate input, the transformation is applied columnwise and the
#' transformed columns are concatenated.
#'
#' @param X A numeric vector or matrix. Rows correspond to observations and
#'   columns to variables.
#' @param c1 Positive scaling constant used before the hyperbolic tangent
#'   transformation. See details.
#' @param c2 Positive scaling constant applied to the first transformed
#'   coordinate. See details.
#' @param do_scale Logical; if \code{TRUE}, standardize each variable before
#'   applying the biloop transformation.
#' @param location Optional centering vector used when `do_scale = TRUE`.
#'   Defaults to [stats::median()]. Ignored otherwise.
#' @param scale Optional scaling vector used when `do_scale = TRUE`.
#'   Defaults to [stats::mad()]. Ignored otherwise.
#'
#' @return A numeric matrix with \code{nrow(X)} rows and \code{2 * ncol(X)}
#'   columns. If \code{X} is a vector, the result has two columns.
#'
#' @details
#' 
#' \eqn{u(x)} and \eqn{v(x)} are defined as:
#'
#' \eqn{
#' u(x) =
#' \begin{cases}
#' c_2(1 + \cos(2\pi \tanh(x/c_1) + \pi)), & x \ge 0 \\
#' -c_2(1 + \cos(2\pi \tanh(x/c_1) - \pi)), & x < 0
#' \end{cases}
#' }
#'
#' \eqn{v(x) = \sin(2\pi \tanh(x/c_1))}
#' 
#' By default, no scaling is performed. If \code{do_scale = TRUE} and neither
#' \code{location} nor \code{scale} is supplied, the columns are centered using
#' [stats::median()] and scaled using [stats::mad()]. Alternatively,
#' user-supplied location and scale vectors can be provided.
#'
#' The biloop transformation can be useful when computing dependence measures on
#' heavy-tailed or contaminated data.
#'
#' @references
#' Leyder, S., Raymaekers, J., and Rousseeuw, P. J. (2026).
#' Robust Distance Covariance.
#' \emph{International Statistical Review}, 94, 1--25.
#' <doi:10.1111/insr.70005>
#'
#' @examples
#' x <- rnorm(10)
#' biloop_transform(x)
#'
#' X <- cbind(rnorm(10), rt(10, df = 3))
#' colnames(X) <- c("A","B")
#' biloop_transform(X, do_scale = TRUE)
#'
#' # Robust scaling using median and Qn from robustbase
#' if (requireNamespace("robustbase", quietly = TRUE)) {
#'   loc <- apply(X, 2, median)
#'   scl <- apply(X, 2, robustbase::Qn)
#'   biloop_transform(X, do_scale = TRUE, location = loc, scale = scl)
#' }
#'
#' # Robust distance correlation after the biloop transformation
#' if (requireNamespace("robustbase", quietly = TRUE) &&
#'     requireNamespace("energy", quietly = TRUE)) {
#'   x <- rt(100, df = 3)
#'   y <- x^2 + rnorm(100, sd = 0.2)
#'   loc_x <- median(x)
#'   scl_x <- robustbase::Qn(x)
#'   loc_y <- median(y)
#'   scl_y <- robustbase::Qn(y)
#'
#'   xb <- biloop_transform(x, do_scale = TRUE, location = loc_x, scale = scl_x)
#'   yb <- biloop_transform(y, do_scale = TRUE, location = loc_y, scale = scl_y)
#'   energy::dcor(xb, yb)
#' }
#'
#' @export
biloop_transform <- function(X,
                             c1 = 4,
                             c2 = 4,
                             do_scale = FALSE,
                             location = NULL,
                             scale = NULL) {
  X <- .as_matrix_input(X)
  
  if (do_scale) {
    if (is.null(location) && is.null(scale)) {
      location <- apply(X, 2, stats::median)
      scale <- apply(X, 2, stats::mad)
      .check_scale_positive(scale)
      X <- scale(X, center = location, scale = scale)
    } else {
      if (is.null(location)) {
        location <- apply(X, 2, stats::median)
      }
      if (is.null(scale)) {
        scale <- apply(X, 2, stats::mad)
      }
      .check_scale_positive(scale)
      X <- scale(X, center = location, scale = scale)
    }
  }
  
  n <- nrow(X)
  p <- ncol(X)
  out <- matrix(0, nrow = n, ncol = 2L * p)
  
  for (i in seq_len(p)) {
    out[, (2L * i - 1L):(2L * i)] <- .biloop_transform_univariate_scaled(
      X[, i],
      c1 = c1,
      c2 = c2
    )
  }
  
  base_names <- colnames(X)
  if (is.null(base_names)) {
    base_names <- paste0("V", seq_len(p))
  }
  
  colnames(out) <- as.vector(rbind(
    paste0(base_names, "_biloop1"),
    paste0(base_names, "_biloop2")
  ))
  
  out
}