The framework provided in this package can be easily extended to work with other charts, other data models and other estimation procedures. The following are some examples of extensions.
Consider a CUSUM chart for observations with a normal distribution with unknown mean and standard deviation. We now illustrate how to change the estimation procedure to use the median and the mean absolute deviation (MAD) instead of the mean and the sample standard deviation. To achieve this we only need to modify one function of the existing data model, the method Pofdata that estimates the parameters .
library(spcadjust) model <- SPCModelNormal(Delta=1) model$Pofdata model$Pofdata <- function(data){ list(mu= median(data), sd= mad(data), m=length(data)) }
Properties of this chart can then be computed as before:
X <- rnorm(100) chartrobust <- new("SPCCUSUM",model=model) SPCproperty(data=X,nrep=50,property="calARL", chart=chartrobust,params=list(target=100),quiet=TRUE)
In this example we construct a CUSUM chart for exponentially
distributed data with unknown rate $\lambda$. Such a CUSUM can be constructed by defining the
updates as the log-likelihood ratio between an out of control model
with rate $\Delta\lambda$ and an in-control model with
rate $\lambda$. I.e.
$$
S_0=0, \quad S_t=\max\left(0,S_{t-1}+R_t \right)
$$
where
$$
R_t=\log\left(
\frac{\lambda\Delta\exp(-\lambda\Delta X_t)}{\lambda\exp(-\lambda X_t)}
\right)
=\log(\Delta)-\lambda(\Delta-1)X_t
$$
defines the updates.
We choose to use parametric bootstrapping, and
the rate is estimated from past data
$X_{-n},\dots,X_{-1}$ by the maximum likelihood estimator
$\hat\lambda=n/\sum_{i=1}^nX_{-i}$.
To implement this CUSUM we need to define a new data model object of class SPCDataModel. The following code does this.
SPCModelExponential=function(Delta=1.25){ structure( list( Pofdata=function(data){ list(lambda=1/mean(data), n=length(data)) }, xiofP=function(P) P$lambda, resample=function(P) rexp(P$n,rate=P$lambda), getcdfupdates=function(P, xi) { function(x){ if(Delta<1) pmax(0,1-exp(-P$lambda*(x-log(Delta))/(xi*(1-Delta)))) else pmin(1,exp(-P$lambda*(log(Delta)-x)/(xi*(Delta-1)))) } }, updates=function(xi,data) log(Delta)-xi*(Delta-1)*data), class="SPCDataModel") }
In the above code Pofdata estimates the probability model, xiofP computes the parameter needed to run the chart, resample specifies the parametric bootstrap and getcdfupdates specifies the cdf of the updates when the parameter estimate xi is used in the updates.
To create a CUSUM chart with this data model we first initialise the chart.
ExpCUSUMchart=new("SPCCUSUM",model=SPCModelExponential(Delta=1.25))
set.seed(223819)
The following creates some past observations, estimates that chart parameters and runs a chart for 100 steps with new observations.
X <- rexp(1000) plot(runchart(ExpCUSUMchart, newdata=rexp(100),xi=xiofdata(ExpCUSUMchart,X)), ylab=expression(S[t]),xlab="t",type="b")
The following computes various properties of the chart. For real applications increase the number of bootstrap replications (the argument nrep) and consider using parallel processing (the argument parallel).
SPCproperty(data=X,nrep=100,property="hitprob", chart=ExpCUSUMchart,params=list(threshold=3,nsteps=100), covprob=c(0.5,0.9),quiet=TRUE) SPCproperty(data=X,nrep=100,property="ARL", chart=ExpCUSUMchart,params=list(threshold=3),covprob=c(0.5,0.9),quiet=TRUE) SPCproperty(data=X,nrep=100,property="calARL",chart=ExpCUSUMchart, params=list(target=1000),covprob=c(0.5,0.9),quiet=TRUE)
To make a CUSUM plot with thresholds:
cal <- SPCproperty(data=X,nrep=1000,property="calARL",chart=ExpCUSUMchart, params=list(target=1000),quiet=TRUE,parallel=1) S <- runchart(ExpCUSUMchart, newdata=rexp(100),xi=xiofdata(ExpCUSUMchart,X)) par(mfrow=c(1,1),mar=c(4,5,0,0)) 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)
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.