Nothing
##' Calculates a sequence of Cusum statistics
##'
##' Calculates a sequence of one-sided upper Cusum statistics given the reference value and
##' the control limit.
##'
##' @details Cusum is assumed to be of the form: \emph{C[i] = max(0, C[i-1] + X[i] - k)},
##' where the signal occurs when \emph{C[i] > h}. Note that \code{X} can be the Cusum scores, or weights,
##' given by the log-likelihood ratio, in which case \code{k = 0} would make sense.
##'
##' @export
##' @param X A numeric vector.
##'
##' @param k The reference value.
##'
##' @param h The upper control limit.
##'
##' @param initial The starting value of the Cusum (\emph{C[0]}).
##'
##' @param reset Logical indicating whether the Cusum is reset to 0 after crossing the control limit.
##'
##' @param x Object of class \code{cusum}
##'
##' @param object Object of class \code{cusum}
##'
##' @param \dots Additional arguments to \code{\link{print.default}} or \code{\link{plot.default}}. Ignored by the \code{signal} method.
##'
##' @return A object of class \code{cusum}, which is a vector of the Cusum statistics, along with the following attributes:
##' \code{X}, \code{k}, \code{h}, \code{initial}, and \code{reset} (which correspond to the original arguments provided to
##' the function) and \code{resetCounter}, a vector of integers corresponding to \code{cusum} that indicates when the
##' Cusum resets.
##'
##' @references Hawkins DM and Olwell DH. (1998) Cumulative Sum Charts and Charting for Quality Improvement. Springer.
##'
##' @examples
##' y <- cusum(rnorm(50), 0.2, 2)
##' y
##'
##' # Plot the cusum
##' plot(y)
##'
##' # Show the indexes where the chart signaled
##' signal(y)
##'
##' # A look at the attributes
##' attributes(y)
# A Wrapper for the c method 'cusum'
# Calculates the Cusum for a vector of data
# Cusum of the form: C[i] = max(0, C[i-1] + X[i] - k)
# Signal occurs when C[i] > h
cusum <- function(X, k, h, initial = 0, reset = TRUE) {
# Verify inputs
stopifnotMsg(# X
is.numeric(X) & (length(X) > 0),
"'X' must be numeric and have positive length",
# k
if (is.numeric(k) & (length(k) == 1)) {
k >= 0
} else FALSE,
"'k' must be a single, non-negative number",
# h
if (is.numeric(h) & (length(h) == 1)) {
h > 0
} else FALSE,
"'h' must be a single, positive number",
# initial
is.numeric(initial) & (length(initial) == 1),
"'initial' must be a single number",
# reset
is.logical(reset) & (length(reset) == 1),
"'reset' must be TRUE or FALSE")
if (initial > h) {
warning("Cusum initialized above the control limit")
}
# Calculate the cusum
cusum <- .C("cusum_c",
X = as.double(X),
k = as.double(k),
h = as.double(h),
initial = as.double(initial),
reset = as.integer(reset),
upper = 1L, # upper set to TRUE
n = as.integer(length(X)),
cusum = double(length(X)),
# stagger = double(length(X)),
resetCounter = integer(length(X)))[c(1:5, 8, 9)]
# Add in names if needed
if (!is.null(names(X))) {
names(cusum$cusum) <- names(cusum$resetCounter) <- names(cusum$X) <- names(X)
}
out <- cusum$cusum
class(out) <- c("cusum", class(out))
attributes(out) <- c(attributes(out),
cusum[c("X", "k", "h", "initial")],
list(reset = reset),
cusum["resetCounter"])
# Return the list
return(out)
} # cusum
##' @method print cusum
##'
##' @describeIn cusum Prints the \code{cusum} object by only showing the Cusum statistics and suppressing the attributes.
##'
##' @export
print.cusum <- function(x, ...) {
printWithoutAttributes(x, ...)
} # print.cusum
##' @method plot cusum
##'
##' @describeIn cusum Plots the \code{cusum} object.
##'
##' @param indexes A vector of indexes that select the elements of the cusum statistics that will be plotted.
##'
##' @param emphOOC A logical indicating whether out of control points should be emphasized in red.
##'
##' @export
plot.cusum <- function(x, indexes = NULL, emphOOC = TRUE, ...) {
# Check inputs
stopifnotMsg(if (!is.null(indexes)) {
if (is.numeric(indexes) & (length(indexes) <= length(x)) & (length(indexes) > 0)) {
all(indexes %in% 1:length(x))
} else FALSE
} else TRUE,
"'indexes' must be whole numbers in the set '1:length(x)', or NULL",
is.logical(emphOOC) & (length(emphOOC) == 1),
"'emphOOC' must be TRUE or FALSE")
# If indexes are NULL, set to all possible to make selection easier
if (is.null(indexes)) {
indexes <- 1:length(x)
}
# Get the set of data points that will be plotted
xvals <- indexes
yvals <- x[indexes]
# Get the control limit, as it is used often
h <- attributes(x)$h
# Set the default plot args
defaultArgs <- list(x = indexes,
y = yvals,
xlab = "Index",
ylab = "Cusum Statistic",
ylim = range(yvals, h),
type = "b",
col = "Blue",
pch = 1,
cex = 1)
# Blend in the supplied plot arguments from ...
finalArgs <- blendArgs(defaultArgs, ...)
# Create the plot if no reset
if (!attributes(x)$reset) {
do.call(plot, finalArgs)
}
else {
# Make a blank plot and then add in the lines for each cusum run
finalArgsReset <- finalArgs
finalArgsReset$type <- "n"
do.call(plot, finalArgsReset)
# Add a connected line segment for each run
cnt <- attributes(x)$resetCounter[indexes]
for (i in unique(cnt)) {
runInd <- cnt == i
# Add data points first if they were requested
do.call(lines, c(list(x = indexes[runInd], y = yvals[runInd]), finalArgs[c("type", "col", "pch", "cex")]))
}
}
# Add in the control limit
abline(h = h, col = "Red")
# Color the signal points red
if (emphOOC) {
signalInd <- which(yvals > h)
points(indexes[signalInd], yvals[signalInd], col = "Red", pch = finalArgs$pch, cex = finalArgs$cex)
}
} # plot.cusum
##' @method signal cusum
##'
##' @describeIn cusum Prints the indexes in a \code{cusum} object that exceed the control limit
##'
##' @export
signal.cusum <- function(object, ...) {
which(object > attributes(object)$h)
} # signal.cusum
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.