R/decor_multi.R

In XiulinXie/SPCmonitor2: Process Control Chart for Monitoring Short-Range Serially Correlated Data

##' The function decorrelates a Multivariate Data Set so that the new Multivariate data will have no correlation or very small correlations
##'
##' The input x must be a matrix. decor_multi() decorrelate a Multivariate Data Set based on Cholesky decomposition. It is assumed that the data observations are covariance stationary and the serial correlation exists only when two observations are within bmax > 0 in their observation indices. More specifically, it is assumed that γ(q) = Cov(Xi,Xi+q) only depends on q when i changes, and γ(q) = 0 when q > bmax, where Xi and Xi+q are two multivariate process observations obtained at times i and i + q. In practice, the autocorrelation between Xi and Xi+q usually decays when q increases. In such cases, γ(q) is small when q is large, and thus a proper value of bmax can be chosen such that γ(q) ≈ 0 when q > bmax. The default value of bmax is 10, which indicates that the assumption is γ(q) = Cov(Xi,Xi+q) = 0 when q is greater than 10
##' @title Decorrelate Multivariate data observations
##' @param x Multivariate data
##' @param bmax The smallest number that any data points have no or small correlations with previous data points which means γ(q) ≈ 0 when q > bmax
##' @return The Multivariate data after decorrelation
##' @author Xiulin Xie
##' @export
##' @import matrixcalc
##' @examples
##' decor(matrix(rnorm(900,0,1),nrow = 3))
decor_multi <- function(x, bmax = 10) {
    cov_vector <- MVcov(x)$cov_vector
    cov_matrix <- MVcovmatrix(x)
    sigma_ <- MVcov(x)$sigma_
    avg <- rowMeans(x)
    nx <- matrix(rep(0, nrow(x) * ncol(x)), nrow = nrow(x))
    for (i in 1:ncol(x)) {
        if (i == 1) {
            nx[, i] <- Nhalf_power(cov_vector[, 1:nrow(x)]) %*% (x[, i] - avg)
            b = 1
        } else {
            l = ncol(sigma_) - (b * nrow(x) - 1) - nrow(x)
            l2 <- b * nrow(x)
            a <- cov_vector[, 1:nrow(x)] - sigma_[, l:(ncol(sigma_) - nrow(x))] %*% matrix.inverse(cov_matrix[1:(l2), 
                1:(l2)]) %*% t(sigma_[, l:(ncol(sigma_) - nrow(x))])
            e <- matrix((x[, (i - b):(i - 1)] - avg), nrow = b * nrow(x))
            nx[, i] <- (Nhalf_power(a) %*% (x[, i] - avg - sigma_[, l:(ncol(sigma_) - nrow(x))] %*% matrix.inverse(cov_matrix[1:(l2), 
                1:(l2)]) %*% e))
            b <- min(b + 1, bmax)
        }
        
    }
    return(nx)
}