R/normBoxCox.R

In fsdaR: Robust Data Analysis Through Monitoring and Dynamic Visualization

######
##  VT::18.05.2026
##
##
##  roxygen2::roxygenise("C:/users/valen/onedrive/myrepo/R/fsdaR", load_code=roxygen2:::load_installed)
##
#'  Computes (normalized) Box-Cox transformation
#'
#' @description Computes (normalized) Box-Cox transformation
#' 
#' @details When \eqn{\lambda \ne 0} if \code{Jacobian=TRUE},
#'  \deqn{ytra = (y^\lambda-1)/ (G^{(\lambda-1)} \lambda)}
#'  else if \code{Jacobian=FALSE}
#'  \deqn{ytra = (y^\lambda-1)/ \lambda}
#'  where G is the geometric mean of the observations.
#'
#'  When \eqn{\lambda = 0} if \code{Jacobian=TRUE},
#'  \deqn{ytra = G log(y)}
#'  else if \code{Jacobian=FALSE}
#'  \deqn{ytra = log(y)}
#'  where G is the geometric mean of the observations.
#'
#' @param X The data matrix: \code{n} observations and \code{p} variables. 
#'  The rows of X represent observations, and the columns represent variables.
#'  Missing values (NA's) and infinite values (Inf's) are allowed, since observations (rows) 
#'  with missing or infinite values will automatically be excluded from the computations.
#'
#' @param Col2Tra Which variables to transform. An integer vector of length \code{k} 
#'  specifying the variables which must be transformed. If it is missing and 
#'  \code{length(la)=p} all variables are transformed
#'
#' @param la Transformation parameters. A vector of lenghth \code{k} containing set 
#'  of transformation parameters for the \code{k} \code{Col2Tra} variables.
#'
#' @param bsb Units to be used in the computation of the Jacobian, a vector of length \code{m} or 
#'  or a logical vector of length \code{n}. The default value is \code{bsb=1:n}, i.e. 
#'  all units are used to compute the Jacobian. Note that this option takes effect 
#'  only if \code{Jacobian=TRUE}.
#'
#' @param inverse Wheather to return the inverse transformation. The default is \code{inverse=FALSE}.
#'
#' @param Jacobian Requested Jacobian of transformed values. If \code{Jacobian=TRUE} 
#'  the default, the transformation is normalized to have Jacobian equal to 1. 
#'  Note that this optional argument is ignored if \code{inverse=TRUE}.
#'
#' @param trace Whether to print intermediate results. Default is \code{trace=FALSE}.
#'
#' @param ... potential further arguments passed to lower level functions.
#'
#' @return The transformed data matrix. 
#'
#' @references
#'   Box, G.E.P. and Cox, D.R. (1964), An analysis of transformations (with
#'   Discussion), "Journal of the Royal Statistical Society Series B", 
#'   Vol. 26, pp. 211-252.
#'
#' @examples
#'
#'  \dontrun{
#'  ##  Mussels data.
#'  data(mussels)
#'  
#'  X <- mussels
#'  la <- c(0.5, 0, 0.5, 0, 0)
#'  
#'  ## Transform all columns of the matrix X according to the values of la
#'  (Y <- normBoxCox(X, la=la))
#'  
#'  ##  Check the inverse tranformation.
#'  data(mussels)
#'  X <- mussels
#'  la <- c(0.5, 0, 0.5, 0, 0)
#'  
#'  ## Transform all columns of matrix Y according to the values of la
#'  Ytra <- normBoxCox(X, la=la, Jacobian=FALSE)
#'  Ychk <- normBoxCox(Ytra, la=la, inverse=TRUE)
#'  
#'  print(max(max(abs(X-Ychk))))
#'  
#'  ## Comparison between Box-Cox and Yeo-Johnson transformation.
#'  y <- seq(from=-2, to=2, by=0.1)
#'  n <- length(y)
#'  la <- seq(from=-1, to=3, by=1) 
#'  nla <- length(la)
#'  YtraYJ <- matrix(0, nrow=n, ncol=nla)
#'  YtraBC <- matrix(NA, nrow=n, ncol=nla)
#'  
#'  posy <- y>0
#'  for(j in 1:nla) {
#'    YtraYJ[,j] <- normYJ(y, 1, la[j], Jacobian=FALSE)
#'    YtraBC[posy, j] <- normBoxCox(y[posy], 1, la[j],Jacobian=FALSE)
#'  }
#'  
#'  oldpar <- par(mfrow=c(1,2))
#'  plot(y, YtraYJ[,1], type="n", xlab="Original values", 
#'      ylab="Transformed values", main="Yeo-Johnson transformation")
#'  for(j in 1:nla)
#'      lines(y, YtraYJ[,j], col=j)
#'  for(j in 1:nla) {
#'      text(y[1], YtraYJ[1,j], paste0("la=", la[j]))
#'  }
#'      
#'  plot(y, YtraBC[,1], type="n", xlab="Original values", 
#'      ylab="Transformed values", main="Box-Cox transformation")
#'  for(j in 1:nla)
#'      lines(y, YtraBC[,j], col=j)
#'  for(j in 1:nla) {
#'      text(y[16], YtraBC[22,j], paste0("la=", la[j]))
#'  }
#'  
#'  par(oldpar)
#'
#'  }
#'
#' @export
#' @author FSDA team, \email{valentin.todorov@@chello.at}

