The following is an example of an application to a basic CUSUM chart, assuming that all observations are normally distributed.

Based on $n$ past in-control observations $X_{-n},\dots,X_{-1}$, the in-control mean can be estimated by $\hat \mu = \frac{1}{n}\sum_{i=-n}^{-1} X_i$ and the in-control variance by $\hat \sigma^2=\frac{1}{n-1}\sum_{i=-n}^{-1} (X_i-\hat \mu)^2$. For new observations $X_1,X_2,\dots$, a CUSUM chart based on these estimated parameters is then defined by $$ S_0=0, \quad S_t=\max\left(0,\frac{S_{t-1}+X_t-\hat \mu-\Delta/2}{\hat \sigma}\right). $$

set.seed(12381900)

The following generates a data set of past observations (replace this with your observed past data).

X <- rnorm(250)

par(mar=c(4,5,0,0)) plot(-(250:1),X,xlab="t",ylab=expression(X[t]))

Next, we initialise the chart and compute the estimates needed for running the chart - in this case $\hat \mu$ and $\hat \sigma$.

library(spcadjust) chart <- new("SPCCUSUM",model=SPCModelNormal(Delta=1)); xihat <- xiofdata(chart,X) str(xihat)

We now compute a threshold that with roughly 90% probability results in an average run length of at least 100 in control. This is based on parametric resampling assuming normality of the observations.

cal <- SPCproperty(data=X,nrep=50, property="calARL",chart=chart,params=list(target=100),quiet=TRUE) cal

The number of bootstrap replications (the argument nrep) shoud be increased for real applications. Use the parallell option to speed up the bootstrap by parallel processing.

Next, we run the chart with new observations that are in-control.

newX <- rnorm(100) S <- runchart(chart, newdata=newX,xi=xihat)

Then we plot the data and the chart.

par(mfrow=c(1,2),mar=c(4,5,0.1,0.1)) plot(newX,xlab="t") plot(S,ylab=expression(S[t]),xlab="t",type="b",ylim=range(S,cal@res+1,cal@raw)) lines(c(0,100),rep(cal@res,2),col="red") lines(c(0,100),rep(cal@raw,2),col="blue") legend("topleft",c("Adjusted Threshold","Unadjusted Threshold"),col=c("red","blue"),lty=1)

In the next example, the chart is run with data which are out-of-control from time 51 and onwards.

newX <- rnorm(100,mean=c(rep(0,50),rep(1,50))) S <- runchart(chart, newdata=newX,xi=xihat)

par(mfrow=c(1,2),mar=c(4,5,0.1,0.1)) plot(newX,xlab="t") plot(S,ylab=expression(S[t]),xlab="t",type="b",ylim=pmin(range(S,cal@res,cal@raw),15)) lines(c(0,100),rep(cal@res,2),col="red") lines(c(0,100),rep(cal@raw,2),col="blue") legend("topleft",c("Adjusted Threshold","Unadjusted Threshold"),col=c("red","blue"),lty=1)

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

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.