R/ecoDesign.R In edcc: Economic Design of Control Charts

Documented in contour.edccechCusumechEwmaechXbarecoCusumecoEwmaecoXbarupdate.edcc

##' Calculate the optimum parameters, n(sample size), h(sampling
##' interval) and L(number of s.d. from control limits to center line)
##' for Economic Design of the X-bar control chart .
##'
##' When parameter \code{par} is specified, optimization algorithms
##' would be used as default. \code{par} can be specified as:
##' \code{par = c(h, L)} where \code{h} and \code{L} are the initial
##' values of smapling interval and control limit when \code{n} is
##' specified; or \code{par = c(h, L, n)}. Good inital values may lead
##' to good optimum results.
##'
##' When parameters \code{h}, \code{L}, \code{n} are all undefined,
##' \code{ecoXbar} function will try to find the global optimum point
##' to minimize the ECH (Expected Cost per Hour) using optimization
##' algorithms (\code{optim} function), but in this case \code{n}
##' would not be integer. It is usually helpful for the experimenter
##' to find the region where the optimum point may exist quickly. When
##' \code{h} and \code{L} are undefined but n is given as an integer
##' vector, \code{ecoXbar} function will try to find the optimum point
##' for each \code{n} value using optimization algorithms. When
##' \code{h}, \code{L} and \code{n} are all given, \code{ecoXbar}
##' function will use the "grid method" to calculate the optimum point,
##' that is ECH for all the combinations of the parameters will be
##' calculated. The "grid method" way is much slower than using
##' optimization algorithms, but it would be a good choice when
##' optimization algorithms fail to work well.
##'
##' For cost parameters either {P0, P1} or {C0, C1} is needed.  If P0
##' and P1 are given, they will be used first, else C0 and C1 will be
##' used.  For economic design of the X-bar chart, when \code{d1} and
##' \code{d2} are both 1, only if the difference between P0 and P1
##' keeps the same, the results are identical. If the difference
##' between C0 and C1 keeps the same, the optimum parameters are
##' almost the same but the ECH(Expected Cost per Hour) values will
##' change.
##'
##' \code{echXbar} is used to calculate the ECH (Expected Cost per
##' Hour) for one given design point.
##'
##' @title Economic design for the X-bar control chart
##' @param h sampling interval. It can be a numeric vector or left
##' undefined. See 'Details'
##' @param L number of standard deviations from control limits to
##' center line. It can be a numeric vector or left undefined. See
##' 'Details'
##' @param n sample size. It can be an integer vector or left
##' undefined. See 'Details'
##' @param lambda we assume the in-control time follows an exponential
##' distribution with mean 1/lambda. Default value is 0.05.
##' @param delta shift in process mean in standard deviation units
##' when assignable cause occurs (delta = (mu1 - mu0)/sigma), where
##' sigma is the standard deviation of observations; mu0 is the
##' in-control process mean; mu1 is the out-of-control process mean.
##' Default value is 2.
##' @param P0 profit per hour earned by the process operating in
##' control. See 'Details'.
##' @param P1 profit per hour earned by the process operating out of
##' control(P0 > P1).
##' @param C0 cost per hour due to nonconformities produced while the
##' process is in control.
##' @param C1 cost per hour due to nonconformities produced while the
##' process is out of control.(C1 > C0)
##' @param Cr cost for searching and repairing the assignable cause,
##' including any downtime.
##' @param Cf cost per false alarm, including the cost of searching
##' for the cause and the cost of downtime if production ceases during
##' search.
##' @param T0 time to sample and chart one item.
##' @param Tc expected time to discover the assignable cause.
##' @param Tf expected search time when false alarm occurs.
##' @param Tr expected time to repair the process.
##' @param a fixed cost per sample.
##' @param b cost per unit sampled.
##' @param d1 flag for whether production continues during searches
##' (1-yes, 0-no). Default value is 1.
##' @param d2 flag for whether production continues during repairs
##' (1-yes, 0-no). Default value is 1.
##' @param nlevels number of contour levels desired. Default value is
##' 30. It works only when \code{contour.plot} is TRUE.
##' @param sided distinguish between one- and two-sided X-bar chart by
##' choosing one'' or two'' respectively.  When \code{sided =
##' one''}, \code{delta > 0} means the control chart for detecting a
##' positive shift, and vice versa. Default is two''.
##' @param par initial values for the parameters to be optimized
##' over. It can be a vector of length 2 or 3. See 'Details'
##' @param contour.plot a logical value indicating whether a contour
##' plot should be drawn. Default is FALSE. Only works when the
##' parameters h, L and n are all specified.
##' @param call.print a logical value indicating whether the "call"
##' should be printed on the contour plot. Default is TRUE
##' @param ... other arguments to be passed to \code{optim} function.
##' @return The \code{ecoXbar} function returns an object of class
##' "edcc", which is a list of elements \code{optimum},
##' \code{cost.frame}, \code{FAR} and \code{ATS}. \code{optimum} is a
##' vector with the optimum parameters and the corresponding ECH
##' value; \code{cost.frame} is a dataframe with the optimum
##' parameters and the corresponding ECH values for all given
##' \code{n}(if \code{n} is not specified, \code{cost.frame} won't be
##' returned); \code{FAR} indicates the false alarm rate during the
##' in-control time, which is calculated as lambda*(average number of
##' false alarm); \code{ATS} indicates the average time to signal
##' after the occurrence of an assignable cause, calculated as h*ARL2
##' - tau, where tau is the expected time of occurrence of the
##' assignable cause given it occurs between the i-th and (i+1)st
##' samples.
##' The \code{echXbar} function returns the calculated ECH value only.
##' @references Weicheng Zhu, Changsoon Park (2013), {edcc: An R
##' Package for the Economic Design of the Control Chart}. \emph{Journal
##' of Statistical Software}, 52(9),
##' 1-24. \url{http://www.jstatsoft.org/v52/i09/}
##'
##' Douglas (2009). \emph{Statistical quality control: a modern
##' introduction}, sixth edition,  463-471.
##'
##' Lorenzen and Vance (1986). The economic design of control charts:
##' a unified approach, \emph{Technometrics}, 28. 3-10.
##' @examples
##' # Douglas (2009). Statistical quality control: a modern introduction, sixth edition, p470.
##' # In control profit per hour is 110, out of control profit per hour is 10
##' ecoXbar(P0=110,P1=10)
##' # In control profit per hour is 150, out of control profit per hour
##' # is 50, the result is identical with the previous one, because the
##' #difference between P0 and P1 are the same
##' ecoXbar(P0=150,P1=50)
##' # In control cost per hour is 0, out of control cost per hour is 100.
##' # The result is the same with the previous one
##' ecoXbar(C0=0,C1=100)
##' # The optimum parameters are the same with the previous one,
##' # but Cost values are different. See 'details'
##' ecoXbar(C0=10,C1=110)
##' # Based on the global optimum above, we specify the range of the
##' # parameters like this
##' x <- ecoXbar(h=seq(0.7,0.9,by=.01),L=seq(2.8,3.3,by=.01),n=4:6,P0=110,
##' P1=10,nlevels=50,contour.plot=TRUE)
##' x
##' # Modify the contour plot
##' contour(x,nlevels=20,lty=2,col=2,call.print=FALSE)
##' # update the parameters
##' update(x,P0=NULL,P1=NULL,C0=10,C1=110)
##' @export
ecoXbar <- function(h, L, n, delta = 2, lambda = .05, P0 = NULL, P1 = NULL, C0 = NULL, C1 = NULL, Cr = 25, Cf = 50, T0 = 0.0167, Tc = 1, Tf = 0, Tr = 0, a = 1, b = .1, d1 = 1, d2 = 1, nlevels = 30, sided = "two", par = NULL, contour.plot = FALSE, call.print = TRUE, ...){
if(!is.null(par)){
if(!is.vector(par)) stop("par should be a vector of length 2 or 3!") else
if(!(length(par)==3 || length(par)==2)) stop("par should be a vector of length 2 or 3!") else
if(length(par)==3){
opt <- optim(par,.echXbar1,lambda=lambda,delta=delta,P0=P0,P1=P1,C0=C0,C1=C1,Cr=Cr,Cf=Cf,T0=T0,Tc=Tc,Tf=Tf,Tr=Tr,a=a,b=b,d1=d1,d2=d2,sided=sided,...)
optimum <- as.data.frame(t(c(opt$par,opt$value)))
gn <- round(as.numeric(optimum[3])) # global optimum of n
cost.frame <- NULL
for(k in c(gn-1, gn, gn+1)){
opt <- optim(par[1:2], .echXbar2,n=k,lambda=lambda,delta=delta,P0=P0,P1=P1,C0=C0,C1=C1,Cr=Cr,Cf=Cf,T0=T0,Tc=Tc,Tf=Tf,Tr=Tr,a=a,b=b,d1=d1,d2=d2,sided=sided,...)
cost.frame <- rbind(cost.frame,c(opt$par,k,opt$value))
}
optimum <- cost.frame[which(cost.frame[,4]==min(cost.frame[,4])),]
names(optimum) <- c("Optimum h","Optimum L","Optimum n","ECH")
optXbar <- list(optimum=optimum)
} else{
cost.frame <- NULL
for(k in n){
opt <- optim(par,.echXbar2,n=k,lambda=lambda,delta=delta,P0=P0,P1=P1,C0=C0,C1=C1,Cr=Cr,Cf=Cf,T0=T0,Tc=Tc,Tf=Tf,Tr=Tr,a=a,b=b,d1=d1,d2=d2,sided=sided,...)
cost.frame <- rbind(cost.frame,c(opt$par,k,opt$value))
}
optimum <- cost.frame[which(cost.frame[,4]==min(cost.frame[,4])),]
names(optimum) <- c("Optimum h","Optimum L","Optimum n","ECH")
colnames(cost.frame) <- c("Optimum h","Optimum L","Optimum n","ECH")
rownames(cost.frame) <- rep("",length(n))
optXbar <- list(optimum=optimum,cost.frame=cost.frame)
}
}else
if(missing(h) && missing(L) && missing(n)){
opt <- optim(c(1,3,5),.echXbar1,lambda=lambda,delta=delta,P0=P0,P1=P1,C0=C0,C1=C1,Cr=Cr,Cf=Cf,T0=T0,Tc=Tc,Tf=Tf,Tr=Tr,a=a,b=b,d1=d1,d2=d2,sided=sided,...)
#    cat('The optimum parameters maybe around h=', opt$par[1],' L=', opt$par[2],' n=', opt$par[3],' and the minimum ECH is about', opt$value,'\n')
optimum <- as.data.frame(t(c(opt$par,opt$value)))
gn <- round(as.numeric(optimum[3])) # global optimum of n
cost.frame <- NULL
for(k in c(gn-1, gn, gn+1)){
opt <- optim(c(1,3),.echXbar2,n=k,delta=delta,lambda=lambda,P0=P0,P1=P1,C0=C0,C1=C1,Cr=Cr,Cf=Cf,T0=T0,Tc=Tc,Tf=Tf,Tr=Tr,a=a,b=b,d1=d1,d2=d2,sided=sided,...)
cost.frame <- rbind(cost.frame,c(opt$par,n=k,opt$value))
}
optimum <- cost.frame[which(cost.frame[,4]==min(cost.frame[,4])),]
names(optimum) <- c("Optimum h","Optimum L","Optimum n","ECH")
optXbar <- list(optimum=optimum)
}else
if(missing(h) && missing(L)){
cost.frame <- NULL
for(k in n){
opt <- optim(c(1,3),.echXbar2,n=k,delta=delta,lambda=lambda,P0=P0,P1=P1,C0=C0,C1=C1,Cr=Cr,Cf=Cf,T0=T0,Tc=Tc,Tf=Tf,Tr=Tr,a=a,b=b,d1=d1,d2=d2,sided=sided,...)
cost.frame <- rbind(cost.frame,c(opt$par,n=k,opt$value))
}
optimum <- cost.frame[which(cost.frame[,4]==min(cost.frame[,4])),]
names(optimum) <- c("Optimum h","Optimum L","Optimum n","ECH")
colnames(cost.frame) <- c("Optimum h","Optimum L","Optimum n","ECH")
rownames(cost.frame) <- rep("",length(n))
optXbar <- list(optimum=optimum,cost.frame=cost.frame)
}else{
cost.frame <- NULL
opt.mat=outer(h,L,FUN=echXbar,n=n[1],delta=delta,lambda=lambda,P0=P0,P1=P1,C0=C0,C1=C1,Cr=Cr,Cf=Cf,T0=T0,Tc=Tc,Tf=Tf,Tr=Tr,a=a,b=b,d1=d1,d2=d2,sided=sided)
aa <- which(opt.mat==min(opt.mat),arr.ind=TRUE)
opt.ech <- min(opt.mat)
cost.frame <- rbind(cost.frame,c(h[aa[1,1]],L[aa[1,2]],n[1],opt.ech))
for(k in n[-1]){
mat=outer(h,L,FUN=echXbar,n=k,delta=delta,lambda=lambda,P0=P0,P1=P1,C0=C0,C1=C1,Cr=Cr,Cf=Cf,T0=T0,Tc=Tc,Tf=Tf,Tr=Tr,a=a,b=b,d1=d1,d2=d2,sided=sided)
aa <- which(mat==min(mat),arr.ind=TRUE)
cost.frame <- rbind(cost.frame,c(h[aa[1,1]],L[aa[1,2]],k, ech <- min(mat)))
if(ech < opt.ech){
opt.ech <- ech; opt.mat <- mat
}
}
optimum <- cost.frame[which(cost.frame[,4]==min(cost.frame[,4])),]
names(optimum) <- c("Optimum h","Optimum L","Optimum n","ECH")
colnames(cost.frame) <- c("Optimum h","Optimum L","Optimum n","ECH")
rownames(cost.frame) <- rep("",length(n))
## FIXME
if(contour.plot){
if(call.print){
ca <- match.call()
ca[[1]] <- NULL
par(mar=c(6.1,4.1,4.1,2.1))
contour(h,L,opt.mat,main=strwrap(paste("ecoXbar(",paste(names(ca),"=",unlist(ca),sep="",collapse=", "),")")),xlab="h",ylab="L",nlevels=nlevels)
} else{
par(mar=c(6.1,4.1,2.1,2.1))
contour(h,L,opt.mat,xlab="h",ylab="L",nlevels=nlevels)
}
mtext(sprintf('Opt h=%s   Opt L=%s   n=%s   ECH=%s',optimum[1], optimum[2], optimum[3],round(optimum[4],digits=4)),side=1,line=4.5)
points(optimum[1],optimum[2],pch=3)
}
optXbar <- list(optimum=optimum,cost.frame=cost.frame, opt.mat=opt.mat)
}
## FIXME: alpha~sided
opth <- as.numeric(optXbar$optimum[1]) optL <- as.numeric(optXbar$optimum[2])
optn <- as.numeric(optXbar$optimum[3]) if(sided == "one"){ ARL1 <- 1/pnorm(-optL) ARL2 <- 1/pnorm(-optL + abs(delta*sqrt(optn))) } if(sided == "two"){ alpha <- 2*pnorm(-optL) beta <- pnorm(optL - delta*sqrt(optn))-pnorm(-optL- delta*sqrt(optn)) ARL1 <- 1/alpha ARL2 <- 1/(1-beta) } s <- 1/(exp(lambda*opth)-1) tau <- (1-(1+lambda*opth)*exp(-lambda*opth))/(lambda*(1-exp(-lambda*opth))) optXbar$FAR <- lambda*s/ARL1
optXbar$ATS <- ARL2*opth - tau optXbar$call <- match.call()
return(structure(optXbar,class="edcc"))
}
##' @rdname ecoXbar
##' @export
# calculate the ECH for given parameters
echXbar <- function(h,L,n,delta=2,lambda=.05,P0=NULL,P1=NULL,C0=NULL,C1=NULL,Cr=25,Cf=50,T0=0.0167,Tc=1,Tf=0,Tr=0,a=1,b=.1,d1=1,d2=1, sided = "two"){
delta.std <- delta*sqrt(n)
if(sided == "one"){
ARL1 <- 1/pnorm(-L)
ARL2 <- 1/pnorm(-L + abs(delta.std))
}
if(sided == "two"){
alpha <- 2*pnorm(-L)
beta <- pnorm(L - delta.std)-pnorm(- L- delta.std)
ARL1 <- 1/alpha
ARL2 <- 1/(1-beta)
}
tau <- (1-(1+lambda*h)*exp(-lambda*h))/(lambda*(1-exp(-lambda*h)))
s <- 1/(exp(lambda*h)-1)
if(!is.null(P0)&!is.null(P1)){
ECT <- 1/lambda+(1-d1)*s*Tf/ARL1 - tau + n*T0 + h*ARL2 + Tc + Tr
ECP <- P0/lambda + P1*(-tau+n*T0+h*ARL2+d1*Tc+d2*Tr) - s*Cf/ARL1 - Cr - (a+b*n)*(1/lambda-tau+n*T0+h*ARL2+d1*Tc+d2*Tr)/h
ECH <- P0 - ECP/ECT
}else
if(!is.null(C0)&!is.null(C1)){
ECT <- 1/lambda+(1-d1)*s*Tf/ARL1 - tau + n*T0 + h*ARL2 + Tc + Tr
ECC <- C0/lambda + C1*(-tau+n*T0+h*ARL2+d1*Tc+d2*Tr) + s*Cf/ARL1+Cr+(a+b*n)*(1/lambda-tau+n*T0+h*ARL2+d1*Tc+d2*Tr)/h
ECH <- ECC/ECT
}else
stop("You should at least give a pair of value to P0,P1 or C0,C1")
return(ECH)
}

