# CUSUM function in qcc package

I would like to know the method of calculation for the positive and negative cusum's which is used in the CUSUM function part of the qcc package.
By usuing the method for tabular CUSUM provided in the book "Introduction to statistical quality control 6th edition" by Douglas C.Montgomery which consistes of using the following formulas:

But was unable to obtain the same results as the ones provided by the cusum function.

When using the cusum (qcc)

cusum(data, center = 0.0122, std.dev = 0.006, decision.interval = 4, se.shift = 0.5)

Samples

Group [,1]
1 0.030000000
2 0.014277869
3 0.011239112
4 0.010664081
5 0.008917331
6 0.010301526
7 0.009042467
8 0.010976948
9 0.009000000
10 0.010000000
11 0.010300000

cusumchart[["statistics"]]
1 2 3 4 5 6 7 8 9 10 11
0.030000000 0.014277869 0.011239112 0.010664081 0.008917331 0.010301526 0.009042467 0.010976948 0.009000000 0.010000000 0.010300000
cusumchart[["pos"]]
[1] 2.7166667 2.8129782 2.4028302 1.8968438 1.0997323 0.5333201 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
cusumchart[["neg"]]
[1] 0.000000000 0.000000000 0.000000000 -0.005986428 -0.303097912 -0.369510173 -0.645765623 -0.599607555 -0.882940888 -0.999607555 -1.066274222
'''

You can get the code for the function by typing cusum at the console prompt. Here is what I got.

> cusum
function (data, sizes, center, std.dev, head.start = 0, decision.interval = 5,
se.shift = 1, data.name, labels, newdata, newsizes, newlabels,
plot = TRUE, ...)
{
call <- match.call()
if (missing(data))
stop("'data' argument is not specified")
if (missing(data.name))
data.name <- deparse(substitute(data))
data <- data.matrix(data)
if (missing(sizes)) {
sizes <- apply(data, 1, function(x) sum(!is.na(x)))
}
else {
if (length(sizes) == 1)
sizes <- rep(sizes, nrow(data))
else if (length(sizes) != nrow(data))
stop("sizes length doesn't match with data")
}
if (decision.interval <= 0)
stop("decision.interval must be positive")
stop("head.start must be non-negative and less than decision.interval")
type <- if (any(sizes == 1))
"xbar.one"
else "xbar"
if (ncol(data) == 1 & any(sizes > 1) & missing(std.dev))
stop("sizes larger than 1 but data appears to be single samples. In this case you must provide also the std.dev")
if (missing(labels)) {
if (is.null(rownames(data)))
labels <- 1:nrow(data)
else labels <- rownames(data)
}
stats <- paste("stats.", type, sep = "")
if (!exists(stats, mode = "function"))
stop(paste("function", stats, "is not defined"))
stats <- do.call(stats, list(data, sizes))
statistics <- stats\$statistics
if (missing(center))
center <- stats\$center
sd <- paste("sd.", type, sep = "")
if (!exists(sd, mode = "function"))
stop(paste("function", sd, "is not defined!"))
if (missing(std.dev)) {
std.dev <- switch(type, xbar = {
if (any(sizes > 25)) "RMSDF" else "UWAVE-R"
}, NULL)
std.dev <- do.call(sd, list(data, sizes, std.dev))
}
else {
if (is.character(std.dev)) {
std.dev <- do.call(sd, list(data, sizes, std.dev))
}
else {
if (!is.numeric(std.dev))
stop("if provided the argument 'std.dev' must be a method available or a numerical value. See help(qcc).")
}
}
names(statistics) <- rownames(data) <- labels
names(dimnames(data)) <- list("Group", "Samples")
object <- list(call = call, type = "cusum", data.name = data.name,
data = data, statistics = statistics, sizes = sizes,
center = center, std.dev = std.dev)
if (!missing(newdata)) {
newdata.name <- deparse(substitute(newdata))
newdata <- data.matrix(newdata)
if (missing(newsizes)) {
newsizes <- apply(newdata, 1, function(x) sum(!is.na(x)))
}
else {
if (length(newsizes) == 1)
newsizes <- rep(newsizes, nrow(newdata))
else if (length(newsizes) != nrow(newdata))
stop("newsizes length doesn't match with newdata")
}
stats <- paste("stats.", type, sep = "")
if (!exists(stats, mode = "function"))
stop(paste("function", stats, "is not defined"))
newstats <- do.call(stats, list(newdata, newsizes))\$statistics
if (missing(newlabels)) {
if (is.null(rownames(newdata))) {
start <- length(statistics)
newlabels <- seq(start + 1, start + length(newstats))
}
else {
newlabels <- rownames(newdata)
}
}
names(newstats) <- newlabels
object\$newstats <- newstats
object\$newdata <- newdata
object\$newsizes <- newsizes
object\$newdata.name <- newdata.name
statistics <- c(statistics, newstats)
sizes <- c(sizes, newsizes)
}
n <- length(statistics)
z <- (statistics - center)/(std.dev/sqrt(sizes))
ldb <- -decision.interval
udb <- decision.interval
z.f <- z - se.shift/2
cusum.pos <- rep(NA, n)
cusum.pos[1] <- max(0, head.start + z.f[1])
for (i in 2:n) cusum.pos[i] <- max(0, cusum.pos[i - 1] +
z.f[i])
z.f <- z + se.shift/2
cusum.neg <- rep(NA, n)
cusum.neg[1] <- max(0, head.start - z.f[1])
for (i in 2:n) cusum.neg[i] <- max(0, cusum.neg[i - 1] -
z.f[i])
cusum.neg <- -cusum.neg
violations <- list(lower = which(cusum.neg < ldb), upper = which(cusum.pos >
udb))
object\$type <- "cusum"
object\$pos <- cusum.pos
object\$neg <- cusum.neg
object\$decision.interval <- decision.interval
object\$se.shift <- se.shift
object\$violations <- violations
class(object) <- "cusum.qcc"
if (plot)
plot(object, ...)
return(object)
}
<bytecode: 0x000001caa0dfec48>
<environment: namespace:qcc>

This topic was automatically closed 42 days after the last reply. New replies are no longer allowed.

If you have a query related to it or one of the replies, start a new topic and refer back with a link.