normBoxCox <- function(X, Col2Tra, la, Jacobian=TRUE, inverse=FALSE, bsb, trace=FALSE, ...) {
    if(missing(X))
        stop("Input data matrix is missing")

    if(is.data.frame(X))
      X <- data.matrix(X)
    else if(!is.matrix(X))
      X <- matrix(X, length(X), 1,
                  dimnames = list(names(X), deparse(substitute(X))))
    if(!is.numeric(X)) stop("X is not a numeric matrix")

    storage.mode(X) <- "double"

    dx <- dim(X)
    xn <- (dnx <- dimnames(X))[[2]]
    xn <- if(!is.null(xn))
        xn
    else if (dx[2] > 1)
        paste("X", 1:dx[2], sep = "")
    else if(dx[2])
        "X"
    dimnames(X) <- list(dnx[[1]], xn)

    n <- nrow(X)
    p <- ncol(X)

    if(missing(la))
        stop("Vector 'la' specifying how to transforme the variables is missing!")

    if(missing(Col2Tra)) {
        if(length(la) == p)
            Col2Tra <- 1:p
        else
            stop("Vector 'Col2Tra' specifying which variables to transform is missing!")
    }
    
    control <- list()
    control$Jacobian <- ifelse(Jacobian, 1, 0)
    control$inverse <- ifelse(inverse, 1, 0)
    if(!missing(bsb))
        control$bsb <- bsb

    parlist = c(.jarray(X, dispatch=TRUE), .jarray(Col2Tra, dispatch=TRUE), .jarray(la, dispatch=TRUE))
    paramNames <- names(control)
    if(trace)
        print(control)

    if(length(paramNames) > 0) {
        for(i in 1:length(paramNames)) {
            paramName = paramNames[i]
            paramValue = control[[i]]

            matlabValue = rType2MatlabType(paramName, paramValue)
            parlist = c(parlist, .jnew("java/lang/String", paramName), matlabValue)
        }
    }

    out <- callFsdaFunction("normBoxCox", "[Ljava/lang/Object;", 1, parlist)
    if(is.null(out))
        return(NULL)

##    arr1 = .jcast(out[[1]], "com/mathworks/toolbox/javabuilder/MWStructArray")
##    arr = .jnew("org/jrc/ipsc/globesec/sitaf/fsda/FsdaMWStructArray", arr1)

##    if(trace) {
##        cat("\nReturning from MATLAB normBoxCox().  Fields returned by MATLAB: \n")
##        print(arr$fieldNames())
##    }

##    Xtra <- if(as.integer(arr$hasField("Ytra", as.integer(1))) != 1) NULL
##            else as.matrix(.jevalArray(arr$get("Ytra", as.integer(1)), "[[D", simplify = TRUE))

    arr <- out[[1]]$toDoubleArray()
    Xtra <- .jevalArray(arr, simplify=TRUE)

    freeMatlabResources(out)

    return (Xtra)
}