R/rbmacc.R

In rbcc: Risk-Based Control Charts

#-----------------------------------------------------------------------------#
#                                                                             #
#               RISK-BASED CONTROL CHARTS                                     #
#                                                                             #
#  Written by: Aamir Saghir, Attila I. Katona, Zsolt T. Kosztyan              #
#              Department of Quantitative Methods                             #
#              University of Pannonia, Hungary                                #
#              kzst@gtk.uni-pannon.hu                                         #
#                                                                             #
# Last modified: January 2025                                                 #
#-----------------------------------------------------------------------------#
#' @export
rbmacc <- function(X, UC, C, n=1, w=2, K=3)
{
  if (missing(X))
    stop("data vector/matrix is not specified")
  if (missing(UC))
    stop("Meaurement error vector/matrix is not specified")
  if (missing(C))
  {stop("Cost vector argument is missing")}
  if(!(length(C)==4))
  {stop("Cost should be a vector of length 4!")}
  if(missing(n))
  {n <- 1}
  if(missing(w))
  {w <- 2}
  if(missing(K))
  {K <- 3}
    X <- stats:: na.omit(X)
    UC<- stats:: na.omit(UC)
    n_int <- n*(floor(length(X)/n))
    X <- X[1:n_int]
    UC <- UC[1:n_int]
    x <- matrix(X,ncol=n) #  Data with subgroups
    xrm <- rowMeans(x)
    ma_x <-  pracma::movavg(xrm, w, type="s") # real values of MA statistic
    T11=numeric(length(x))
        for (k in 1:(w-1))  T11[k] <- mean(ma_x)-3*(pracma::std(ma_x)/sqrt(n*k))       # LCL of MA chart
       for (k in w:length(x)) T11[k] <- mean(ma_x)-3*(pracma::std(ma_x)/sqrt(n*w))
    T22=numeric(length(x))
        for (k in 1:(w-1))   T22[k] <- mean(ma_x)+3*(pracma::std(ma_x)/sqrt(n*k))       # UCL of MA chart
        for (k in w:length(x)) T22[k] <- mean(ma_x)+3*(pracma::std(ma_x)/sqrt(n*w))
    T1=T11
    T2=T22
    Y <- X+UC                      # measurement error data matrix
    y <- matrix(Y,ncol=n)
    yrm <- rowMeans(y)
    may <- pracma::movavg(yrm, w, type="s")         # Observed values of MA with measurement errors
    T21=numeric(length(y))
    for (k in 1:(w-1))  T21[k] <- mean(may)-K*(pracma::std(may)/sqrt(n*k))       # LCL of MA chart
    for (k in w:length(y)) T21[k] <- mean(may)-K*(pracma::std(may)/sqrt(n*w))
    T31=numeric(length(y))
    for (k in 1:(w-1))   T31[k] <- mean(may)+K*(pracma::std(may)/sqrt(n*k))       # UCL of MA chart
    for (k in w:length(y)) T31[k] <- mean(may)+K*(pracma::std(may)/sqrt(n*w))
  
    T3 <- T21  #mean(ma_x)-K*(pracma::std(ma_x)/sqrt(n*w))  # LCL of MA chart
    T4 <- T31   #mean(ma_x)+K*(pracma::std(ma_x)/sqrt(n*w))  # UCL of MA chart
    P1 <- ((T1 < ma_x & ma_x < T2) & (T3< may & may<T4))*1  # correct acceptance
    P2 <- ((T1 < ma_x & ma_x < T2) & (T4< may | may<T3))*1 # type I error
    P3 <- ((T2 < ma_x|ma_x < T1) & (T3< may & may<T4))*1 # type II error
    P4 <- ((T2 < ma_x|ma_x < T1) & (T4< may| may<T3))*1 # correct rejecting
  C0 <- sum(P1)*C[1]+sum(P2)*C[2]+sum(P3)*C[3]+sum(P4)*C[4] # calculation of total cost during the process
  C1 <- sum(P1)*C[1]    # total cost related to decision 1 (c11)
  C2 <- sum(P2)*C[2]    # total cost related to decision 2 (c10)
  C3 <- sum(P3)*C[3]    # total cost related to decision 3 (c01)
  C4 <- sum(P4)*C[4]    # total cost related to decision 4 (c00)
 output <- list(cost0=C0, cost1= C1, cost2= C2, cost3= C3, cost4= C4, LCLx=T1, UCLx=T2, LCLy=T3, UCLy=T4, real=ma_x, Observed= may)
 class(output) <- "rbcc"
 return(output)
}