R/update_cov.R

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

##' The function updates the variance and covariance γ(q) = Cov(Xi,Xi+q) of a dataset when we have a new data
##'
##' update_cov() can be used to update the variance and covariance of a large data set without recalculating the whole data set to make compuation easier by a recursive formula. The function can be used to update univariate or multivariate data set.
##' @title Update the Covariance of a dataset
##' @param oldx The previous dataset
##' @param old_mean the mean of previous dataset
##' @param old_cov the covariance of previous dataset
##' @param newx value of your new data
##' @return new covariance γ(q) = Cov(Xi,Xi+q)
##' @author Xiulin Xie
##' @export
##' @examples
##' update_cov(c(1:13),mean(c(1:13)),Vcov(c(1:13)),14)

update_cov <- function(oldx, old_mean, old_cov, newx) {
    if (length(newx) == 1) {
        new_mean <- update_mean(old_mean, length(oldx), newx)
        bmax <- length(old_cov) - 1
        m <- length(oldx)
        cov_vector <- rep(0, (bmax + 1))
        for (s in 0:(bmax)) {
            u = s + 1
            if (s == 0) {
                cov_vector[u] <- 1/(m + 1 - s) * (newx - new_mean) * (newx - new_mean) + (m - s)/(m + 1 - 
                  s) * (old_cov[1] + (old_mean - new_mean)^2)
            } else {
                cov_vector[u] <- 1/(m + 1 - s) * (newx - new_mean) * (oldx[m - s + 1] - new_mean) + (m - s)/(m + 
                  1 - s) * (old_cov[u] + (old_mean - new_mean)^2) + (old_mean - new_mean)/(m - s + 1) * (2 * 
                  s * old_mean - sum(oldx[1:s]) - sum(oldx[(m - s + 1):m]))
            }
        }
    }
    
    if (length(newx) != 1) {
        new_mean <- update_mean(old_mean, ncol(oldx), newx)
        bmax <- (ncol(old_cov) - nrow(old_cov))/nrow(old_cov)
        m <- ncol(oldx)
        x <- oldx
        p <- nrow(x)
        cov_vector = matrix(rep(0, p * p * (bmax + 1)), nrow = nrow(x))
        for (s in 0:(bmax)) {
            u = s + 1
            l = ((u - 1) * p + 1)
            if (s == 0) {
                cov_vector[, l:(l + p - 1)] <- 1/(m + 1 - s) * (newx - new_mean) %*% t(newx - new_mean) + 
                  (m - s)/(m + 1 - s) * (old_cov[, l:(l + p - 1)] + (old_mean - new_mean) %*% t(old_mean - 
                    new_mean))
            }
            
            if (s == 1) {
                cov_vector[, l:(l + p - 1)] <- 1/(m + 1 - s) * (newx - new_mean) %*% t(oldx[, (m - s + 1)] - 
                  new_mean) + (m - s)/(m + 1 - s) * (old_cov[, l:(l + p - 1)] + (old_mean - new_mean) %*% 
                  t(old_mean - new_mean)) + matrix((old_mean - new_mean)/(m - s + 1), nrow = p) %*% t(2 * 
                  s * old_mean - (oldx[, 1:s]) - (oldx[, (m - s + 1):m]))
            }
            
            if (s > 1) {
                
                cov_vector[, l:(l + p - 1)] <- 1/(m + 1 - s) * (newx - new_mean) %*% t(oldx[, (m - s + 1)] - 
                  new_mean) + (m - s)/(m + 1 - s) * (old_cov[, l:(l + p - 1)] + (old_mean - new_mean) %*% 
                  t(old_mean - new_mean)) + matrix((old_mean - new_mean)/(m - s + 1), nrow = p) %*% t(2 * 
                  s * old_mean - rowSums(oldx[, 1:s]) - rowSums(oldx[, (m - s + 1):m]))
            }
        }
    }
    return(cov_vector)
}