<!--
%\VignetteEngine{knitr::knitr}
%\VignetteIndexEntry{Extensions}
-->
Extensions 
===========================================
Extensions to other data models and estimation procedures
-------------------------------------------------------


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.

Robust estimation
------------------
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 . 
```{r}
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:
```{r}
X <-  rnorm(100)
chartrobust <- new("SPCCUSUM",model=model)
SPCproperty(data=X,nrep=50,property="calARL",
            chart=chartrobust,params=list(target=100),quiet=TRUE)
```


Parametric exponential CUSUM chart
------------------------------------------


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. 
```{r}
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. 
```{r}
ExpCUSUMchart=new("SPCCUSUM",model=SPCModelExponential(Delta=1.25))
```

```{r,echo=FALSE}
set.seed(223819)
```

The following creates some past observations, estimates that chart parameters and runs a chart for
100 steps with new observations.
```{r,fig=TRUE,fig.width=8,fig.height=4}
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). 
```{r}
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: 
```{r,fig=TRUE,fig.width=10,fig.height=4}
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)
```