R/MVcov.R

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

##' The function MVcov computes the Covariance of a multivariate data set x and the covariance or correlation between x(i) and x(i+b). Covariance γ(q) = Cov(Xi,Xi+q) of a Data Set where X is a multivariate data set
##'
##' The input x must be a matrix. MVcov() computes covariance matrix between X(i) and X(i+q),where X(i) and X(i+q) are multivariate process observations obtained at times i and i + q where q =0,1,2.....bmax. Covariance γ(q) = Cov(Xi,Xi+q), where X is a p-dimensinal vector and γ(q) will be a p by p matrix. The default value of bmax is 10, which means the output will be [γ(0),γ(1),.....γ(10)] where γ(q) = Cov(X(i),X(i+q)) a p by p matrix, so the output will be a p by 11*p matrix.
##' @title Correlation, Variance and Covariance of a Multivariate Data Set
##' @param x Multivariate data
##' @param bmax γ(q) = Cov(Xi,Xi+q), the maximun value of q
##' @return Multivariate Covariance Matrix [γ(0),γ(1),.....γ(bmax)] where γ(q) = Cov(Xi,Xi+q) a p by p matrixs
##' @author Xiulin Xie
##' @export
##' @examples
##' MVcov(matrix(rnorm(900,0,1),nrow = 3))$cov_vector
MVcov <- function(x, bmax = 10) {
    avg <- rowMeans(x)
    p <- nrow(x)
    cov_vector = matrix(rep(0, p * p * (bmax + 1)), nrow = nrow(x))
    sigma_ <- matrix(rep(0, p * p * (bmax + 1)), nrow = nrow(x))
    for (i in 1:(bmax + 1)) {
        l = ((i - 1) * p + 1)
        cov_vector[, l:(l + p - 1)] <- (x[, 1:(ncol(x) - i + 1)] - rowMeans(x)) %*% t(x[, i:ncol(x)] - rowMeans(x))/(ncol(x) - 
            i + 1)
        l2 <- ncol(sigma_)
        sigma_[, (l2 - l - (p - 2)):(l2 - l + 1)] <- t(cov_vector[, l:(l + p - 1)])
    }
    return(list(cov_vector = cov_vector, sigma_ = sigma_))
}