##' Calculate the optimum parameters of n(sample size), h(sampling
##' interval), k(reference value) and H(decision interval) for
##' the reference value see 'Details'.
##'
##' When parameter \code{par} is specified, optimization algorithms
##' would be used as default. \code{par} can be specified as:
##' \code{par = c(h, H)} where \code{h} and \code{H} are the initial
##' values of smapling interval and decision interval when \code{n} is
##' specified; or \code{par = c(h, H, n)}. Good inital values may lead
##' to good optimum results.
##'
##' When parameters \code{h}, \code{H}, \code{n} are all undefined,
##' \code{ecoCusum} function will try to find the global optimum point
##' to minimize the ECH (Expected Cost per Hour) using optimization
##' algorithms (\code{optim} function), but in this case \code{n}
##' would not be integer. It is usually helpful for the experimenter
##' to find the region where the optimum point may exist quickly. When
##' \code{h} and \code{H} are undefined but \code{n} is given as an
##' integer vector, \code{ecoCusum} function will try to find the
##' optimum point for each \code{n} value using optimization
##' algorithms. When \code{h}, \code{H} and \code{n} are all given,
##' \code{ecoCusum} function will use a "grid method" way to calculate the
##' optimum point, that is ECH for all the combinations of the
##' parameters will be calculated. The "grid method" way is much slower
##' than using optimization algorithms, but it would be a good choice
##' when optimization algorithms fail to work well.
##'
##' There is strong numerical and theoretical evidence that for given
##' L1, the value of L0 approaches its maximum when k(reference value)
##' is chosen mid-way the between AQL and the RQL: $k = mu0 + ##' 0.5*delta*sigma (Appl. Statist.(1974) 23, No. 3, p. 420). For this ##' reason we treat k as a constant value and optimize n, h and H. ##' For cost parameters either {P0, P1} or {C0, C1} is needed. If P0 ##' and P1 are given, they will be used first, else C0 and C1 will be ##' used. For economic design of the CUSUM chart, when \code{d1} and ##' \code{d2} are both 1, only if the difference between P0 and P1 ##' keeps the same, the results are identical. If the difference ##' between C0 and C1 keeps the same, the optimum parameters are ##' almost the same but the ECH(Expected Cost per Hour) values will ##' change. ##' ##' \code{echCusum} is used to calculate the ECH (Expected Cost per ##' Hour) for one given design point. ##' ##' @title Economic design for the CUSUM control chart ##' @param h sampling interval. It can be a numeric vector or left ##' undefined. See 'Details' ##' @param H decision interval. It can be a numeric vector or left ##' undefined. See 'Details' ##' @param n sample size. It can be an integer vector or left ##' undefined. See 'Details' ##' @param delta shift in process mean in standard deviation units ##' when assignable cause occurs (delta = |mu1 - mu0|/sigma), where ##' sigma is the standard deviation of observations; mu0 is the ##' in-control process mean; mu1 is the out-of-control process mean. ##' Default value is 2. ##' @param lambda we assume the in-control time follows a exponential ##' distribution with mean 1/lambda. Default value is 0.05. ##' @param P0 profit per hour earned by the process operating in ##' control. See 'Details'. ##' @param P1 profit per hour earned by the process operating out of ##' control ##' @param C0 cost per hour due to nonconformities produced while the ##' process is in control. ##' @param C1 cost per hour due to nonconformities produced while the ##' process is out of control.(C1 > C0) ##' @param Cr cost for searching and repairing the assignable cause, ##' including any downtime. ##' @param Cf cost per false alarm, including the cost of searching ##' for the cause and the cost of downtime if production ceases during ##' search. ##' @param T0 time to sample and chart one item. ##' @param Tc expected time to discover the assignable cause. ##' @param Tf expected search time when false alarm occurs. ##' @param Tr expected time to repair the process. ##' @param a fixed cost per sample. ##' @param b cost per unit sampled. ##' @param d1 flag for whether production continues during searches ##' (1-yes, 0-no). Default value is 1. ##' @param d2 flag for whether production continues during repairs ##' (1-yes, 0-no). Default value is 1. ##' @param nlevels ##' 30. It works only when \code{contour.plot} is TRUE. ##' @param sided distinguish between one-, two-sided and Crosier's ##' modified two-sided CUSUM scheme by choosing "one", "two", and ##' "Crosier", respectively. See details in \code{\link[spc]{xcusum.arl}} ##' @param par initial values for the parameters to be optimized ##' over. It can be a vector of length 2 or 3. See 'Details' ##' @param contour.plot a logical value indicating wether a contour ##' plot should be drawn. Default is FALSE. ##' @param call.print a logical value indicating whether the "call" ##' should be printed on the contour plot. Default is TRUE ##' @param ... other arguments to be passed to \code{optim} function. ##' @return The \code{ecoCusum} function returns an object of class ##' "edcc", which is a list of elements \code{optimum}, ##' \code{cost.frame}, \code{FAR} and \code{ATS}. \code{optimum} is a ##' vector with the optimum parameters and the corresponding ECH ##' value; \code{cost.frame} is a dataframe with the optimum ##' parameters and the corresponding ECH values for all given ##' \code{n}(if \code{n} is not specified, \code{cost.frame} won't be ##' returned); \code{FAR} indicates the false alarm rate during the ##' in-control time, which is calculated as lambda*(average number of ##' false alarm); \code{ATS} indicates the average time to signal ##' after the occurrence of an assignable cause, calculated as h*ARL2 ##' - tau, where tau is the expected time of occurrence of the ##' assignable cause given it occurs between the i-th and (i+1)st ##' samples. ##' The \code{echCusum} function returns the calculated ECH value only. ##' @seealso \code{\link{ecoXbar}}, \code{\link{ecoEwma}}, ##' \code{\link[spc]{xcusum.arl}}, \code{\link[stats]{optim}}, ##' \code{\link{update.edcc}}, \code{\link{contour}} ##' @references Weicheng Zhu, Changsoon Park (2013), {edcc: An R ##' Package for the Economic Design of the Control ##' Chart}. \emph{Journal of Statistical Software}, 52(9), ##' 1-24. \url{http://www.jstatsoft.org/v52/i09/} ##' ##' Lorenzen and Vance (1986). The economic design of ##' control charts: a unified approach, \emph{Technometrics}, 28. 3-10. ##' ##' Chiu, W.K. (1974). The economic design of CUSUM charts for ##' controlling normal means, \emph{Journal of the Royal Statistical ##' Society. Series C (Applied Statistics)}, 23(3), 420-433. ##' @examples ##' #Chiu, W.K. (1974). Applied Statistics, 23, p427 Table3, row 1-2,14 ##' ## LINE 1 ##' ## global optimization to h, H and n, when lambda = 0.01, "Nelder-Mead" optimization algorithm doesn't work ##' #(y <- ecoCusum( P0=150,P1=50,Cr=30,d1=0,d2=0)) ##' ## we can try other algorithms: ##' (y1 <- ecoCusum( P0=150,P1=50,Cr=30,d1=0,d2=0,method="BFGS")) ##' # Based on the global optimum above, we specify the range of the ##' # parameters like this ##' (yy1 <- ecoCusum( h=seq(1.3,1.45,by=.01), H=seq(.5,0.6,by=.01),n=4:6, ##' P0=150,P1=50,Cr=30,d1=0,d2=0)) ##' ## LINE 2 ##' (y2 <- ecoCusum( P0=150,P1=50,Cr=30,d1=0,d2=0,lambda=0.05)) ##' (yy2 <- ecoCusum( h=seq(.6,0.7,by=.01), H=seq(.5,0.6,by=.01),n=3:6, ##' P0=150,P1=50,Cr=30,d1=0,d2=0,lambda=0.05)) ##' contour(yy2) ##' ## LINE 14 ##' (y14 <- ecoCusum(n=30,P0=150,P1=50,Cr=30,delta=0.5,d1=0,d2=0,method="L-BFGS-B")) ##' (yy14 <- ecoCusum(h=seq(2.55,2.65,by=0.01),H=seq(0.3,0.4,by=0.01), ##' n=28:30,P0=150,P1=50,Cr=30,delta=0.5,d1=0,d2=0)) ##' #Douglas (2009). Statistical quality control: a modern introduction, sixth edition, p470. ##' ecoCusum(lambda=.05,P0=110,P1=10,Cr=25,Cf=50,Tr=0,Tf=0,Tc=1,T0=.0167,a=1) ##' ecoCusum(h=seq(0.75,0.85,by=.01),H=seq(.55,0.65,by=.01),n=4:6,lambda=.05, ##' P0=110,P1=10,Cr=25,Cf=50,Tr=0,Tf=0,Tc=1,T0=.0167,a=1) ##' @export ecoCusum <- function( h, H,n,delta = 2,lambda = .01, P0 = NULL, P1 = NULL,C0 = NULL,C1 = NULL, Cr = 20, Cf = 10,T0 = 0, Tc = .1,Tf = .1,Tr = 0.2, a = .5, b = .1, d1 = 1, d2 = 1, nlevels=30, sided = "one", par=NULL, contour.plot=FALSE, call.print=TRUE,...){ if(!is.null(par)){ if(!is.vector(par)) stop("par should be a vector of length 2 or 3!") else if(!(length(par)==3 || length(par)==2)) stop("par should be a vector of length 2 or 3!") else if(length(par)==3){ opt <- optim(par,.echCusum1,delta = delta,lambda = lambda, P0 = P0, P1 = P1, C0 = C0, C1 = C1, Cr = Cr, Cf = Cf,T0 = T0, Tc = Tc,Tf = Tf,Tr = Tr, a = a, b = b, d1 = d1, d2 = d2,sided = sided,...) # cat('The optimum parameters maybe around h=', opt$par[1],' H=', opt$par[2],' n=', opt$par[3],' and the minimum ECH is about', opt$value,'\n') optimum <- as.data.frame(t(c(opt$par,opt$value))) gn <- round(as.numeric(optimum[3])) # global optimum of n cost.frame <- NULL for(k in c(gn-1, gn, gn+1)){ opt <- optim(par[1:2],.echCusum2,n=k,delta = delta,lambda = lambda, P0 = P0, P1 = P1, C0 = C0, C1 = C1, Cr = Cr, Cf = Cf,T0 = T0, Tc = Tc,Tf = Tf,Tr = Tr, a = a, b = b, d1 = d1, d2 = d2,sided = sided,...) cost.frame <- rbind(cost.frame,c(opt$par,k,opt$value)) } optimum <- cost.frame[which(cost.frame[,4]==min(cost.frame[,4])),] names(optimum) <- c("Optimum h","Optimum H","Optimum n","ECH") optCusum <- list(optimum=optimum) } else{ cost.frame <- NULL for(k in n){ opt <- optim(par,.echCusum2,n=k,delta = delta,lambda = lambda, P0 = P0, P1 = P1, C0 = C0, C1 = C1, Cr = Cr, Cf = Cf,T0 = T0, Tc = Tc,Tf = Tf,Tr = Tr, a = a, b = b, d1 = d1, d2 = d2,sided = sided,...) cost.frame <- rbind(cost.frame,c(opt$par,k,opt$value)) } optimum <- cost.frame[which(cost.frame[,4]==min(cost.frame[,4])),] names(optimum) <- c("Optimum h","Optimum H","Optimum n","ECH") colnames(cost.frame) <- c("Optimum h","Optimum H","Optimum n","ECH") rownames(cost.frame) <- rep("",length(n)) optCusum <- list(optimum=optimum,cost.frame=cost.frame) } }else if(missing(h) && missing(H) && missing(n)){ opt <- optim(c(1,1,5),.echCusum1,delta = delta,lambda = lambda, P0 = P0, P1 = P1, C0 = C0, C1 = C1, Cr = Cr, Cf = Cf,T0 = T0, Tc = Tc,Tf = Tf,Tr = Tr, a = a, b = b, d1 = d1, d2 = d2,sided = sided,...) # cat('The optimum parameters maybe around h=', opt$par[1],' H=', opt$par[2],' n=', opt$par[3],' and the minimum ECH is about', opt$value,'\n') optimum <- as.data.frame(t(c(opt$par,opt$value))) gn <- round(as.numeric(optimum[3])) # global optimum of n cost.frame <- NULL for(k in c(gn-1, gn, gn+1)){ opt <- optim(c(1,1),.echCusum2,n=k,delta = delta,lambda = lambda, P0 = P0, P1 = P1, C0 = C0, C1 = C1, Cr = Cr, Cf = Cf,T0 = T0, Tc = Tc,Tf = Tf,Tr = Tr, a = a, b = b, d1 = d1, d2 = d2,sided = sided,...) cost.frame <- rbind(cost.frame,c(opt$par,k,opt$value)) } optimum <- cost.frame[which(cost.frame[,4]==min(cost.frame[,4])),] names(optimum) <- c("Optimum h","Optimum H","Optimum n","ECH") optCusum <- list(optimum=optimum) }else if(missing(h) && missing(H)){ cost.frame <- NULL for(k in n){ opt <- optim(c(1,1),.echCusum2,n=k,delta = delta,lambda = lambda, P0 = P0, P1 = P1, C0 = C0, C1 = C1, Cr = Cr, Cf = Cf,T0 = T0, Tc = Tc,Tf = Tf,Tr = Tr, a = a, b = b, d1 = d1, d2 = d2,sided = sided,...) cost.frame <- rbind(cost.frame,c(opt$par,n=k,opt$value)) } optimum <- cost.frame[which(cost.frame[,4]==min(cost.frame[,4])),] names(optimum) <- c("Optimum h","Optimum H","Optimum n","ECH") colnames(cost.frame) <- c("Optimum h","Optimum H","Optimum n","ECH") rownames(cost.frame) <- rep("",length(n)) optCusum <- list(optimum=optimum,cost.frame=cost.frame) }else{ cost.frame <- NULL opt.mat <- matrix(NA,length(h),length(H)) for(i in 1:length(h)) for(j in 1:length(H)) opt.mat[i,j] <- echCusum(h[i],H[j],n[1],delta = delta,lambda = lambda, P0 = P0, P1 = P1, C0 = C0, C1 = C1, Cr = Cr, Cf = Cf,T0 = T0, Tc = Tc,Tf = Tf,Tr = Tr, a = a, b = b, d1 = d1, d2 = d2,sided = sided) aa <- which(opt.mat==min(opt.mat),arr.ind=TRUE) opt.ech <- min(opt.mat) cost.frame <- rbind(cost.frame,c(h[aa[1,][1]],H[aa[1,][2]],n[1],opt.ech)) for(nk in n[-1]){ mat <- matrix(NA,length(h),length(H)) for(i in 1:length(h)) for(j in 1:length(H)) mat[i,j] <- echCusum(h[i],H[j],nk,delta = delta,lambda = lambda, P0 = P0, P1 = P1, C0 = C0, C1 = C1, Cr = Cr, Cf = Cf,T0 = T0, Tc = Tc,Tf = Tf,Tr = Tr, a = a, b = b, d1 = d1, d2 = d2,sided = sided) aa <- which(mat==min(mat),arr.ind=TRUE) cost.frame <- rbind(cost.frame,c(h[aa[1,][1]],H[aa[1,][2]],nk, ech <- min(mat))) if(ech < opt.ech){ opt.ech <- ech; opt.mat <- mat } } optimum <- cost.frame[which(cost.frame[,4]==min(cost.frame[,4])),] names(optimum) <- c("Optimum h","Optimum H","Optimum n","ECH") colnames(cost.frame) <- c("Optimum h","Optimum H","Optimum n","ECH") rownames(cost.frame) <- rep("",length(n)) if(contour.plot){ if(call.print){ ca <- match.call() ca[[1]] <- NULL par(mar=c(6.1,4.1,4.1,2.1)) contour(h,H,opt.mat,main=strwrap(paste("ecoCusum(",paste(names(ca),"=",unlist(ca),sep="",collapse=", "),")")),xlab="h",ylab="H",nlevels=nlevels) } else{ par(mar=c(6.1,4.1,2.1,2.1)) contour(h,H,opt.mat,xlab="h",ylab="H",nlevels=nlevels) } mtext(sprintf('Opt h=%s Opt H=%s n=%s ECH=%s',optimum[1], optimum[2], optimum[3],round(optimum[4],digits=4)),side=1,line=4.5) points(optimum[1],optimum[2],pch=3) } optCusum <- list(optimum=optimum, cost.frame=cost.frame, opt.mat=opt.mat) } opth <- as.numeric(optCusum$optimum[1])
optH <- as.numeric(optCusum$optimum[2]) optn <- as.numeric(optCusum$optimum[3])
delta.std <- sqrt(as.numeric(optn))*delta
ARL1 <- as.numeric(xcusum.arl(0.5*delta.std, optH, 0, sided=sided))
ARL2 <- as.numeric(xcusum.arl(0.5*delta.std, optH, delta.std, sided=sided))
s <- 1/(exp(lambda*opth)-1)
tau <- (1-(1+lambda*opth)*exp(-lambda*opth))/(lambda*(1-exp(-lambda*opth)))
optCusum$FAR <- lambda*s/ARL1 optCusum$ATS <- ARL2*opth - tau
optCusum$call <- match.call() return(structure(optCusum,class="edcc")) } ##' @rdname ecoCusum ##' @export # calculate the ECH for given parameters echCusum <- function(h,H,n,delta = 2,lambda = .01, P0 = NULL, P1 = NULL,C0 = NULL,C1 = NULL, Cr = 20, Cf = 10,T0 = 0, Tc = .1,Tf = .1,Tr = 0.2, a = .5, b = .1, d1 = 1, d2 = 1,sided = "one"){ delta.std <- sqrt(n)*delta #standardization for delta k <- delta.std/2 ARL1 <- as.numeric(xcusum.arl(k,H,0,sided=sided)) ARL2 <- as.numeric(xcusum.arl(k,H,delta.std,sided=sided)) tau <- (1-(1+lambda*h)*exp(-lambda*h))/(lambda*(1-exp(-lambda*h))) s <- 1/(exp(lambda*h)-1) if(!is.null(P0)&!is.null(P1)){ ECT <- 1/lambda+(1-d1)*s*Tf/ARL1 - tau + n*T0 + h*ARL2 + Tc + Tr ECP <- P0/lambda + P1*(-tau+n*T0+h*ARL2+d1*Tc+d2*Tr) - s*Cf/ARL1 - Cr - (a+b*n)*(1/lambda-tau+n*T0+h*ARL2+d1*Tc+d2*Tr)/h ECH <- P0 - ECP/ECT }else if(!is.null(C0)&!is.null(C1)){ ECT <- 1/lambda+(1-d1)*s*Tf/ARL1 - tau + n*T0 + h*ARL2 + Tc + Tr ECC <- C0/lambda + C1*(-tau+n*T0+h*ARL2+d1*Tc+d2*Tr) + s*Cf/ARL1+Cr+(a+b*n)*(1/lambda-tau+n*T0+h*ARL2+d1*Tc+d2*Tr)/h ECH <- ECC/ECT }else stop("You should at least give a pair of value to P0,P1 or C0,C1") return(ECH) } ##' Calculate the optimum parameters, n(sample size), h(sampling ##' interval), w(weight to the present sample) and k(number of ##' s.d. from control limits to center line) for economic Design of ##' the EWMA control chart . ##' ##' Parameter \code{w} should always be given, because the range of ##' \code{w} is so restricted that optimization algorithms usually ##' don't converge. ##' ##' When parameter \code{par} is specified, optimization algorithms ##' would be used as default. \code{par} can be specified as: ##' \code{par = c(h, k)} where \code{h} and \code{k} are the initial ##' values of smapling interval and control limit when \code{n} is ##' specified; or \code{par = c(h, k, n)}. Good inital values may lead ##' to good optimum results. ##' ##' When parameters \code{h}, \code{k}, \code{n} are ##' all undefined, \code{ecoEwma} function will try to find the global ##' optimum point to minimize the ECH (Expected Cost per Hour) using ##' optimization algorithms (\code{optim} function), but in this case ##' \code{n} would not be integer. It is usually helpful for the ##' experimenter to find the region where the optimum point may exist ##' quickly. When \code{h} and \code{k} are undefined but \code{n} is ##' given as an integer vector, \code{ecoEwma} function will try to ##' find the optimum point for each \code{n} value using optimization ##' algorithms. When \code{h}, \code{k} and \code{n} are all given, ##' \code{ecoEwma} function will use a "grid method" way to calculate the ##' optimum point, that is ECH for all the combinations of the ##' parameters will be calculated. The "grid method" way is much slower ##' than using optimization algorithms, but it would be a good choice ##' when optimization algorithms fail to work well. ##' ##' For cost parameters either {P0, P1} or {C0, C1} is needed. If P0 ##' and P1 are given, they will be used first, else C0 and C1 will be ##' used. For economic design of the EWMA chart, when \code{d1} and ##' \code{d2} are both 1, only if the difference between P0 and P1 ##' keeps the same, the results are identical. If the difference ##' between C0 and C1 keeps the same, the optimum parameters are ##' almost the same but the ECH(Expected Cost per Hour) values will ##' change. ##' ##' \code{echEwma} is used to calculate the ECH (Expected Cost per ##' Hour) for one given design point. ##' ##' @title Economic design for the EWMA control chart ##' @param h sampling interval. It can be a numeric vector or left ##' undefined. See 'Details' ##' @param w the weight value between 0 and 1 given to the latest ##' sample. It must be specified. ##' @param k control limit coefficient. It can be a numeric vector or ##' left undefined. See 'Details' ##' @param n sample size. It can be an integer vector or left ##' undefined. See 'Details' ##' @param delta shift in process mean in standard deviation units ##' when assignable cause occurs (delta = |mu1 - mu0|/sigma), where ##' sigma is the standard deviation of observations; mu0 is the ##' in-control process mean; mu1 is the out-of-control process mean. ##' Default value is 2. ##' @param lambda we assume the in-control time follows a exponential ##' distribution with mean 1/lambda. Default value is 0.05. ##' @param P0 profit per hour earned by the process operating in ##' control. See 'Details'. ##' @param P1 profit per hour earned by the process operating out of ##' control ##' @param C0 cost per hour due to nonconformities produced while the ##' process is in control. ##' @param C1 cost per hour due to nonconformities produced while the ##' process is out of control.(C1 > C0) ##' @param Cr cost for searching and repairing the assignable cause, ##' including any downtime. ##' @param Cf cost per false alarm, including the cost of searching ##' for the cause and the cost of downtime if production ceases during ##' search. ##' @param T0 time to sample and chart one item. ##' @param Tc expected time to discover the assignable cause. ##' @param Tf expected search time when false alarm occurs. ##' @param Tr expected time to repair the process. ##' @param a fixed cost per sample. ##' @param b cost per unit sampled. ##' @param d1 flag for whether production continues during searches ##' (1-yes, 0-no). Default value is 1. ##' @param d2 flag for whether production continues during repairs ##' (1-yes, 0-no). Default value is 1. ##' @param nlevels number of contour levels desired. Default value is ##' 30. It works only when \code{contour.plot} is TRUE. ##' @param sided distinguish between one- and two-sided EWMA control ##' chart by choosing "one" and "two", respectively. See details in ##' \code{\link[spc]{xewma.arl}} ##' @param par initial values for the parameters to be optimized ##' over. It can be a vector of length 2 or 3. See 'Details' ##' @param contour.plot a logical value indicating whether a contour ##' plot should be drawn. Default is FALSE. ##' @param call.print a logical value indicating whether the "call" ##' should be printed on the contour plot. Default is TRUE ##' @param ... other arguments to be passed to contour function. ##' @return The \code{ecoEwma} function returns an object of class ##' "edcc", which is a list of elements \code{optimum}, ##' \code{cost.frame}, \code{FAR} and \code{ATS}. \code{optimum} is a ##' vector with the optimum parameters and the corresponding ECH ##' value; \code{cost.frame} is a dataframe with the optimum ##' parameters and the corresponding ECH values for all given ##' \code{n}(if \code{n} is not specified, \code{cost.frame} won't be ##' returned); \code{FAR} indicates the false alarm rate during the ##' in-control time, which is calculated as lambda*(average number of ##' false alarm); \code{ATS} indicates the average time to signal ##' after the occurrence of an assignable cause, calculated as h*ARL2 ##' - tau, where tau is the expected time of occurrence of the ##' assignable cause given it occurs between the i-th and (i+1)st ##' samples. ##' The \code{echEwma} function returns the calculated ECH value only. ##' @seealso \code{\link{ecoXbar}}, \code{\link{ecoCusum}}, ##' \code{\link[spc]{xewma.arl}}, \code{\link{update.edcc}}, ##' \code{\link[stats]{optim}},\code{\link{contour}} ##' @references Weicheng Zhu, Changsoon Park (2013), {edcc: An R ##' Package for the Economic Design of the Control ##' Chart}. \emph{Journal of Statistical Software}, 52(9), ##' 1-24. \url{http://www.jstatsoft.org/v52/i09/} ##' ##' Lorenzen and Vance (1986). The economic design of ##' control charts: a unified approach, \emph{Technometrics}, 28. 3-10. ##' @examples ##' #Douglas (2009). Statistical quality control: a modern introduction, sixth edition, p470. ##' ## Set w from 0.1 to 1 by 0.1 to catch the trend. ##' ecoEwma(w=seq(0.1,1,by=0.1),P0=110,P1=10,Cf=50) ##' #yy = ecoEwma(h = seq(.7,1,by=.01), w = seq(0.8,1,by=.01),k = seq(2.9,3.3, by = 0.01), n = 4:5, P0 = 110, P1 = 10, Cf = 50, contour.plot = TRUE) ##' ##' ##$optimum
##' ##Optimum h Optimum k Optimum n Optimum w       ECH
##' ##  0.81000   2.99000   5.00000   0.95000  10.36482
##' #contour(yy)
##' @export
ecoEwma <- function( h, w, k, n, delta = 2,lambda = .05, P0 = NULL, P1 = NULL,C0 = NULL,C1 = NULL, Cr = 25, Cf = 10,T0 = 0.0167,Tc = 1, Tf = 0, Tr = 0, a = 1, b = .1,d1=1,d2=1, nlevels=30, sided="two",par=NULL,contour.plot=FALSE,call.print=TRUE,...){
if(!is.null(par)){
if(!is.vector(par)) stop("par should be a vector of length 2 or 3!") else
if(!(length(par)==2 || length(par)==3)) stop("par should be a vector of length 2 or 3!") else
if(length(par)==3){
cost.frame <- NULL
for(wi in w){
opt <- optim(par,.echEwma1,w=wi,delta = delta,lambda = lambda, P0 = P0, P1 = P1, C0 = C0, C1 = C1, Cr = Cr, Cf = Cf,T0 = T0, Tc = Tc,Tf = Tf,Tr = Tr, a = a, b = b, d1 = d1, d2 = d2,sided = sided,...)
cost.frame <- rbind(cost.frame,c(opt$par,wi,opt$value))
}
optimum <- cost.frame[which(cost.frame[,5]==min(cost.frame[,5])),]
names(optimum) <- c("Optimum h","Optimum k","Optimum n", "Optimum w","ECH")
colnames(cost.frame) <- c("Optimum h","Optimum k","Optimum n", "Optimum w","ECH")
rownames(cost.frame) <- rep("",length(w))
optEwma <- list(optimum=optimum, cost.frame=cost.frame)
} else{
cost.frame <- NULL
for(wi in w){
for(nk in n){
opt <- optim(par,.echEwma2,w=wi,n=nk,delta = delta,lambda = lambda, P0 = P0, P1 = P1, C0 = C0, C1 = C1, Cr = Cr, Cf = Cf,T0 = T0, Tc = Tc,Tf = Tf,Tr = Tr, a = a, b = b, d1 = d1, d2 = d2,sided = sided,...)
cost.frame <- rbind(cost.frame,c(opt$par,nk,wi,opt$value))
}
}
optimum <- cost.frame[which(cost.frame[,5]==min(cost.frame[,5])),]
names(optimum) <- c("Optimum h","Optimum k","Optimum n", "Optimum w","ECH")
colnames(cost.frame) <- c("Optimum h","Optimum k","Optimum n","Optimum w","ECH")
rownames(cost.frame) <- rep("",length(n)*length(w))
optEwma <- list(optimum=optimum, cost.frame=cost.frame)
}
}else
if( missing(w))
stop("At least 'w' should be given!") else
if(missing(h) && missing(k) && missing(n)){
cost.frame <- NULL
for(wi in w){
opt <- optim(c(1,3,5),.echEwma1,w=wi,delta = delta,lambda = lambda, P0 = P0, P1 = P1, C0 = C0, C1 = C1, Cr = Cr, Cf = Cf,T0 = T0, Tc = Tc,Tf = Tf,Tr = Tr, a = a, b = b, d1 = d1, d2 = d2,sided = sided,...)
cost.frame <- rbind(cost.frame,c(opt$par,wi,opt$value))
}
optimum <- cost.frame[which(cost.frame[,5]==min(cost.frame[,5])),]
names(optimum) <- c("Optimum h","Optimum k","Optimum n", "Optimum w","ECH")
colnames(cost.frame) <- c("Optimum h","Optimum k","Optimum n", "Optimum w","ECH")
rownames(cost.frame) <- rep("",length(w))
optEwma <- list(optimum=optimum, cost.frame=cost.frame)
}else
if(missing(h) && missing(k)){
cost.frame <- NULL
for(wi in w){
for(nk in n){
opt <- optim(c(1,3),.echEwma2,w=wi,n=nk,delta = delta,lambda = lambda, P0 = P0, P1 = P1, C0 = C0, C1 = C1, Cr = Cr, Cf = Cf,T0 = T0, Tc = Tc,Tf = Tf,Tr = Tr, a = a, b = b, d1 = d1, d2 = d2,sided = sided,...)
cost.frame <- rbind(cost.frame,c(opt$par,nk,wi,opt$value))
}
}
optimum <- cost.frame[which(cost.frame[,5]==min(cost.frame[,5])),]
names(optimum) <- c("Optimum h","Optimum k","Optimum n", "Optimum w","ECH")
colnames(cost.frame) <- c("Optimum h","Optimum k","Optimum n","Optimum w","ECH")
rownames(cost.frame) <- rep("",length(n)*length(w))
optEwma <- list(optimum=optimum,cost.frame=cost.frame)
}else{
cost.frame <- NULL
opt.mat <- array(NA,c(length(h),length(w),length(k)))
for(i in 1:length(k))
for(j in 1:length(w))
opt.mat[,j,i] = sapply(h,FUN=echEwma,w=w[j],k=k[i],n=n[1],delta = delta,lambda = lambda, P0 = P0, P1 = P1,C0 = C0,C1 = C1, Cr = Cr, Cf = Cf,T0 = T0,Tc = Tc, Tf = Tf, Tr = Tr, a = a, b = b,d1=d1,d2=d2,sided=sided)
aa <- which(opt.mat==min(opt.mat),arr.ind=TRUE)
opt.ech <- min(opt.mat)
cost.frame <- rbind(cost.frame,c(h[aa[1,1]],k[aa[1,3]],n[1],w[aa[1,2]],min(opt.mat)))
for(nk in n[-1]){
mat <- array(NA,c(length(h),length(w),length(k)))
for(i in 1:length(k))
for(j in 1:length(w))
mat[,j,i] = sapply(h,FUN=echEwma,w=w[j],k=k[i],n=nk,delta = delta,lambda = lambda, P0 = P0, P1 = P1,C0 = C0,C1 = C1, Cr = Cr, Cf = Cf,T0 = T0,Tc = Tc, Tf = Tf, Tr = Tr, a = a, b = b,d1=d1,d2=d2,sided=sided)
aa <- which(mat==min(mat),arr.ind=TRUE)
cost.frame <- rbind(cost.frame,c(h[aa[1,1]],k[aa[1,3]],nk, w[aa[1,2]],ech <- min(mat)))
if(ech < opt.ech){
opt.ech <- ech; opt.mat <- mat
}
}
optimum <- cost.frame[which(cost.frame[,5]==min(cost.frame[,5])),]
names(optimum) <- c("Optimum h","Optimum k","Optimum n", "Optimum w","ECH")
colnames(cost.frame) <- c("Optimum h","Optimum k","n","Optimum w","ECH")
rownames(cost.frame) <- rep("",length(n))
mat.hk <- opt.mat[,which(w==optimum[4]),]
mat.hw <- opt.mat[,,which(k==optimum[2])]
mat.wk <- opt.mat[which(h==optimum[1]),,]
if(contour.plot){
if(call.print){
ca <- match.call()
ca[[1]] <- NULL
par(mfrow=c(2,2),mar=c(5.1,4.1,1.1,1.1),oma=c(2,0,3,0))
contour(h,k,mat.hk,xlab="h",ylab="k",nlevels=nlevels)
points(optimum[1],optimum[3],pch=3)
contour(h,w,mat.hw,xlab="h",ylab="w",nlevels=nlevels)
points(optimum[1],optimum[2],pch=3)
contour(w,k,mat.wk,xlab="w",ylab="k",nlevels=nlevels)
points(optimum[2],optimum[3],pch=3)
mtext(sprintf('Opt h=%s   Opt w=%s   Opt k=%s   n=%s   ECH=%s',optimum[1], optimum[4], optimum[2], optimum[3], round(optimum[5], digits=4)),outer=T,line=.3,side=1)
title(strwrap(paste("ecoEwma(",paste(names(ca),"=",unlist(ca),sep="",collapse=", "),")")),outer=T)
} else{
par(mfrow=c(2,2),mar=c(5.1,4.1,1.1,1.1),oma=c(2,0,1,0))
contour(h,k,mat.hk,xlab="h",ylab="k",nlevels=nlevels)
points(optimum[1],optimum[2],pch=3)
contour(h,w,mat.hw,xlab="h",ylab="w",nlevels=nlevels)
points(optimum[1],optimum[4],pch=3)
contour(w,k,mat.wk,xlab="w",ylab="k",nlevels=nlevels)
points(optimum[4],optimum[2],pch=3)
mtext(sprintf('Opt h=%s   Opt w=%s   Opt k=%s   n=%s   ECH=%s',optimum[1], optimum[4], optimum[2],optimum[3],round(optimum[5],digits=4)),outer=T,line=.3,side=1)
}
}
optEwma <- list(optimum=optimum,cost.frame=cost.frame,opt.mat=list(mat.hk=mat.hk,mat.hw=mat.hw,mat.wk=mat.wk))
}
opth <- as.numeric(optEwma$optimum[1]) optk <- as.numeric(optEwma$optimum[2])
optn <- as.numeric(optEwma$optimum[3]) optw <- as.numeric(optEwma$optimum[4])
delta.std <- sqrt(as.numeric(optn))*delta
ARL1 <- as.numeric(xewma.arl(optw,optk,0,sided=sided))
ARL2 <- as.numeric(xewma.arl(optw,optk,delta.std,sided=sided))
s <- 1/(exp(lambda*opth)-1)
tau <- (1-(1+lambda*opth)*exp(-lambda*opth))/(lambda*(1-exp(-lambda*opth)))
optEwma$FAR <- lambda*s/ARL1 optEwma$ATS <- ARL2*opth - tau
optEwma$call <- match.call() return(structure(optEwma,class="edcc")) } ##' @rdname ecoEwma ##' @export # calculate the ECH for given parameters echEwma <- function(h,w,k,n,delta = 2,lambda = .05, P0 = NULL, P1 = NULL,C0 = NULL,C1 = NULL, Cr = 25, Cf = 10,T0 = 0.0167,Tc = 1, Tf = 0, Tr = 0, a = 1, b = .1,d1=1,d2=1,sided="two"){ delta.std <- sqrt(n)*delta #standardization fordelta ARL1 <- as.numeric(xewma.arl(w,k,0,sided=sided)) ARL2 <- as.numeric(xewma.arl(w,k,delta.std,sided=sided)) tau <- (1-(1+lambda*h)*exp(-lambda*h))/(lambda*(1-exp(-lambda*h))) s <- 1/(exp(lambda*h)-1) if(!is.null(P0)&!is.null(P1)){ ECT <- 1/lambda+(1-d1)*s*Tf/ARL1 - tau + n*T0 + h*ARL2 + Tc + Tr ECP <- P0/lambda + P1*(-tau+n*T0+h*ARL2+d1*Tc+d2*Tr) - s*Cf/ARL1 - Cr - (a+b*n)*(1/lambda-tau+n*T0+h*ARL2+d1*Tc+d2*Tr)/h ECH <- P0 - ECP/ECT }else if(!is.null(C0)&!is.null(C1)){ ECT <- 1/lambda+(1-d1)*s*Tf/ARL1 - tau + n*T0 + h*ARL2 + Tc + Tr ECC <- C0/lambda + C1*(-tau+n*T0+h*ARL2+d1*Tc+d2*Tr) + s*Cf/ARL1+Cr+(a+b*n)*(1/lambda-tau+n*T0+h*ARL2+d1*Tc+d2*Tr)/h ECH <- ECC/ECT }else stop("You should at least give a pair of value to P0,P1 or C0,C1") return(ECH) } ##' @S3method print edcc ##' @method print edcc print.edcc <- function(x,...){ xx <- x xx$opt.mat <- NULL
xx$call <- NULL class(xx) <- NULL print.default(xx,...) } ##' 'update' will update and (by default) re-fit a model. It does ##' this by extracting the call stored in the object, updating the ##' call and (by default) evaluating that call. ##' ##' S3 method for update. ##' @title Update for an "edcc" class object ##' @param object an object of "edcc" class ##' @param ... additional arguments to the call, or arguments with ##'anged values. ##' @param evaluate If true evaluate the new call else return the ##' call. ##' @return the fitted object ##' @S3method update edcc ##' @method update edcc ##' @examples x <- ecoXbar(P0=110,P1=10) ##' update(x,P0=NULL,P1=NULL,C0=10,C1=110) update.edcc <- function(object,...,evaluate=TRUE){ if (is.null(call <- object$call))
stop("need an object with call component")
extras <- match.call(expand.dots = FALSE)$... if (length(extras)) { existing <- !is.na(match(names(extras), names(call))) for (a in names(extras)[existing]) call[[a]] <- extras[[a]] if (any(!existing)) { call <- c(as.list(call), extras[!existing]) call <- as.call(call) } } if (evaluate) eval(call, parent.frame()) else call } ##' contour plot of an "edcc" class object ##' ##' S3 method of contour plot for "edcc" class object ##' @title Contour plot of an "edcc" class object ##' @param x an object of "edcc" class ##' @param call.print a logical value indicating whether the the R ##' command should be printed on the contour plot. Default is TRUE ##' @param ... arguments to be passed to contour plot, see ##' \code{\link[graphics]{contour}} for details ##' @return a contour plot ##' @S3method contour edcc ##' @method contour edcc ##' @seealso \code{\link{ecoXbar}}, \code{\link{ecoCusum}}, ##' \code{\link{ecoEwma}}, \code{\link{update.edcc}}, ##' \code{\link[graphics]{contour}} ##' @examples ##' x <- ecoXbar(h=seq(0.7,0.9,by=.01),L=seq(2.8,3.3,by=.01),n=4:6,P0=110, ##' P1=10,nlevels=50,contour.plot=TRUE) ##' contour(x,nlevels=20,lty=2,col=2,call.print=FALSE) contour.edcc <- function(x,call.print=TRUE,...){ if (is.null(call <- x$call))
stop("need an object with call component")
extras <- match.call(expand.dots = FALSE)$... if(call[[1]]==quote(ecoXbar)){ h <- call[["h"]] L <- call[["L"]] nlevels <- call[["nlevels"]] if(is.null(nlevels)) nlevels <- 30 if(missing(call.print)){ call.print <- call[["call.print"]] if(is.null(call.print)) call.print <- TRUE } opt.mat <- x$opt.mat
call1 <- list(quote(contour.default),h,L,quote(opt.mat),nlevels=nlevels,xlab="h",ylab="L")
if (length(extras)) {
existing <- !is.na(match(names(extras), c("nlevels")))
call1[["nlevels"]] <- extras[[names(extras)[existing]]]
if (any(!existing)) {
call1 <- c(call1, extras[!existing])
}
}
call1 <- as.call(call1)
if(call.print){
par(mar=c(6.1,4.1,4.1,2.1))
eval(call1)
call[[1]] <- NULL
title(strwrap(paste("ecoXbar(",paste(names(call),"=",unlist(call),sep="",collapse=", "),")")))
} else{
par(mar=c(6.1,4.1,2.1,2.1))
eval(call1)
}
optimum <- x$optimum mtext(sprintf('Opt h=%s Opt L=%s n=%s ECH=%s',optimum[1], optimum[2], optimum[3],round(optimum[4],digits=4)),side=1,line=4.5) points(optimum[1],optimum[2],pch=3) }else if(call[[1]]==quote(ecoCusum)){ h <- call[["h"]] H <- call[["H"]] nlevels <- call[["nlevels"]] if(is.null(nlevels)) nlevels <- 30 if(missing(call.print)){ call.print <- call[["call.print"]] if(is.null(call.print)) call.print <- TRUE } opt.mat <- x$opt.mat
call1 <- list(quote(contour.default),h,H,quote(opt.mat),nlevels=nlevels,xlab="h",ylab="H")
if (length(extras)) {
existing <- !is.na(match(names(extras), c("nlevels")))
call1[["nlevels"]] <- extras[[names(extras)[existing]]]
if (any(!existing)) {
call1 <- c(call1, extras[!existing])
}
}
call1 <- as.call(call1)
if(call.print){
par(mar=c(6.1,4.1,4.1,2.1))
eval(call1)
call[[1]] <- NULL
title(strwrap(paste("ecoCusum(",paste(names(call),"=",unlist(call),sep="",collapse=", "),")")))
} else{
par(mar=c(6.1,4.1,2.1,2.1))
eval(call1)
}
optimum <- x$optimum mtext(sprintf('Opt h=%s Opt H=%s n=%s ECH=%s',optimum[1], optimum[2], optimum[3],round(optimum[4],digits=4)),side=1,line=4.5) points(optimum[1],optimum[2],pch=3) }else if(call[[1]]==quote(ecoEwma)){ h <- call[["h"]] w <- call[["w"]] k <- call[["k"]] nlevels <- call[["nlevels"]] if(is.null(nlevels)) nlevels <- 30 if(missing(call.print)){ call.print <- call[["call.print"]] if(is.null(call.print)) call.print <- TRUE } opt.mat <- x$opt.mat
call1 <- list(quote(contour.default),h,k,quote(opt.mat$mat.hk),nlevels=nlevels,xlab="h",ylab="k") call2 <- list(quote(contour.default),h,w,quote(opt.mat$mat.hw),nlevels=nlevels,xlab="h",ylab="w")
call3 <- list(quote(contour.default),w,k,quote(opt.mat$mat.wk),nlevels=nlevels,xlab="w",ylab="k") if (length(extras)) { existing <- !is.na(match(names(extras), c("nlevels"))) call1[["nlevels"]] <- extras[[names(extras)[existing]]] call2[["nlevels"]] <- extras[[names(extras)[existing]]] call3[["nlevels"]] <- extras[[names(extras)[existing]]] if (any(!existing)) { call1 <- c(call1, extras[!existing]) call2 <- c(call2, extras[!existing]) call3 <- c(call3, extras[!existing]) } } call1 <- as.call(call1) call2 <- as.call(call2) call3 <- as.call(call3) optimum <- x$optimum
if(call.print){
par(mfrow=c(2,2),mar=c(5.1,4.1,1.1,1.1),oma=c(2,0,3,0))
eval(call1)
points(optimum[1],optimum[3],pch=3)
eval(call2)
points(optimum[1],optimum[2],pch=3)
eval(call3)
points(optimum[2],optimum[3],pch=3)
call[[1]] <- NULL
title(strwrap(paste("ecoEwma(",paste(names(call),"=",unlist(call),sep="",collapse=", "),")")),outer=T)
} else{
par(mfrow=c(2,2),mar=c(5.1,4.1,1.1,1.1),oma=c(2,0,1,0))
eval(call1)
points(optimum[1],optimum[3],pch=3)
eval(call2)
points(optimum[1],optimum[2],pch=3)
eval(call3)
points(optimum[2],optimum[3],pch=3)
}
mtext(sprintf('Opt h=%s   Opt w=%s   Opt k=%s   n=%s   ECH=%s',optimum[1], optimum[2], optimum[3],optimum[4],round(optimum[5],digits=4)),outer=T,line=.3,side=1)
}
}


Try the edcc package in your browser

Any scripts or data that you put into this service are public.

edcc documentation built on May 2, 2019, 1:40 p.m.