R/online_monitor.R

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

##' online_monitor() estimats the correlation structure nonparametrically from an IC dataset, and decorrelating the original process observations. After data decorrelation, a univariate 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.
##'
##' Traditional statistical process control charts are based on the assumptions that process observations are independent and identically normally distributed when the related process is in-control (IC). online_monitor() applied a general charting scheme for monitoring serially correlated process observations with short-memory data dependence and unknown process distributions. The method focus on Phase II online monitoring of process observations X1,X2,...,Xn, where n ≥ 1 is the current time point during process monitoring. The IC process distribution is assumed to be unknown, and the process observations are serially correlated. The output of the function 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 14, the result indicates that at X14, the function detects a mean shift.
##' @title Online Monitoring of Serially Correlated Data
##' @param x1 the in control data set
##' @param xx the new data needed to monitor
##' @param h the control lomit
##' @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
##' online_monitor(x1 = rnorm(200,0,1),xx = rnorm(100,10,1) ,h=20,k=0.01)
online_monitor <- function(x1, xx, h = 10, k = 0.01, bmax = 10) {
    xd1 <- decor(x1, bmax)
    Q <- quantile(xd1, c(1:9)/10)
    mu <- mean(x1)
    sig = sqrt(var(x1))
    p <- 10
    n = length(xx)
    C = 0
    j = 1
    bmax = 10
    Sobs <- rep(0, p)
    Sexp <- rep(0, p)
    T_ = 0
    while (C <= h & j <= n) {
        if (j == 1 | T_ == 0) {
            obs = (xx[j] - mu)/sig
        }
        if (T_ == 1) {
            sig1 <- cov_matrix[1:T_, (T_ + 1)]
            e <- xx[(j - T_):(j - 1)] - mu
            a <- sig^2 - sig1/sig1_ * sig1
            
            obs <- (xx[j] - mu - sig1/sig1_ * e)/sqrt(a)
        }
        if (T_ > 1) {
            sig1 <- cov_matrix[1:T_, (T_ + 1)]
            e <- xx[(j - T_):(j - 1)] - mu
            a <- sig^2 - t(sig1) %*% matrix.inverse((sig1_)) %*% sig1
            l <- t(sig1) %*% matrix.inverse(sig1_) %*% e
            obs <- (xx[j] - mu - l)/sqrt(a)
        }
        r = 1
        Y <- rep(0, p)
        if (obs > Q[9]) {
            Y[10] = Y[10] + 1
        } else {
            while (obs > Q[r] & r < 9) {
                r = r + 1
            }
            
            Y[r] = Y[r] + 1
        }
        
        
        f <- rep(1/p, p)
        D <- t((Sobs - Sexp + Y - f)) %*% matrix.inverse((diag(Sexp + f))) %*% (Sobs - Sexp + Y - f)
        C <- max(0, (D - k))
        
        xd1 <- append(xd1, obs)
        Q <- quantile(xd1, c(1:(p - 1))/(p))
        x1 <- append(x1, xx[j])
        
        if (j == 1 & C != 0) {
            T_ = 1
        } else {
            T_ = min((T_ + 1), bmax)
        }
        
        mu <- update_mean(mu, length(x1), xx[j])
        sig = sqrt(var(x1))
        if (T_ != 0 & D > k) {
            cov_matrix <- Vcovmatrix(x1, bmax)
            sig1_ <- cov_matrix[1:T_, 1:T_]
            l = (D - k)/D
            Sobs <- c(Sobs + Y) * c(l)
            Sexp <- c(Sexp + f) * c(l)
        } else {
            cov_matrix <- Vcovmatrix(x1, bmax)
            sig1_ <- cov_matrix[1:T_, 1:T_]
            Sobs <- rep(0, p)
            Sexp <- rep(0, p)
        }
        
        j = j + 1
        
    }
    if (j < n) {
        return(j)
    } else {
        print("No mean shift dected")
    }
}