R/multi_online_monitor.R

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

##' The function estimats the correlation structure nonparametrically from an IC dataset, and decorrelating the original multivariate process observations. After data decorrelation, a multivariate nonparametric CUSUM chart based on data categorization was applied to monitor the new data set to detect if there is any mean shift occurs, and give us a signal as soon as possible.
##'
##' Early MSPC methods require the assumptions that multivariate process observations are independent at different time points and the in-control (IC) process distribution is multivariate normal or another parametric distribution. In practice, the assumed parametric distribution and independence would be hardly valid and consequently the traditional MSPC charts would be unreliable to use for many MSPC applications. multi_online_monitor() applied a methods focuses mainly on online monitoring of p-dimensional processes, where p > 0 is an interger. Let the sequence of process observations under sequential monitoring be {Xn = (Xn1,Xn2,...,Xnp)′,n ≥ 1}. One assumption needed by the proposed method is that the serial data correlation in the observed data is stationary when the process is IC, which implies that the covariance matrix γ(s) = Cov(Xi,Xi+s), for any i, depends only on s in such cases. This assumption is often reasonable because it is believed that the IC process distribution, including the IC serial data correlation, does not change over time in many applications, including the production lines in manufacturing industries. Regarding the serial data correlation of an IC process, it is also assumed that there is an integer w ≥ 1 such that γ(s) = 0 when s > w. This assumption implies that the correlation between two observations would disappear if the related observation times are far away, which should be (approximately) true in many applications. The IC process distribution is assumed to be unknown, and the process observations are multivariate and serially correlated. The output of the function multi_online_monitor() is the index of the obsevation in the new data that gives a out of control signal. For example, if the output is 20, the result indicates that at X20, the function detects a mean shift.
##' @title Online Monitoring of Serially Correlated Multivariate Data
##' @param x the in control Multivariate data
##' @param xx the Multivariate data needed to monitor
##' @param h the control limit
##' @param k the allowance value
##' @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 index of the obsevation in the new data that gives a mean shift signal
##' @author Xiulin Xie
##' @import stats
##' @import matrixcalc
##' @export
##' @examples
##' multi_online_monitor(matrix(rnorm(900,0,1),nrow = 3),matrix(rnorm(900,2,1),nrow = 3))
multi_online_monitor <- function(x, xx, h = 20, k = 0.01, bmax = 10) {
    summarytable <- function(x) {
        nx <- decor_multi(x)
        med <- apply(nx, 1, median)
        f <- matrix(as.numeric(nx > med), nrow = nrow(x))
        data <- as.data.frame(t(f))
        table1 <- table(data)
        table1.df <- as.data.frame(as.table(table1))
        p <- nrow(x)
        table1.df[, -(p + 1)] <- lapply(table1.df[, -(p + 1)], relevel, ref = 1)
        result <- cbind(table1.df, expected = table1.df$Freq/sum(table1.df$Freq))
        return(result)
    }
    xd1 <- decor_multi(x, bmax)
    mu <- rowMeans(x)
    cov_vector <- MVcov(x)$cov_vector
    sigma_ <- MVcov(x)$sigma_
    result <- summarytable(x)
    cov_matrix <- MVcovmatrix(x)
    p <- nrow(x)
    n = ncol(xx)
    C = 0
    j = 1
    bmax = 10
    Sobs <- rep(0, 2^p)
    Sexp <- rep(0, 2^p)
    T_ = 0
    while (C <= h & j <= n) {
        if (j == 1 | T_ == 0) {
            obs = Nhalf_power(cov_vector[, 1:p]) %*% (xx[, j] - mu)
        } else {
            l = ncol(sigma_) - (T_ * p - 1) - p
            l2 <- T_ * p
            a <- cov_vector[, 1:p] - sigma_[, l:(ncol(sigma_) - p)] %*% matrix.inverse(cov_matrix[1:(l2), 
                1:(l2)]) %*% t(sigma_[, l:(ncol(sigma_) - p)])
            e <- matrix((xx[, (j - T_):(j - 1)] - mu), nrow = T_ * nrow(x))
            obs <- Nhalf_power(a) %*% (xx[, j] - mu - sigma_[, l:(ncol(sigma_) - p)] %*% matrix.inverse(cov_matrix[1:(l2), 
                1:(l2)]) %*% e)
        }
        
        p <- nrow(x)
        
        
        Y <- rep(0, (2^p))
        med <- apply(xd1, 1, median)
        out <- as.numeric(obs > med)
        all_y <- matrix(as.numeric(unlist(result[, 1:p])) - 1, ncol = p)
        for (z in 1:length(Y)) {
            if (identical(out, all_y[z, ])) {
                Y[z] <- 1
            }
        }
        
        f <- result$expect
        D <- t((Sobs - Sexp + Y - f)) %*% matrix.inverse(diag(Sexp + f)) %*% (Sobs - Sexp + Y - f)
        C <- max(0, (D - k))
        x <- cbind(x, xx[j])
        mu <- rowMeans(x)
        cov_vector <- MVcov(x)$cov_vector
        sigma_ <- MVcov(x)$sigma_
        cov_matrix <- MVcovmatrix(x)
        if (j == 1 & C != 0) {
            T_ = 1
        } else {
            T_ = min((T_ + 1), bmax)
        }
        
        
        if (T_ != 0 & D > k) {
            l = (D - k)/D
            Sobs <- c(Sobs + Y) * c(l)
            Sexp <- c(Sexp + f) * c(l)
            
            
        } else {
            Sobs <- rep(0, 2^p)
            Sexp <- rep(0, 2^p)
        }
        j <- j + 1
        
    }
    if (j < n) {
        return(j)
    } else {
        print("No mean shift detected")
    }
}