Percent defective for normal samples

Description

Density, distribution function and quantile function for the sample percent defective calculated on normal samples with mean equal to mu and standard deviation equal to sigma.

Usage

dphat(x, n, mu=0, sigma=1, type="known", LSL=-3, USL=3, nodes=30)

pphat(q, n, mu=0, sigma=1, type="known", LSL=-3, USL=3, nodes=30)

qphat(p, n, mu=0, sigma=1, type="known", LSL=-3, USL=3, nodes=30)

Arguments

Details

Bruhn-Suhr/Krumbholz (1990) derived the cumulative distribution function of the sample percent defective calculated on normal samples to applying them for a new variables sampling plan. These results were heavily used in Krumbholz/Zöller (1995) for Shewhart and in Knoth/Steinmetz (2013) for EWMA control charts. For algorithmic details see, essentially, Bruhn-Suhr/Krumbholz (1990). Two design variants are treated: The simple case, type="known", with known normal variance and the presumably much more relevant and considerably intricate case, type="estimated", where both parameters of the normal distribution are unknown. Basically, given lower and upper specification limits and the normal distribution, one estimates the expected yield based on a normal sample of size n.

Value

Returns vector of pdf, cdf or qf values for the statistic phat.

Author(s)

Sven Knoth

References

M. Bruhn-Suhr and W. Krumbholz (1990), A new variables sampling plan for normally distributed lots with unknown standard deviation and double specification limits, Statistical Papers 31(1), 195-207.

W. Krumbholz and A. Zöller (1995), p-Karten vom Shewhartschen Typ für die messende Prüfung, Allgemeines Statistisches Archiv 79, 347-360.

S. Knoth and S. Steinmetz (2013), EWMA p charts under sampling by variables, International Journal of Production Research 51(13), 3795-3807.

See Also

phat.ewma.arl for routines using the herewith considered phat statistic.

Examples

# Figures 1 (c) and (d) from Knoth/Steinmetz (2013)
n      <-  5
LSL    <- -3
USL    <-  3

par(mar=c(5, 5, 1, 1) + 0.1)

p.star <- 2*pnorm( (LSL-USL)/2 ) # for p <= p.star pdf and cdf vanish

p_ <- seq(p.star+1e-10, 0.07, 0.0001) # define support of Figure 1

# Figure 1 (c)
pp_ <- pphat(p_, n)
plot(p_, pp_, type="l", xlab="p", ylab=expression(P( hat(p) <= p )),
     xlim=c(0, 0.06), ylim=c(0,1), lwd=2)
abline(h=0:1, v=p.star, col="grey")

# Figure 1 (d)
dp_ <- dphat(p_, n)
plot(p_, dp_, type="l", xlab="p", ylab="f(p)", xlim=c(0, 0.06),
     ylim=c(0,50), lwd=2)
abline(h=0, v=p.star, col="grey")

Compute ARLs of Poisson NCS-EWMA control charts

Description

Computation of the (zero-state) Average Run Length (ARL) at given Poisson mean mu.

Usage

euklid.ewma.arl(gX, gY, kL, kU, mu, y0, r0=0)

Arguments

Details

A new idea of applying EWMA smoothing to count data based on integer divison with remainders. It is highly recommended to read the corresponding paper (see below).

Value

Return single value which resemble the ARL.

Author(s)

Sven Knoth

References

A. C. Rakitzis, P. Castagliola, P. E. Maravelakis (2015), A new memory-type monitoring technique for count data, Computers and Industrial Engineering 85, 235-247.

See Also

later.

Examples

# RCM (2015), Table 12, page 243, first NCS column
gX <- 5
gY <- 24
kL <- 16
kU <- 24
mu0 <- 20
#L0 <- euklid.ewma.arl(gX, gY, kL, kU, mu0, mu0)
# should be 1219.2
 

Compute control limits of MR charts for normal data

Description

Computation of control limits of standalone MR charts.

Usage

imr.RuRl_alone(L0, N=30, qm=30, M0=12, eps=1e-3)

imr.RuRl_alone_s3(L0, N=30, qm=30, M0=12)

imr.RuRl_alone_tail(L0, N=30, qm=30, M0=12)

imr.Ru_Rlgiven(Rl, L0, N=30, qm=30, M0=12)

imr.Rl_Rugiven(Ru, L0, N=30, qm=30, M0=12)

Arguments

Details

Crowder (1987a) provided some math to determine the ARL of the so-called individual moving range (IMR) chart, which consists of the mean X chart and the standard deviation MR chart. Making the alarm threshold, M0, huge (default value here is 12) for the X chart allows us to utilize Crowder's setup for standalone MR charts. For details about the IMR numerics see imr.arl. The three different versions of imr.RuRl_alone determine limits that form an ARL unbiased design, follow the restriction Rl = 1/Ru^3 and feature equal probability tails for the MR's half-normal distribution, respectively in the order given above). The other two functions are helper routines for imr.RuRl_alone. Note that the elegant approach given in Acosta-Mejia/Pignatiello (2000) is only an approximation, because the MR series is not Markovian.

Value

Returns control limit factors (alias multiples).

Author(s)

Sven Knoth

References

S. V. Crowder (1987a) Computation of ARL for Combined Individual Measurement and Moving Range Charts, Journal of Quality Technology 19(2), 98-102.

S. V. Crowder (1987b) A Program for the Computation of ARL for Combined Individual Measurement and Moving Range Charts, Journal of Quality Technology 19(2), 103-106.

D. Radson, L. C. Alwan (1995) Detecting Variance Reductions Using the Moving Range, Quality Engineering 8(1), 165-178.

C. A. Acosta-Mejia, J. J. Pignatiello (2000) Monitoring process dispersion without subgrouping, Journal of Quality Technology 32(2), 89-102.

See Also

later.

Examples

## Radson, Alwan (1995), Table 2 (Monte Carlo based), half-normal, known parameter case
## two-sided MR-alone chart, hence the ARL results has to be decreased by 1
## Here: a large M0=12 (default of the functions above) is deployed to mimic Inf
alpha <- 0.00915
Ru <- sqrt(2) * qnorm(1-alpha/4)
Rl <- sqrt(2) * qnorm(0.5+alpha/4)
M0 <- 12
## Not run: 
ARL0 <- imr.arl(M0, Ru, 0, 1, vsided="two", Rl=Rl)
RRR1995 <- imr.RuRl_alone_tail(ARL0)
RRRs <- imr.RuRl_alone_s3(ARL0)
RRR <- imr.RuRl_alone(ARL0)
results <- rbind(c(Rl, Ru), RRR1995, RRRs, RRR)
results
## End(Not run)

Compute ARLs and control limit factors for I-MR combos in case of normal data

Description

Computation of the (zero-state) Average Run Length (ARL) at given mean mu and sigma etc.

Usage

imr.arl(M, Ru, mu, sigma, vsided="upper", Rl=0, cmode="coll", N=30, qm=30)

imr.Ru_Mgiven(M, L0, N=30, qm=30)

imr.Rl_Mgiven(M, L0, N=30, qm=30)

imr.MandRu(L0, N=30, qm=30)

imr.MandRuRl(L0, N=30, qm=30)

Arguments

Details

Crowder (1987a) provided some math to determine the ARL of the so-called individual moving range (IMR) chart. The given integral equation was approximated by a linear equation system applying trapezoidal quadratures. Interestingly, Crowder did not recognize the specific behavior of the solution for Ru >= M (which is the more common case), where the resulting function L() is constant in the central part of the considered domain. In addition, by performing collocation on two (Ru > M) or three (Ru < M) subintervals piecewise, one obtains highly accurate ARL numbers. Note that imr.MandRu yields M and Ru for the upper MR trace, whereas imr.MandRuRl provides in addition the lower factor Rl for IMR consisting of two two-sided control charts. Note that the underlying ARL unbiased condition suppresses the upper limit Ru in all considered cases so far. This is not completely surprising, because the mean chart is already quite sensitive for increases in the variance. The two functions imr.Ru_Mgiven and imr.Rl_Mgiven deliver the single upper and lower limit, respectively, if a one-sided MR design is utilized and the control lmit factor M of the mean chart is set already. Note that for Ru > 2*M, the upper MR limit is not operative anymore. If it was initially an upper MR chart, then it reduces to a single mean chart. If it was originally a two-sided MR design, then it becomes a two-sided mean/lower variance chart combo. Within the latter scheme, the mean chart signals variance increases (a well-known pattern), whereas the MR subchart delivers only decreasing variance signals. However, these simple Shewhart charts exhibit in all configurations week variance decreases detection behavior. Eventually, we should note that the scientific control chart community mainly recommends to ignore MR charts, see, for example, Vardeman and Jobe (2016), whereas standards (such as ISO), commercial SPC software and many training manuals provide the IMR scheme with completely wrong upper limits for the MR chart.

Value

Returns either the ARL or control limit factors (alias multiples).

Author(s)

Sven Knoth

References

S. V. Crowder (1987a) Computation of ARL for Combined Individual Measurement and Moving Range Charts, Journal of Quality Technology 19(2), 98-102.

S. V. Crowder (1987b) A Program for the Computation of ARL for Combined Individual Measurement and Moving Range Charts, Journal of Quality Technology 19(2), 103-106.

K. C. B. Roes, R. J. M. M. Does, Y. Schurink, Shewhart-Type Control Charts for Individual Observations, Journal of Quality Technology 25(3), 188-198.

S. E. Rigdon, E. N. Cruthis, C. W. Champ (1994) Design Strategies for Individuals and Moving Range Control Charts, Journal of Quality Technology 26(4), 274-287.

D. Radson, L. C. Alwan (1995) Detecting Variance Reductions Using the Moving Range, Quality Engineering 8(1), 165-178.

S. R. Adke, X. Hong (1997) A Supplementary Test Based on the Control Chart for Individuals, Journal of Quality Technology 29(1), 16-20.

R. W. Amin, R. A. Ethridge (1998) A Note on Individual and Moving Range Control Charts, Journal of Quality Technology 30(1), 70-74.

C. A. Acosta-Mejia, J. J. Pignatiello (2000) Monitoring process dispersion without subgrouping, Journal of Quality Technology 32(2), 89-102.

N. B. Marks, T. C. Krehbiel (2011) Design And Application Of Individuals And Moving Range Control Charts, Journal of Applied Business Research (JABR) 25(5), 31-40.

D. Rahardja (2014) Comparison of Individual and Moving Range Chart Combinations to Individual Charts in Terms of ARL after Designing for a Common “All OK” ARL, Journal of Modern Applied Statistical Methods 13(2), 364-378.

S. B. Vardeman, J. M. Jobe (2016) Statistical Methods for Quality Assurance, Springer, 2nd edition.

See Also

later.

Examples

## Crowder (1987b), Output Listing 1, trapezoidal quadrature (less accurate)

M <- 2
Ru <- 3
mu <- seq(0, 2, by=0.25)
LL <- LL2 <- rep(NA, length(mu))
for ( i in 1:length(mu) ) {
  LL[i] <- round( imr.arl(M, Ru, mu[i], 1), digits=4)
  LL2[i] <- round( imr.arl(M, Ru, mu[i], 1, cmode="Crowder", N=80), digits=4)
}
LL1987b <- c(18.2164, 16.3541, 12.4282, 8.7559, 6.1071, 4.3582, 3.2260, 2.4878, 1.9989)
print( data.frame(mu, LL2, LL1987b, LL) )

## Crowder (1987a), Table 1, trapezoidal quadrature (less accurate)

M <- 4
Ru <- 3
mu <- seq(0, 2, by=0.25)
LL <- rep(NA, length(mu))
for ( i in 1:length(mu) ) LL[i] <- round( imr.arl(M, Ru, mu[i], 1), digits=4)
LL1987a <- c(34.44, 34.28, 34.07, 33.81, 33.45, 32.82, 31.50, 28.85, 24.49)
print( data.frame(mu, LL1987a, LL) )

## Rigdon, Cruthis, Champ (1994), Table 1, Monte Carlo based

M <- 2.992
Ru <- 4.139
icARL <- imr.arl(M, Ru, 0, 1)
icARL1994 <- 200
print( data.frame(icARL1994, icARL) )

M <- 3.268
Ru <- 4.556
icARL <- imr.arl(M, Ru, 0, 1)
icARL1994 <- 500
print( data.frame(icARL1994, icARL) )

## ..., Table 2, Monte Carlo based

M <- 2.992
Ru <- 4.139
tau <- c(seq(1, 1.3, by=0.05), seq(1.4, 2, by=0.1))
LL <- rep(NA, length(tau))
for ( i in 1:length(tau) ) LL[i] <- round( imr.arl(M, Ru, 0, tau[i]), digits=2)
LL1994 <- c(200.54, 132.25, 90.84, 65.66, 49.35, 38.92, 31.11, 21.35, 15.47,
12.04, 9.81, 8.21, 7.03, 6.14)
print( data.frame(tau, LL1994, LL) )

## Radson, Alwan (1995), Table 2 (Monte Carlo based), half-normal, known parameter case
## two-sided (!) MR-alone (!) chart, hence the ARL results has to be decreased by 1
## Here: a large M (=12) is deployed to mimic Inf
alpha <- 0.00915
Ru <- sqrt(2) * qnorm(1-alpha/4)
Rl <- sqrt(2) * qnorm(0.5+alpha/4)
k <- 1.5 - (0:7)/10
LL <- rep(NA, length(k))
for ( i in 1:length(k) )
  LL[i] <- round( imr.arl(12, Ru, 0, k[i], vsided="two", Rl=Rl), digits=2) - 1
RA1995 <- c(18.61, 24.51, 34.21, 49.74, 75.08, 113.14, 150.15, 164.54)
print( data.frame(k, RA1995, LL) )

## Amin, Ethridge (1998), Table 2, column sigma/sigma_0 = 1.00

M <- 3.27
Ru <- 4.56
#M <- 3.268
#Ru <- 4.556
mu <- seq(0, 2, by=0.25)
LL <- rep(NA, length(mu))
for ( i in 1:length(mu) ) LL[i] <- round( imr.arl(M, Ru, mu[i], 1), digits=1)
LL1998 <- c(505.3, 427.6, 276.7, 156.2, 85.0, 46.9, 26.9, 16.1, 10.1)
print( data.frame(mu, LL1998, LL) )

## ..., column sigma/sigma_0 = 1.05

for ( i in 1:length(mu) ) LL[i] <- round( imr.arl(M, Ru, mu[i], 1.05), digits=1)
LL1998 <- c(296.8, 251.6, 169.6, 101.6, 58.9, 34.5, 20.9, 13.2, 8.7)
print( data.frame(mu, LL1998, LL) )

## Acosta-Mejia, Pignatiello (2000), Table 2
## AMP utilized Markov chain approximation
## However, the MR series is not Markovian!
## MR-alone (!) chart, hence the ARL results has to be decreased by 1
## Here: a large M (=8) is deployed to mimic Inf
Ru <- 3.93
sigma <- c(1, 1.05, 1.1, 1.15, 1.2, 1.3, 1.4, 1.5, 1.75)
LL <- rep(NA, length(sigma))
for ( i in 1:length(sigma) ) LL[i] <- round( imr.arl(8, Ru, 0, sigma[i], N=30), digits=1) - 1
AMP2000 <- c(201.0, 136.8, 97.9, 73.0, 56.3, 36.4, 25.6, 19.1, 11.0)
print( data.frame(sigma, AMP2000, LL) )

## Mark, Krehbiel (2011), Table 2, deployment of Crowder (1987b), nominal ic ARL 500

M <- c(3.09, 3.20, 3.30, 3.50, 4.00)
Ru <- c(6.00, 4.67, 4.53, 4.42, 4.36)
LL0 <- rep(NA, length(M))
for ( i in 1:length(M) ) LL0[i] <- round( imr.arl(M[i], Ru[i], 0, 1), digits=1)
print( data.frame(M, Ru, LL0) )

Compute ARLs of EWMA ln S^2 control charts (variance charts)

Description

Computation of the (zero-state) Average Run Length (ARL) for different types of EWMA control charts (based on the log of the sample variance S^2) monitoring normal variance.

Usage

lns2ewma.arl(l,cl,cu,sigma,df,hs=NULL,sided="upper",r=40)

Arguments

Details

lns2ewma.arl determines the Average Run Length (ARL) by numerically solving the related ARL integral equation by means of the Nystroem method based on Gauss-Legendre quadrature.

Value

Returns a single value which resembles the ARL.

Author(s)

Sven Knoth

References

S. V. Crowder and M. D. Hamilton (1992), An EWMA for monitoring a process standard deviation, Journal of Quality Technology 24, 12-21.

S. Knoth (2005), Accurate ARL computation for EWMA-S^2 control charts, Statistics and Computing 15, 341-352.

See Also

xewma.arl for zero-state ARL computation of EWMA control charts for monitoring normal mean.

Examples

lns2ewma.ARL <- Vectorize("lns2ewma.arl", "sigma")

## Crowder/Hamilton (1992)
## moments of ln S^2
E_log_gamma <- function(df) log(2/df) + digamma(df/2)
V_log_gamma <- function(df) trigamma(df/2)
E_log_gamma_approx <- function(df) -1/df - 1/3/df^2 + 2/15/df^4
V_log_gamma_approx <- function(df) 2/df + 2/df^2 + 4/3/df^3 - 16/15/df^5

## results from Table 3 ( upper chart with reflection at 0 = log(sigma0=1) )
## original entries are (lambda = 0.05, K = 1.06, df=n-1=4)
# sigma   ARL
# 1       200
# 1.1      43
# 1.2      18
# 1.3      11
# 1.4       7.6
# 1.5       6.0
# 2         3.2

df <- 4
lambda <- .05
K <- 1.06
cu <- K * sqrt( lambda/(2-lambda) * V_log_gamma_approx(df) )

sigmas <- c(1 + (0:5)/10, 2)
arls <- round(lns2ewma.ARL(lambda, 0, cu, sigmas, df, hs=0, sided="upper"), digits=1)
data.frame(sigmas, arls)

## Knoth (2005)
## compare with Table 3 (p. 351)
lambda <- .05
df <- 4
K <- 1.05521
cu <- 1.05521 * sqrt( lambda/(2-lambda) * V_log_gamma_approx(df) )

## upper chart with reflection at sigma0=1 in Table 4
## original entries are
# sigma   ARL_0    ARL_-.267
# 1       200.0    200.0
# 1.1      43.04    41.55
# 1.2      18.10    19.92
# 1.3      10.75    13.11
# 1.4       7.63     9.93
# 1.5       5.97     8.11
# 2         3.17     4.67

M <- -0.267
cuM <- lns2ewma.crit(lambda, 200, df, cl=M, hs=M, r=60)[2]
arls1 <- round(lns2ewma.ARL(lambda, 0, cu, sigmas, df, hs=0, sided="upper"), digits=2)
arls2 <- round(lns2ewma.ARL(lambda, M, cuM, sigmas, df, hs=M, sided="upper", r=60), digits=2)
data.frame(sigmas, arls1, arls2)

Compute critical values of EWMA ln S^2 control charts (variance charts)

Description

Computation of the critical values (similar to alarm limits) for different types of EWMA control charts (based on the log of the sample variance S^2) monitoring normal variance.

Usage

lns2ewma.crit(l,L0,df,sigma0=1,cl=NULL,cu=NULL,hs=NULL,sided="upper",mode="fixed",r=40)

Arguments

Details

lns2ewma.crit determines the critical values (similar to alarm limits) for given in-control ARL L0 by applying secant rule and using lns2ewma.arl(). In case of sided="two" and mode="unbiased" a two-dimensional secant rule is applied that also ensures that the maximum of the ARL function for given standard deviation is attained at sigma0. See Knoth (2010) and the related example.

Value

Returns the lower and upper control limit cl and cu.

Author(s)

Sven Knoth

References

C. A. Acosta-Mej\'ia and J. J. Pignatiello Jr. and B. V. Rao (1999), A comparison of control charting procedures for monitoring process dispersion, IIE Transactions 31, 569-579.

S. V. Crowder and M. D. Hamilton (1992), An EWMA for monitoring a process standard deviation, Journal of Quality Technology 24, 12-21.

S. Knoth (2005), Accurate ARL computation for EWMA-S^2 control charts, Statistics and Computing 15, 341-352.

S. Knoth (2010), Control Charting Normal Variance – Reflections, Curiosities, and Recommendations, in Frontiers in Statistical Quality Control 9, H.-J. Lenz and P.-T. Wilrich (Eds.), Physica Verlag, Heidelberg, Germany, 3-18.

See Also

lns2ewma.arl for calculation of ARL of EWMA ln S^2 control charts.

Examples

## Knoth (2005)
## compare with 1.05521 mentioned on page 350 third line from below
L0 <- 200
lambda <- .05
df <- 4
limits <- lns2ewma.crit(lambda, L0, df, cl=0, hs=0)
limits["cu"]/sqrt( lambda/(2-lambda)*(2/df+2/df^2+4/3/df^3-16/15/df^5) )

Compute ARLs of MEWMA control charts

Description

Computation of the (zero-state) Average Run Length (ARL) for multivariate exponentially weighted moving average (MEWMA) charts monitoring multivariate normal mean.

Usage

mewma.arl(l, cE, p, delta=0, hs=0, r=20, ntype=NULL, qm0=20, qm1=qm0)

mewma.arl.f(l, cE, p, delta=0, r=20, ntype=NULL, qm0=20, qm1=qm0)

mewma.ad(l, cE, p, delta=0, r=20, n=20, type="cond", hs=0, ntype=NULL, qm0=20, qm1=qm0)

Arguments

Details

Basically, this is the implementation of different numerical algorithms for solving the integral equation for the MEWMA in-control (delta = 0) ARL introduced in Rigdon (1995a) and out-of-control (delta != 0) ARL in Rigdon (1995b). Most of them are nothing else than the Nystroem approach – the integral is replaced by a suitable quadrature. Here, the Gauss-Legendre (more powerful), Radau (used by Rigdon, 1995a), Clenshaw-Curtis, and Simpson rule (which is really bad) are provided. Additionally, the collocation approach is offered as well, because it is much better for small odd values for p. FORTRAN code for the Radau quadrature based Nystroem of Rigdon (1995a) was published in Bodden and Rigdon (1999) – see also https://lib.stat.cmu.edu/jqt/31-1. Furthermore, FORTRAN code for the Markov chain approximation (in- and out-ot-control) could be found at https://lib.stat.cmu.edu/jqt/33-4/. The related papers are Runger and Prabhu (1996) and Molnau et al. (2001). The idea of the Clenshaw-Curtis quadrature was taken from Capizzi and Masarotto (2010), who successfully deployed a modified Clenshaw-Curtis quadrature to calculate the ARL of combined (univariate) Shewhart-EWMA charts. It turns out that it works also nicely for the MEWMA ARL. The version mewma.arl.f() without the argument hs provides the ARL as function of one (in-control) or two (out-of-control) arguments.

Value

Returns a single value which is simply the zero-state ARL.

Author(s)

Sven Knoth

References

Kevin M. Bodden and Steven E. Rigdon (1999), A program for approximating the in-control ARL for the MEWMA chart, Journal of Quality Technology 31(1), 120-123.

Giovanna Capizzi and Guido Masarotto (2010), Evaluation of the run-length distribution for a combined Shewhart-EWMA control chart, Statistics and Computing 20(1), 23-33.

Sven Knoth (2017), ARL Numerics for MEWMA Charts, Journal of Quality Technology 49(1), 78-89.

Wade E. Molnau et al. (2001), A Program for ARL Calculation for Multivariate EWMA Charts, Journal of Quality Technology 33(4), 515-521.

Sharad S. Prabhu and George C. Runger (1997), Designing a multivariate EWMA control chart, Journal of Quality Technology 29(1), 8-15.

Steven E. Rigdon (1995a), An integral equation for the in-control average run length of a multivariate exponentially weighted moving average control chart, J. Stat. Comput. Simulation 52(4), 351-365.

Steven E. Rigdon (1995b), A double-integral equation for the average run length of a multivariate exponentially weighted moving average control chart, Stat. Probab. Lett. 24(4), 365-373.

George C. Runger and Sharad S. Prabhu (1996), A Markov Chain Model for the Multivariate Exponentially Weighted Moving Averages Control Chart, J. Amer. Statist. Assoc. 91(436), 1701-1706.

See Also

mewma.crit for getting the alarm threshold to attain a certain in-control ARL.

Examples

# Rigdon (1995a), p. 357, Tab. 1
p <- 2
r <- 0.25
h4 <- c(8.37, 9.90, 11.89, 13.36, 14.82, 16.72)
for ( i in 1:length(h4) ) cat(paste(h4[i], "\t", round(mewma.arl(r, h4[i], p, ntype="ra")), "\n"))

r <- 0.1
h4 <- c(6.98, 8.63, 10.77, 12.37, 13.88, 15.88)
for ( i in 1:length(h4) ) cat(paste(h4[i], "\t", round(mewma.arl(r, h4[i], p, ntype="ra")), "\n"))


# Rigdon (1995b), p. 372, Tab. 1
## Not run: 
r <- 0.1
p <- 4
h <- 12.73
for ( sdelta in c(0, 0.125, 0.25, .5, 1, 2, 3) )
  cat(paste(sdelta, "\t",
      round(mewma.arl(r, h, p, delta=sdelta^2, ntype="ra", r=25), digits=2), "\n"))

p <- 5
h <- 14.56
for ( sdelta in c(0, 0.125, 0.25, .5, 1, 2, 3) )
  cat(paste(sdelta, "\t",
      round(mewma.arl(r, h, p, delta=sdelta^2, ntype="ra", r=25), digits=2), "\n"))

p <- 10
h <- 22.67
for ( sdelta in c(0, 0.125, 0.25, .5, 1, 2, 3) )
  cat(paste(sdelta, "\t",
      round(mewma.arl(r, h, p, delta=sdelta^2, ntype="ra", r=25), digits=2), "\n"))

## End(Not run)

# Runger/Prabhu (1996), p. 1704, Tab. 1
## Not run: 
r <- 0.1
p <- 4
H <- 12.73
cat(paste(0, "\t", round(mewma.arl(r, H, p, delta=0, ntype="mc", r=50), digits=2), "\n"))
for ( delta in c(.5, 1, 1.5, 2, 3) )
  cat(paste(delta, "\t",
      round(mewma.arl(r, H, p, delta=delta, ntype="mc", r=25), digits=2), "\n"))
# compare with Fortran program (MEWMA-ARLs.f90) from Molnau et al. (2001) with m1 = m2 = 25
# H4      P     R   DEL  ARL
# 12.73  4.  0.10  0.00 199.78
# 12.73  4.  0.10  0.50  35.05
# 12.73  4.  0.10  1.00  12.17
# 12.73  4.  0.10  1.50   7.22
# 12.73  4.  0.10  2.00   5.19
# 12.73  4.  0.10  3.00   3.42

p <- 20
H <- 37.01
cat(paste(0, "\t",
    round(mewma.arl(r, H, p, delta=0, ntype="mc", r=50), digits=2), "\n"))
for ( delta in c(.5, 1, 1.5, 2, 3) )
  cat(paste(delta, "\t",
      round(mewma.arl(r, H, p, delta=delta, ntype="mc", r=25), digits=2), "\n"))
# compare with Fortran program (MEWMA-ARLs.f90) from Molnau et al. (2001) with m1 = m2 = 25
# H4      P     R   DEL  ARL
# 37.01 20.  0.10  0.00 199.09
# 37.01 20.  0.10  0.50  61.62
# 37.01 20.  0.10  1.00  20.17
# 37.01 20.  0.10  1.50  11.40
# 37.01 20.  0.10  2.00   8.03
# 37.01 20.  0.10  3.00   5.18

## End(Not run)

# Knoth (2017), p. 85, Tab. 3, rows with p=3
## Not run: 
p <- 3
lambda <- 0.05
h4 <- mewma.crit(lambda, 200, p)
benchmark <- mewma.arl(lambda, h4, p, delta=1, r=50)
  
mc.arl  <- mewma.arl(lambda, h4, p, delta=1, r=25, ntype="mc")
ra.arl  <- mewma.arl(lambda, h4, p, delta=1, r=27, ntype="ra")
co.arl  <- mewma.arl(lambda, h4, p, delta=1, r=12, ntype="co2")
gl3.arl <- mewma.arl(lambda, h4, p, delta=1, r=30, ntype="gl3")
gl5.arl <- mewma.arl(lambda, h4, p, delta=1, r=25, ntype="gl5")
  
abs( benchmark - data.frame(mc.arl, ra.arl, co.arl, gl3.arl, gl5.arl) )

## End(Not run)

# Prabhu/Runger (1997), p. 13, Tab. 3
## Not run: 
p <- 2
r <- 0.1
H <- 8.64
cat(paste(0, "\t",
    round(mewma.ad(r, H, p, delta=0, type="cycl", ntype="mc", r=60), digits=2), "\n"))
for ( delta in c(.5, 1, 1.5, 2, 3) )
  cat(paste(delta, "\t",
      round(mewma.ad(r, H, p, delta=delta, type="cycl", ntype="mc", r=30), digits=2), "\n"))

# better accuracy
for ( delta in c(0, .5, 1, 1.5, 2, 3) )
  cat(paste(delta, "\t",
      round(mewma.ad(r, H, p, delta=delta^2, type="cycl", r=30), digits=2), "\n"))

## End(Not run)

Compute alarm threshold of MEWMA control charts

Description

Computation of the alarm threshold for multivariate exponentially weighted moving average (MEWMA) charts monitoring multivariate normal mean.

Usage

mewma.crit(l, L0, p, hs=0, r=20)

Arguments

Details

mewma.crit determines the alarm threshold of for given in-control ARL L0 by applying secant rule and using mewma.arl() with ntype="gl2".

Value

Returns a single value which resembles the critical value c.

Author(s)

Sven Knoth

References

Sven Knoth (2017), ARL Numerics for MEWMA Charts, Journal of Quality Technology 49(1), 78-89.

Steven E. Rigdon (1995), An integral equation for the in-control average run length of a multivariate exponentially weighted moving average control chart, J. Stat. Comput. Simulation 52(4), 351-365.

See Also

mewma.arl for zero-state ARL computation.

Examples

# Rigdon (1995), p. 358, Tab. 1
p <- 4
L0 <- 500
r <- .25
h4 <- mewma.crit(r, L0, p)
h4
## original value is 16.38.

# Knoth (2017), p. 82, Tab. 2
p <- 3
L0 <- 1e3
lambda <- c(0.25, 0.2, 0.15, 0.1, 0.05)
h4 <- rep(NA, length(lambda) )
for ( i in 1:length(lambda) ) h4[i] <- mewma.crit(lambda[i], L0, p, r=20)
round(h4, digits=2)
## original values are
## 15.82 15.62 15.31 14.76 13.60

Compute steady-state density of the MEWMA statistic

Description

Computation of the (zero-state) steady-state density function of the statistic deployed in multivariate exponentially weighted moving average (MEWMA) charts monitoring multivariate normal mean.

Usage

mewma.psi(l, cE, p, type="cond", hs=0, r=20)

Arguments

Details

Basically, ideas from Knoth (2017, MEWMA numerics) and Knoth (2016, steady-state ARL concepts) are merged. More details will follow.

Value

Returns a function.

Author(s)

Sven Knoth

References

Sven Knoth (2016), The Case Against the Use of Synthetic Control Charts, Journal of Quality Technology 48(2), 178-195.

Sven Knoth (2017), ARL Numerics for MEWMA Charts, Journal of Quality Technology 49(1), 78-89.

Sven Knoth (2018), The Steady-State Behavior of Multivariate Exponentially Weighted Moving Average Control Charts, Sequential Analysis 37(4), 511-529.

See Also

mewma.arl for calculating the in-control ARL of MEWMA.

Examples

lambda <- 0.1
L0 <- 200
p <- 3
h4 <- mewma.crit(lambda, L0, p)
x_ <- seq(0, h4*lambda/(2-lambda), by=0.002)
psi <- mewma.psi(lambda, h4, p)
psi_ <- psi(x_)
# plot(x_, psi_, type="l", xlab="x", ylab=expression(psi(x)), xlim=c(0,1.2))
# cf. to Figure 1 in Knoth (2018), p. 514, p=3

Compute ARLs of binomial EWMA p control charts

Description

Computation of the (zero-state) Average Run Length (ARL) at given rate p.

Usage

p.ewma.arl(lambda, ucl, n, p, z0, sided="upper", lcl=NULL, d.res=1,
r.mode="ieee.round", i.mode="integer")

Arguments

Details

The monitored data follow a binomial distribution with size n and failure/success probability p. The ARL values of the resulting EWMA control chart are determined by Markov chain approximation. Here, the original EWMA values are approximated by multiples of one over d.res. Different ways of rounding (see r.mode) to the next multiple are implemented. Besides Gan's paper nothing is published about the numerical subtleties.

Value

Return single value which resemble the ARL.

Author(s)

Sven Knoth

References

F. F. Gan (1990), Monitoring observations generated from a binomial distribution using modified exponentially weighted moving average control chart, J. Stat. Comput. Simulation 37, 45-60.

S. Knoth and S. Steinmetz (2013), EWMA p charts under sampling by variables, International Journal of Production Research 51, 3795-3807.

See Also

later.

Examples

## Gan (1990)

# Table 1

n <- 150
p0 <- .1
z0 <- n*p0

lambda <- c(1, .51, .165)
hu <- c(27, 22, 18)

p.value <- .1 + (0:20)/200

p.EWMA.arl <- Vectorize(p.ewma.arl, "p")

arl1.value <- round(p.EWMA.arl(lambda[1], hu[1], n, p.value, z0, r.mode="round"), digits=2)
arl2.value <- round(p.EWMA.arl(lambda[2], hu[2], n, p.value, z0, r.mode="round"), digits=2)
arl3.value <- round(p.EWMA.arl(lambda[3], hu[3], n, p.value, z0, r.mode="round"), digits=2)

arls <- matrix(c(arl1.value, arl2.value, arl3.value), ncol=length(lambda))
rownames(arls) <- p.value
colnames(arls) <- paste("lambda =", lambda)
arls

## Knoth/Steinmetz (2013)

n <- 5
p0 <- 0.02
z0 <- n*p0
lambda <- 0.3
ucl <- 0.649169922 ## in-control ARL 370.4 (determined with d.res = 2^14 = 16384)

res.list <- 2^(1:11)
arl.list <- NULL
for ( res in res.list ) {
  arl <- p.ewma.arl(lambda, ucl, n, p0, z0, d.res=res)
  arl.list <- c(arl.list, arl)
}
cbind(res.list, arl.list)

Compute ARLs of EWMA phat control charts

Description

Computation of the (zero-state) Average Run Length (ARL), upper control limit (ucl) for given in-control ARL, and lambda for minimal out-of control ARL at given shift.

Usage

phat.ewma.arl(lambda, ucl, mu, n, z0, sigma=1, type="known", LSL=-3, USL=3, N=15,
qm=25, ntype="coll")

phat.ewma.crit(lambda, L0, mu, n, z0, sigma=1, type="known", LSL=-3, USL=3, N=15, qm=25)

phat.ewma.lambda(L0, mu, n, z0, sigma=1, type="known", max_l=1, min_l=.001, LSL=-3, USL=3,
qm=25)

Arguments

Details

The three implemented functions allow to apply a new type control chart. Basically, lower and upper specification limits are given. The monitoring vehicle then is the empirical probability that an item will not follow these specification given the sequence of sample means. If the related EWMA sequence violates the control limits, then the alarm indicates a significant process deterioration. For details see the paper mentioned in the references. To be able to construct the control charts, see the first example.

Value

Return single values which resemble the ARL, the critical value, and the optimal lambda, respectively.

Author(s)

Sven Knoth

References

S. Knoth and S. Steinmetz (2013), EWMA p charts under sampling by variables, International Journal of Production Research 51, 3795-3807.

See Also

sewma.arl for a further collocation based ARL calculation routine.

Examples

## Simple example to demonstrate the chart.

# some functions
h.mu <- function(mu) pnorm(LSL-mu) + pnorm(mu-USL)
ewma <- function(x, lambda=0.1, z0=0) filter(lambda*x, 1-lambda, m="r", init=z0)

# parameters
LSL <- -3       # lower specification limit
USL <-  3	# upper specification limit
n <- 5		# batch size
lambda <- 0.1	# EWMA smoothing parameter
L0 <- 1000	# in-control Average Run Length (ARL)
z0 <- h.mu(0)	# start at minimal defect level
ucl <- phat.ewma.crit(lambda, L0, 0, n, z0, LSL=LSL, USL=USL)

# data
x0 <- matrix(rnorm(50*n), ncol=5)	# in-control data
x1 <- matrix(rnorm(50*n, mean=0.5), ncol=5)# out-of-control data
x <- rbind(x0,x1)			# all data

# create chart
xbar <- apply(x, 1, mean)
phat <- h.mu(xbar)
z <- ewma(phat, lambda=lambda, z0=z0)
plot(1:length(z), z, type="l", xlab="batch", ylim=c(0,.02))
abline(h=z0, col="grey", lwd=.7)
abline(h=ucl, col="red")


## S. Knoth, S. Steinmetz (2013)

# Table 1

lambdas <- c(.5, .25, .2, .1)
L0 <- 370.4
n <- 5
LSL <- -3
USL <- 3

phat.ewma.CRIT <- Vectorize("phat.ewma.crit", "lambda")
p.star <- pnorm( LSL ) + pnorm( -USL ) ## lower bound of the chart
ucls <- phat.ewma.CRIT(lambdas, L0, 0, n, p.star, LSL=LSL, USL=USL)
print(cbind(lambdas, ucls))

# Table 2

mus <- c((0:4)/4, 1.5, 2, 3)
phat.ewma.ARL <- Vectorize("phat.ewma.arl", "mu")
arls <- NULL
for ( i in 1:length(lambdas) ) {
  arls <- cbind(arls, round(phat.ewma.ARL(lambdas[i], ucls[i], mus,
                n, p.star, LSL=LSL, USL=USL), digits=2))
}
arls <- data.frame(arls, row.names=NULL)
names(arls) <- lambdas
print(arls)

# Table 3

## Not run: 
mus <- c(.25, .5, 1, 2)
phat.ewma.LAMBDA <- Vectorize("phat.ewma.lambda", "mu")
lambdas <- phat.ewma.LAMBDA(L0, mus, n, p.star, LSL=LSL, USL=USL)
print(cbind(mus, lambdas))
## End(Not run)

Compute ARLs of Poisson CUSUM control charts

Description

Computation of the (zero-state) Average Run Length (ARL) at given mean mu.

Usage

pois.cusum.arl(mu, km, hm, m, i0=0, sided="upper", rando=FALSE,
gamma=0, km2=0, hm2=0, m2=0, i02=0, gamma2=0)

Arguments

Details

The monitored data follow a Poisson distribution with mu. The ARL values of the resulting EWMA control chart are determined via Markov chain calculations. We follow the algorithm given in Lucas (1985) expanded with some arithmetic 'tricks' (e.g., by deploying Toeplitz matrix algebra). A paper explaining it is under preparation.

Value

Returns a single value which resembles the ARL.

Author(s)

Sven Knoth

References

J. M. Lucas (1985) Counted data CUSUM's, Technometrics 27(2), 129-144.

C. H. White and J. B. Keats (1996) ARLs and Higher-Order Run-Length Moments for the Poisson CUSUM, Journal of Quality Technology 28(3), 363-369.

C. H. White, J. B. Keats and J. Stanley (1997) Poisson CUSUM versus c chart for defect data, Quality Engineering 9(4), 673-679.

G. Rossi and L. Lampugnani and M. Marchi (1999), An approximate CUSUM procedure for surveillance of health events, Statistics in Medicine 18(16), 2111-2122.

S. W. Han, K.-L. Tsui, B. Ariyajunya, and S. B. Kim (2010), A comparison of CUSUM, EWMA, and temporal scan statistics for detection of increases in poisson rates, Quality and Reliability Engineering International 26(3), 279-289.

M. B. Perry and J. J. Pignatiello Jr. (2011) Estimating the time of step change with Poisson CUSUM and EWMA control charts, International Journal of Production Research 49(10), 2857-2871.

See Also

later.

Examples

## Lucas 1985, upper chart (Tables 2 and 3)
k   <- .25
h   <- 10
m   <- 4
km  <- m * k
hm  <- m * h
mu0 <- 1 * k
ARL <- pois.cusum.arl(mu0, km, hm-1, m)
# Lucas reported 438 (in Table 2, first block, row 10.0 .25 .0 ..., column 1.0
# Recall that Lucas and other trigger an alarm, if the CUSUM statistic is greater than
# or equal to the alarm threshold h
print(ARL)

ARL <- pois.cusum.arl(mu0, km, hm-1, m, i0=round((hm-1)/2))
# Lucas reported 333 (in Table 3, first block, row 10.0 .25 .0 ..., column 1.0
print(ARL)

## Lucas 1985, lower chart (Tables 4 and 5)
ARL <- pois.cusum.arl(mu0, km, hm-1, m, sided="lower")
# Lucas reported 437 (in Table 4, first block, row 10.0 .25 .0 ..., column 1.0
print(ARL)

ARL <- pois.cusum.arl(mu0, km, hm-1, m, i0=round((hm-1)/2), sided="lower")
# Lucas reported 318 (in Table 5, first block, row 10.0 .25 .0 ..., column 1.0
print(ARL)

Compute alarm thresholds and randomization constants of Poisson CUSUM control charts

Description

Computation of the CUSUM upper limit and, if needed, of the randomization probability, given mean mu0.

Usage

pois.cusum.crit(mu0, km, A, m, i0=0, sided="upper", rando=FALSE)

Arguments

Details

The monitored data follow a Poisson distribution with mu (here the in-control level mu0). The ARL values of the resulting EWMA control chart are determined via Markov chain calculations. With some grid search, we obtain the smallest value for the integer threshold component hm so that the resulting ARL is not smaller than A. If equality is needed, then activating rando=TRUE yields the corresponding randomization probability gamma. More details will follow in a paper that will be submitted in 2020.

Value

Returns two single values, integer threshold hm resulting in the final alarm threshold h=hm/m, and the randomization probability.

Author(s)

Sven Knoth

References

J. M. Lucas (1985) Counted data CUSUM's, Technometrics 27(2), 129-144.

C. H. White and J. B. Keats (1996) ARLs and Higher-Order Run-Length Moments for the Poisson CUSUM, Journal of Quality Technology 28(3), 363-369.

C. H. White, J. B. Keats and J. Stanley (1997) Poisson CUSUM versus c chart for defect data, Quality Engineering 9(4), 673-679.

G. Rossi and L. Lampugnani and M. Marchi (1999), An approximate CUSUM procedure for surveillance of health events, Statistics in Medicine 18(16), 2111-2122.

S. W. Han, K.-L. Tsui, B. Ariyajunya, and S. B. Kim (2010), A comparison of CUSUM, EWMA, and temporal scan statistics for detection of increases in poisson rates, Quality and Reliability Engineering International 26(3), 279-289.

M. B. Perry and J. J. Pignatiello Jr. (2011) Estimating the time of step change with Poisson CUSUM and EWMA control charts, International Journal of Production Research 49(10), 2857-2871.

See Also

later.

Examples

## Lucas 1985
mu0 <- 0.25
km <- 1
A <- 430
m  <- 4
#cv <- pois.cusum.crit(mu0, km, A, m)
cv <- c(40, 0)
# Lucas reported h = 10 alias hm = 40 (in Table 2, first block, row 10.0 .25 .0 ..., column 1.0
# Recall that Lucas and other trigger an alarm, if the CUSUM statistic is greater than
# or equal to the alarm threshold h
print(cv)

Compute the CUSUM k and h for given in-control ARL L0 and out-of-control ARL L1, Poisson case

Description

Computation of the reference value k and the alarm threshold h for one-sided CUSUM control charts monitoring Poisson data, if the in-control ARL L0 and the out-of-control ARL L1 are given.

Usage

pois.cusum.crit.L0L1(mu0, L0, L1, sided="upper", OUTPUT=FALSE)

Arguments

Details

pois.cusum.crit.L0L1 determines the reference value k and the alarm threshold h for given in-control ARL L0 and out-of-control ARL L1 by applying grid search and using pois.cusum.arl() and pois.cusum.crit(). These CUSUM design rules were firstly (and quite rarely afterwards) used by Ewan and Kemp. In the Poisson case, Rossi et al. applied them while analyzing three different normal approximations of the Poisson distribution. See the example which illustrates the validity of all these approaches.

Value

Returns a data frame with results for the denominator m of the rational approximation, km as (integer) enumerator of the reference value (approximation), the corresponding out-of-control mean mu1, the final approximation k of the reference value, the threshold values hm (integer) and h (=hm/m), and the randomization constant gamma (the target in-control ARL is exactly matched).

Author(s)

Sven Knoth

References

W. D. Ewan and K. W. Kemp (1960), Sampling inspection of continuous processes with no autocorrelation between successive results, Biometrika 47 (3/4), 363-380.

K. W. Kemp (1962), The Use of Cumulative Sums for Sampling Inspection Schemes, Journal of the Royal Statistical Sociecty C, Applied Statistics 11(1), 16-31.

G. Rossi, L. Lampugnani and M. Marchi (1999), An approximate CUSUM procedure for surveillance of health events, Statistics in Medicine 18(16), 2111-2122.

See Also

pois.cusum.arl for zero-state ARL and pois.cusum.crit for threshold h computation.

Examples

## Table 1 from Rossi et al. (1999) -- one-sided CUSUM
La <- 500 # in-control ARL
Lr <- 7 # out-of-control ARL
m_a <- 0.52 # in-control mean of the Poisson variate
## Not run: kh <- xcusum.crit.L0L1(La, Lr, sided="one")
# kh <- ...: instead of deploying EK1960, one could use more accurate numbers
EK_k <- 0.60 # EK1960 results in
EK_h <- 3.80 # Table 2 on p. 372
eZR <- 2*EK_h # reproduce normal ooc mean from reference value k
m_r <- 1.58 # EK1960 Table 3 on p. 377 for m_a = 0.52
R1 <- round( eZR/sqrt(m_a) + 1, digits=2)
R2 <- round( ( eZR/2/sqrt(m_a) + 1 )^2, digits=2)
R3 <- round(( sqrt(4 + 2*eZR/sqrt(m_a)) - 1 )^2, digits=2)
RS <- round( m_r / m_a, digits=2 )
## Not run: K_hk <- pois.cusum.crit.L0L1(m_a, La, Lr) # 'our' 'exact' approach
K_hk <- data.frame(m=1000, km=948, mu1=1.563777, k=0.948, hm=3832, h=3.832, gamma=0.1201901)
# get k for competing means mu0 (m_a) and mu1 (m_r)
k_m01 <- function(mu0, mu1) (mu1 - mu0) / (log(mu1) - log(mu0))
# get ooc mean mu1 (m_r) for given mu0 (m_a) and reference value k
m1_km0 <- function(mu0, k) {
  zero <- function(x) k - k_m01(mu0,x)
  upper <- mu0 + .5
  while ( zero(upper) > 0 ) upper <- upper + 0.5
  mu1 <- uniroot(zero, c(mu0*1.00000001, upper), tol=1e-9)$root
  mu1
}
K_m_r <- m1_km0(m_a, K_hk$k)
RK <- round( K_m_r / m_a, digits=2 )
cat(paste(m_a, R1, R2, R3, RS, RK, "\n", sep="\t"))

Compute steady-state ARLs of Poisson EWMA control charts

Description

Computation of the steady-state Average Run Length (ARL) at given mean mu.

Usage

pois.ewma.ad(lambda, AL, AU, mu0, mu, sided="two", rando=FALSE, gL=0, gU=0,
mcdesign="classic", N=101)

Arguments

Details

The monitored data follow a Poisson distribution with mu. The ARL values of the resulting EWMA control chart are determined by Markov chain approximation. We follow the algorithm given in Borror, Champ and Rigdon (1998). The function is in an early development phase.

Value

Return single value which resembles the steady-state ARL.

Author(s)

Sven Knoth

References

C. M. Borror, C. W. Champ and S. E. Rigdon (1998) Poisson EWMA control charts, Journal of Quality Technonlogy 30(4), 352-361.

M. C. Morais and S. Knoth (2020) Improving the ARL profile and the accuracy of its calculation for Poisson EWMA charts, Quality and Reliability Engineering International 36(3), 876-889.

See Also

later.

Examples

## Borror, Champ and Rigdon (1998), Table 2, PEWMA column
mu0 <- 20
lambda <- 0.27
A <- 3.319
mu1  <- c(2*(3:15), 35)
ARL1 <- AD1 <- rep(NA, length(mu1))
for ( i in 1:length(mu1) ) {
  ARL1[i] <- round(pois.ewma.arl(lambda,A,A,mu0,mu0,mu1[i],mcdesign="classic"),digits=1)
  AD1[i]  <- round(pois.ewma.ad(lambda,A,A,mu0,mu1[i],mcdesign="classic"),digits=1)
}
print( cbind(mu1, ARL1, AD1) )

## Morais and Knoth (2020), Table 2, lambda = 0.27 column
## randomisation not implemented for pois.ewma.ad()
lambda <- 0.27
AL <- 3.0870
AU <- 3.4870
gL <- 0.001029
gU <- 0.000765
mu2  <- c(16, 18, 19.99, mu0, 20.01, 22, 24)
ARL2 <- AD2 <- rep(NA, length(mu2))
for ( i in 1:length(mu2) ) {
  ARL2[i] <- round(pois.ewma.arl(lambda,AL,AU,mu0,mu0,mu2[i],rando=FALSE), digits=1)
  AD2[i] <- round(pois.ewma.ad(lambda,AL,AU,mu0,mu2[i],rando=FALSE), digits=1)
}
print( cbind(mu2, ARL2, AD2) )

Compute ARLs of Poisson EWMA control charts

Description

Computation of the (zero-state) Average Run Length (ARL) at given mean mu.

Usage

pois.ewma.arl(lambda, AL, AU, mu0, z0, mu, sided="two", rando=FALSE, gL=0, gU=0,
mcdesign="transfer", N=101)

Arguments

Details

The monitored data follow a Poisson distribution with mu. The ARL values of the resulting EWMA control chart are determined by Markov chain approximation. We follow the algorithm given in Borror, Champ and Rigdon (1998). However, by setting mcdesign="transfer" (now the default) from Morais and Knoth (2020), the accuracy is considerably improved.

Value

Return single value which resembles the ARL.

Author(s)

Sven Knoth

References

C. M. Borror, C. W. Champ and S. E. Rigdon (1998) Poisson EWMA control charts, Journal of Quality Technonlogy 30(4), 352-361.

M. C. Morais and S. Knoth (2020) Improving the ARL profile and the accuracy of its calculation for Poisson EWMA charts, Quality and Reliability Engineering International 36(3), 876-889.

See Also

later.

Examples

## Borror, Champ and Rigdon (1998), Table 2, PEWMA column
mu0 <- 20
lambda <- 0.27
A <- 3.319
mu1  <- c(2*(3:15), 35)
ARL1 <- rep(NA, length(mu1))
for ( i in 1:length(mu1) )
  ARL1[i] <- pois.ewma.arl(lambda, A, A, mu0, mu0, mu1[i], mcdesign="classic")
print(cbind(mu1, round(ARL1, digits=1)))

## the same numbers with improved accuracy
ARL2 <- rep(NA, length(mu1))
for ( i in 1:length(mu1) )
  ARL2[i] <- pois.ewma.arl(lambda, A, A, mu0, mu0, mu1[i], mcdesign="transfer")
print(cbind(mu1, round(ARL2, digits=1)))

## Morais and Knoth (2020), Table 2, lambda = 0.27 column
lambda <- 0.27
AL <- 3.0870
AU <- 3.4870
gL <- 0.001029
gU <- 0.000765
mu0 <- 20
mu1  <- c(16, 18, 19.99, mu0, 20.01, 22, 24)
ARL3 <- rep(NA, length(mu1))
for ( i in 1:length(mu1) )
  ARL3[i] <- pois.ewma.arl(lambda,AL,AU,mu0,mu0,mu1[i],rando=TRUE,gL=gL,gU=gU, N=101)
print(cbind(mu1, round(ARL3, digits=1)))

Compute ARLs of Poisson EWMA control charts

Description

Computation of the (zero-state) Average Run Length (ARL) at given mean mu.

Usage

pois.ewma.crit(lambda, L0, mu0, z0, AU=3, sided="two", design="sym", rando=FALSE,
mcdesign="transfer", N=101, jmax=4)

Arguments

Details

The monitored data follow a Poisson distribution with mu. Here we solve the inverse task to the usual ARL calculation. Hence, determine the control limit factors so that the in-control ARL is (roughly) equal to L0. The ARL values underneath the routine are determined by Markov chain approximation. The algorithm is just a grid search that takes care of the discrete ARL behavior.

Value

Return one or two values being he control limit factors.

Author(s)

Sven Knoth

References

C. M. Borror, C. W. Champ and S. E. Rigdon (1998) Poisson EWMA control charts, Journal of Quality Technonlogy 30(4), 352-361.

M. C. Morais and S. Knoth (2020) Improving the ARL profile and the accuracy of its calculation for Poisson EWMA charts, Quality and Reliability Engineering International 36(3), 876-889.

See Also

later.

Examples

## Borror, Champ and Rigdon (1998), page 30, original value is A = 2.8275
mu0 <- 4
lambda <- 0.2
L0 <- 351
A <- pois.ewma.crit(lambda, L0, mu0, mu0, mcdesign="classic")
print(round(A, digits=4))

## Morais and Knoth (2020), Table 2, lambda = 0.27 column
lambda <- 0.27
L0 <- 1233.4
ccgg <- pois.ewma.crit(lambda,1233.4,mu0,mu0,design="unb",rando=TRUE,mcdesign="transfer")
print(ccgg, digits=3)

Calculate quadrature nodes and weights

Description

Computation of the nodes and weights to enable numerical quadrature.

Usage

quadrature.nodes.weights(n, type="GL", x1=-1, x2=1)

Arguments

Details

A more detailed description will follow soon. The algorithm for the Gauss-Legendre quadrature was delivered by Knut Petras to me, while the one for the Radau quadrature was taken from John Burkardt.

Value

Returns two vectors which hold the needed quadrature nodes and weights.

Author(s)

Sven Knoth

References

H. Brass and K. Petras (2011), Quadrature Theory. The Theory of Numerical Integration on a Compact Interval, Mathematical Surveys and Monographs, American Mathematical Society.

John Burkardt (2015), https://people.math.sc.edu/Burkardt/c_src/quadrule/quadrule.c

See Also

Many of the ARL routines use the Gauss-Legendre nodes.

Examples

# GL
n <- 10
qnw <-quadrature.nodes.weights(n, type="GL")
qnw

# Radau
n <- 10
qnw <-quadrature.nodes.weights(n, type="Ra")
qnw

Compute ARLs of CUSUM control charts (variance charts)

Description

Computation of the (zero-state) Average Run Length (ARL) for different types of CUSUM control charts (based on the sample variance S^2) monitoring normal variance.

Usage

scusum.arl(k, h, sigma, df, hs=0, sided="upper", k2=NULL,
h2=NULL, hs2=0, r=40, qm=30, version=2)

Arguments

Details

scusum.arl determines the Average Run Length (ARL) by numerically solving the related ARL integral equation by means of collocation (piecewise Chebyshev polynomials).

Value

Returns a single value which resembles the ARL.

Author(s)

Sven Knoth

References

S. Knoth (2005), Accurate ARL computation for EWMA-S^2 control charts, Statistics and Computing 15, 341-352.

S. Knoth (2006), Computation of the ARL for CUSUM-S^2 schemes, Computational Statistics & Data Analysis 51, 499-512.

See Also

xcusum.arl for zero-state ARL computation of CUSUM control charts for monitoring normal mean.

Examples

## Knoth (2006)
## compare with Table 1 (p. 507)
k <- 1.46 # sigma1 = 1.5
df <- 1
h <- 10

# original values
# sigma coll63       BE     Hawkins  MC 10^9 (s.e.)
# 1     260.7369  260.7546  261.32  260.7399 (0.0081)
# 1.1    90.1319   90.1389   90.31   90.1319 (0.0027)
# 1.2    43.6867   43.6897   43.75   43.6845 (0.0013)
# 1.3    26.2916   26.2932   26.32   26.2929 (0.0007)
# 1.4    18.1231   18.1239   18.14   18.1235 (0.0005)
# 1.5    13.6268   13.6273   13.64   13.6272 (0.0003)
# 2       5.9904    5.9910    5.99    5.9903 (0.0001)
# replicate the column coll63
sigma <- c(1, 1.1, 1.2, 1.3, 1.4, 1.5, 2)
arl <- rep(NA, length(sigma))
for ( i in 1:length(sigma) )
  arl[i] <- round(scusum.arl(k, h, sigma[i], df, r=63, qm=20, version=2), digits=4)
data.frame(sigma, arl)

Compute decision intervals of CUSUM control charts (variance charts)

Description

omputation of the decision intervals (alarm limits) for different types of CUSUM control charts (based on the sample variance S^2) monitoring normal variance.

Usage

scusum.crit(k, L0, sigma, df, hs=0, sided="upper", mode="eq.tails",
k2=NULL, hs2=0, r=40, qm=30)

Arguments

Details

scusum.crit ddetermines the decision interval (alarm limit) for given in-control ARL L0 by applying secant rule and using scusum.arl().

Value

Returns a single value which resembles the decision interval h.

Author(s)

Sven Knoth

References

S. Knoth (2005), Accurate ARL computation for EWMA-S^2 control charts, Statistics and Computing 15, 341-352.

S. Knoth (2006), Computation of the ARL for CUSUM-S^2 schemes, Computational Statistics & Data Analysis 51, 499-512.

See Also

xcusum.arl for zero-state ARL computation of CUSUM control charts monitoring normal mean.

Examples

## Knoth (2006)
## compare with Table 1 (p. 507)
k <- 1.46 # sigma1 = 1.5
df <- 1
L0 <- 260.74
h <- scusum.crit(k, L0, 1, df)
h
# original value is 10

Compute ARLs of CUSUM-Shewhart control charts (variance charts)

Description

Computation of the (zero-state) Average Run Length (ARL) for different types of CUSUM-Shewhart combo control charts (based on the sample variance S^2) monitoring normal variance.

Usage

scusums.arl(k, h, cS, sigma, df, hs=0, sided="upper", k2=NULL,
h2=NULL, hs2=0, r=40, qm=30, version=2)

Arguments

Details

scusums.arl determines the Average Run Length (ARL) by numerically solving the related ARL integral equation by means of collocation (piecewise Chebyshev polynomials).

Value

Returns a single value which resembles the ARL.

Author(s)

Sven Knoth

References

S. Knoth (2006), Computation of the ARL for CUSUM-S^2 schemes, Computational Statistics & Data Analysis 51, 499-512.

See Also

scusum.arl for zero-state ARL computation of standalone CUSUM control charts for monitoring normal variance.

Examples

## will follow

Compute ARLs of EWMA control charts (variance charts)

Description

Computation of the (zero-state) Average Run Length (ARL) for different types of EWMA control charts (based on the sample variance S^2) monitoring normal variance.

Usage

sewma.arl(l,cl,cu,sigma,df,s2.on=TRUE,hs=NULL,sided="upper",r=40,qm=30)

Arguments

Details

sewma.arl determines the Average Run Length (ARL) by numerically solving the related ARL integral equation by means of collocation (Chebyshev polynomials).

Value

Returns a single value which resembles the ARL.

Author(s)

Sven Knoth

References

S. Knoth (2005), Accurate ARL computation for EWMA-S^2 control charts, Statistics and Computing 15, 341-352.

S. Knoth (2006), Computation of the ARL for CUSUM-S^2 schemes, Computational Statistics & Data Analysis 51, 499-512.

See Also

xewma.arl for zero-state ARL computation of EWMA control charts for monitoring normal mean.

Examples

## Knoth (2005)
## compare with Table 1 (p. 347): 249.9997
## Monte Carlo with 10^9 replicates: 249.9892 +/- 0.008
l <- .025
df <- 1
cu <- 1 + 1.661865*sqrt(l/(2-l))*sqrt(2/df)
sewma.arl(l,0,cu,1,df)

## ARL values for upper and lower EWMA charts with reflecting barriers
## (reflection at in-control level sigma0 = 1)
## examples from Knoth (2006), Tables 4 and 5

Ssewma.arl <- Vectorize("sewma.arl", "sigma")

## upper chart with reflection at sigma0=1 in Table 4
## original entries are
# sigma   ARL
# 1       100.0
# 1.01    85.3
# 1.02    73.4
# 1.03    63.5
# 1.04    55.4
# 1.05    48.7
# 1.1     27.9
# 1.2     12.9
# 1.3     7.86
# 1.4     5.57
# 1.5     4.30
# 2       2.11

## Not run: 
l <- 0.15
df <- 4
cu <- 1 + 2.4831*sqrt(l/(2-l))*sqrt(2/df)
sigmas <- c(1 + (0:5)/100, 1 + (1:5)/10, 2)
arls <- round(Ssewma.arl(l, 1, cu, sigmas, df, sided="Rupper", r=100), digits=2)
data.frame(sigmas, arls)
## End(Not run)

## lower chart with reflection at sigma0=1 in Table 5
## original entries are
# sigma   ARL
# 1       200.04
# 0.9     38.47
# 0.8     14.63
# 0.7     8.65
# 0.6     6.31

## Not run: 
l <- 0.115
df <- 5
cl <- 1 - 2.0613*sqrt(l/(2-l))*sqrt(2/df)
sigmas <- c((10:6)/10)
arls <- round(Ssewma.arl(l, cl, 1, sigmas, df, sided="Rlower", r=100), digits=2)
data.frame(sigmas, arls)
## End(Not run)

Compute ARLs of EWMA control charts (variance charts) in case of estimated parameters

Description

Computation of the (zero-state) Average Run Length (ARL) for EWMA control charts (based on the sample variance S^2) monitoring normal variance with estimated parameters.

Usage

sewma.arl.prerun(l, cl, cu, sigma, df1, df2, hs=1, sided="upper",
r=40, qm=30, qm.sigma=30, truncate=1e-10)

Arguments

Details

Essentially, the ARL function sewma.arl is convoluted with the distribution of the sample standard deviation. For details see Jones/Champ/Rigdon (2001) and Knoth (2014?).

Value

Returns a single value which resembles the ARL.

Author(s)

Sven Knoth

References

L. A. Jones, C. W. Champ, S. E. Rigdon (2001), The performance of exponentially weighted moving average charts with estimated parameters, Technometrics 43, 156-167.

S. Knoth (2005), Accurate ARL computation for EWMA-S^2 control charts, Statistics and Computing 15, 341-352.

S. Knoth (2006), Computation of the ARL for CUSUM-S^2 schemes, Computational Statistics & Data Analysis 51, 499-512.

See Also

sewma.arl for zero-state ARL function of EWMA control charts w/o pre run uncertainty.

Examples

## will follow

Compute critical values of EWMA control charts (variance charts)

Description

Computation of the critical values (similar to alarm limits) for different types of EWMA control charts (based on the sample variance S^2) monitoring normal variance.

Usage

sewma.crit(l,L0,df,sigma0=1,cl=NULL,cu=NULL,hs=NULL,s2.on=TRUE,
sided="upper",mode="fixed",ur=4,r=40,qm=30)

Arguments

Details

sewma.crit determines the critical values (similar to alarm limits) for given in-control ARL L0 by applying secant rule and using sewma.arl(). In case of sided="two" and mode="unbiased" a two-dimensional secant rule is applied that also ensures that the maximum of the ARL function for given standard deviation is attained at sigma0. See Knoth (2010) and the related example.

Value

Returns the lower and upper control limit cl and cu.

Author(s)

Sven Knoth

References

H.-J. Mittag and D. Stemann and B. Tewes (1998), EWMA-Karten zur \"Uberwachung der Streuung von Qualit\"atsmerkmalen, Allgemeines Statistisches Archiv 82, 327-338,

C. A. Acosta-Mej\'ia and J. J. Pignatiello Jr. and B. V. Rao (1999), A comparison of control charting procedures for monitoring process dispersion, IIE Transactions 31, 569-579.

S. Knoth (2005), Accurate ARL computation for EWMA-S^2 control charts, Statistics and Computing 15, 341-352.

S. Knoth (2006a), Computation of the ARL for CUSUM-S^2 schemes, Computational Statistics & Data Analysis 51, 499-512.

S. Knoth (2006b), The art of evaluating monitoring schemes – how to measure the performance of control charts? in Frontiers in Statistical Quality Control 8, H.-J. Lenz and P.-T. Wilrich (Eds.), Physica Verlag, Heidelberg, Germany, 74-99.

S. Knoth (2010), Control Charting Normal Variance – Reflections, Curiosities, and Recommendations, in Frontiers in Statistical Quality Control 9, H.-J. Lenz and P.-T. Wilrich (Eds.), Physica Verlag, Heidelberg, Germany, 3-18.

See Also

sewma.arl for calculation of ARL of variance charts.

Examples

## Mittag et al. (1998)
## compare their upper critical value 2.91 that
## leads to the upper control limit via the formula shown below
## (for the usual upper EWMA \eqn{S^2}{S^2}).
## See Knoth (2006b) for a discussion of this EWMA setup and it's evaluation.

l  <- 0.18
L0 <- 250
df <- 4
limits <- sewma.crit(l, L0, df)
limits["cu"]

limits.cu.mittag_et_al <- 1 + sqrt(l/(2-l))*sqrt(2/df)*2.91
limits.cu.mittag_et_al

## Knoth (2005)
## reproduce the critical value given in Figure 2 (c=1.661865) for
## upper EWMA \eqn{S^2}{S^2} with df=1

l  <- 0.025
L0 <- 250
df <- 1
limits <- sewma.crit(l, L0, df)
cv.Fig2 <- (limits["cu"]-1)/( sqrt(l/(2-l))*sqrt(2/df) )
cv.Fig2

## the small difference (sixth digit after decimal point) stems from
## tighter criterion in the secant rule implemented in the R package.

## demo of unbiased ARL curves
## Deploy, please, not matrix dimensions smaller than 50 -- for the
## sake of accuracy, the value 80 was used.
## Additionally, this example needs between 1 and 2 minutes on a 1.6 Ghz box. 

## Not run: 
l  <- 0.1
L0 <- 500
df <- 4
limits <- sewma.crit(l, L0, df, sided="two", mode="unbiased", r=80)
SEWMA.arl <- Vectorize(sewma.arl, "sigma")
SEWMA.ARL <- function(sigma) 
  SEWMA.arl(l, limits[1], limits[2], sigma, df, sided="two", r=80)
layout(matrix(1:2, nrow=1))
curve(SEWMA.ARL, .75, 1.25, log="y")
curve(SEWMA.ARL, .95, 1.05, log="y")
## End(Not run)
# the above stuff needs about 1 minute

## control limits for upper and lower EWMA charts with reflecting barriers
## (reflection at in-control level sigma0 = 1)
## examples from Knoth (2006a), Tables 4 and 5

## Not run: 
## upper chart with reflection at sigma0=1 in Table 4: c = 2.4831
l <- 0.15
L0 <- 100
df <- 4
limits <- sewma.crit(l, L0, df, cl=1, sided="Rupper", r=100)
cv.Tab4 <- (limits["cu"]-1)/( sqrt(l/(2-l))*sqrt(2/df) )
cv.Tab4

## lower chart with reflection at sigma0=1 in Table 5: c = 2.0613
l <- 0.115
L0 <- 200
df <- 5
limits <- sewma.crit(l, L0, df, cu=1, sided="Rlower", r=100)
cv.Tab5 <- -(limits["cl"]-1)/( sqrt(l/(2-l))*sqrt(2/df) )
cv.Tab5
## End(Not run)

Compute critical values of of EWMA (variance charts) control charts under pre-run uncertainty

Description

Computation of quantiles of the Run Length (RL) for EWMA control charts monitoring normal variance.

Usage

sewma.crit.prerun(l,L0,df1,df2,sigma0=1,cl=NULL,cu=NULL,hs=1,sided="upper",
mode="fixed",r=40,qm=30,qm.sigma=30,truncate=1e-10,
tail_approx=TRUE,c.error=1e-10,a.error=1e-9)

Arguments

Details

sewma.crit.prerun determines the critical values (similar to alarm limits) for given in-control ARL L0 by applying secant rule and using sewma.arl.prerun(). In case of sided="two" and mode="unbiased" a two-dimensional secant rule is applied that also ensures that the maximum of the ARL function for given standard deviation is attained at sigma0. See Knoth (2010) for some details of the algorithm involved.

Value

Returns the lower and upper control limit cl and cu.

Author(s)

Sven Knoth

References

H.-J. Mittag and D. Stemann and B. Tewes (1998), EWMA-Karten zur \"Uberwachung der Streuung von Qualit\"atsmerkmalen, Allgemeines Statistisches Archiv 82, 327-338, S. Knoth (2005), Accurate ARL computation for EWMA-S^2 control charts, Statistics and Computing 15, 341-352.

S. Knoth (2010), Control Charting Normal Variance – Reflections, Curiosities, and Recommendations, in Frontiers in Statistical Quality Control 9, H.-J. Lenz and P.-T. Wilrich (Eds.), Physica Verlag, Heidelberg, Germany, 3-18.

See Also

sewma.arl.prerun for calculation of ARL of variance charts under pre-run uncertainty and sewma.crit for the algorithm w/o pre-run uncertainty.

Examples

## will follow

Compute RL quantiles of EWMA (variance charts) control charts

Description

Computation of quantiles of the Run Length (RL) for EWMA control charts monitoring normal variance.

Usage

sewma.q(l, cl, cu, sigma, df, alpha, hs=1, sided="upper", r=40, qm=30)

sewma.q.crit(l,L0,alpha,df,sigma0=1,cl=NULL,cu=NULL,hs=1,sided="upper",
mode="fixed",ur=4,r=40,qm=30,c.error=1e-12,a.error=1e-9)

Arguments

Details

Instead of the popular ARL (Average Run Length) quantiles of the EWMA stopping time (Run Length) are determined. The algorithm is based on Waldmann's survival function iteration procedure. Thereby the ideas presented in Knoth (2007) are used. sewma.q.crit determines the critical values (similar to alarm limits) for given in-control RL quantile L0 at level alpha by applying secant rule and using sewma.sf(). In case of sided="two" and mode="unbiased" a two-dimensional secant rule is applied that also ensures that the minimum of the cdf for given standard deviation is attained at sigma0.

Value

Returns a single value which resembles the RL quantile of order alpha and the lower and upper control limit cl and cu, respectively.

Author(s)

Sven Knoth

References

H.-J. Mittag and D. Stemann and B. Tewes (1998), EWMA-Karten zur \"Uberwachung der Streuung von Qualit\"atsmerkmalen, Allgemeines Statistisches Archiv 82, 327-338,

C. A. Acosta-Mej\'ia and J. J. Pignatiello Jr. and B. V. Rao (1999), A comparison of control charting procedures for monitoring process dispersion, IIE Transactions 31, 569-579.

S. Knoth (2005), Accurate ARL computation for EWMA-S^2 control charts, Statistics and Computing 15, 341-352.

S. Knoth (2007), Accurate ARL calculation for EWMA control charts monitoring simultaneously normal mean and variance, Sequential Analysis 26, 251-264.

S. Knoth (2010), Control Charting Normal Variance – Reflections, Curiosities, and Recommendations, in Frontiers in Statistical Quality Control 9, H.-J. Lenz and P.-T. Wilrich (Eds.), Physica Verlag, Heidelberg, Germany, 3-18.

K.-H. Waldmann (1986), Bounds for the distribution of the run length of geometric moving average charts, Appl. Statist. 35, 151-158.

See Also

sewma.arl for calculation of ARL of variance charts and sewma.sf for the RL survival function.

Examples

## will follow

Compute RL quantiles of EWMA (variance charts) control charts under pre-run uncertainty

Description

Computation of quantiles of the Run Length (RL) for EWMA control charts monitoring normal variance.

Usage

sewma.q.prerun(l,cl,cu,sigma,df1,df2,alpha,hs=1,sided="upper",
r=40,qm=30,qm.sigma=30,truncate=1e-10)

sewma.q.crit.prerun(l,L0,alpha,df1,df2,sigma0=1,cl=NULL,cu=NULL,hs=1,
sided="upper",mode="fixed",r=40, qm=30,qm.sigma=30,truncate=1e-10,
tail_approx=TRUE,c.error=1e-10,a.error=1e-9)

Arguments

Details

Instead of the popular ARL (Average Run Length) quantiles of the EWMA stopping time (Run Length) are determined. The algorithm is based on Waldmann's survival function iteration procedure. Thereby the ideas presented in Knoth (2007) are used. sewma.q.crit.prerun determines the critical values (similar to alarm limits) for given in-control RL quantile L0 at level alpha by applying secant rule and using sewma.sf(). In case of sided="two" and mode="unbiased" a two-dimensional secant rule is applied that also ensures that the minimum of the cdf for given standard deviation is attained at sigma0.

Value

Returns a single value which resembles the RL quantile of order alpha and the lower and upper control limit cl and cu, respectively.

Author(s)

Sven Knoth

References

S. Knoth (2007), Accurate ARL calculation for EWMA control charts monitoring simultaneously normal mean and variance, Sequential Analysis 26, 251-264.

K.-H. Waldmann (1986), Bounds for the distribution of the run length of geometric moving average charts, Appl. Statist. 35, 151-158.

See Also

sewma.q and sewma.q.crit for the version w/o pre-run uncertainty.

Examples

## will follow

Compute the survival function of EWMA run length

Description

Computation of the survival function of the Run Length (RL) for EWMA control charts monitoring normal variance.

Usage

sewma.sf(n, l, cl, cu, sigma, df, hs=1, sided="upper", r=40, qm=30)

Arguments

Details

The survival function P(L>n) and derived from it also the cdf P(L<=n) and the pmf P(L=n) illustrate the distribution of the EWMA run length. For large n the geometric tail could be exploited. That is, with reasonable large n the complete distribution is characterized. The algorithm is based on Waldmann's survival function iteration procedure and on results in Knoth (2007).

Value

Returns a vector which resembles the survival function up to a certain point.

Author(s)

Sven Knoth

References

S. Knoth (2007), Accurate ARL calculation for EWMA control charts monitoring simultaneously normal mean and variance, Sequential Analysis 26, 251-264.

K.-H. Waldmann (1986), Bounds for the distribution of the run length of geometric moving average charts, Appl. Statist. 35, 151-158.

See Also

sewma.arl for zero-state ARL computation of variance EWMA control charts.

Examples

## will follow

Compute the survival function of EWMA run length

Description

Computation of the survival function of the Run Length (RL) for EWMA control charts monitoring normal variance.

Usage

sewma.sf.prerun(n, l, cl, cu, sigma, df1, df2, hs=1, sided="upper",
qm=30, qm.sigma=30, truncate=1e-10, tail_approx=TRUE)

Arguments

Details

The survival function P(L>n) and derived from it also the cdf P(L<=n) and the pmf P(L=n) illustrate the distribution of the EWMA run length. For large n the geometric tail could be exploited. That is, with reasonable large n the complete distribution is characterized. The algorithm is based on Waldmann's survival function iteration procedure and on results in Knoth (2007)...

Value

Returns a vector which resembles the survival function up to a certain point.

Author(s)

Sven Knoth

References

S. Knoth (2007), Accurate ARL calculation for EWMA control charts monitoring simultaneously normal mean and variance, Sequential Analysis 26, 251-264.

K.-H. Waldmann (1986), Bounds for the distribution of the run length of geometric moving average charts, Appl. Statist. 35, 151-158.

See Also

sewma.sf for the RL survival function of EWMA control charts w/o pre-run uncertainty.

Examples

## will follow

Compute ARLs of Poisson TEWMA control charts

Description

Computation of the (zero-state) Average Run Length (ARL) at given Poisson mean mu.

Usage

tewma.arl(lambda, k, lk, uk, mu, z0, rando=FALSE, gl=0, gu=0)

Arguments

Details

A new idea of applying EWMA smoothing to count data. Here, the thinning operation is applied to independent Poisson variates is performed. Moreover, the original thinning principle is expanded to multiples of one over k to allow finer grids and finally better detection perfomance. It is highly recommended to read the corresponding paper (see below).

Value

Return single value which resemble the ARL.

Author(s)

Sven Knoth

References

M. C. Morais, C. H. Weiss, S. Knoth (2019), A thinning-based EWMA chart to monitor counts, submitted.

See Also

later.

Examples

# MWK (2018)
lambda <- 0.1 # (T)EWMA smoothing constant
mu0 <- 5 # in-control mean
k <- 10 # resolution
z0 <- round(k*mu0) # starting value of (T)EWMA sequence
# (i) without randomisation
lk <- 28
uk <- 75
L0 <- tewma.arl(lambda, k, lk, uk, mu0, z0)
# should be 501.9703
# (ii) with randomisation
uk <- 76 # lk is not changed
gl <- 0.5446310
gu <- 0.1375617
L0 <- tewma.arl(lambda, k, lk, uk, mu0, z0, rando=TRUE, gl=gl, gu=gu)
# should be 500

Two-sided tolerance limit factors

Description

For constructing tolerance intervals, which cover a given proportion p of a normal distribution with unknown mean and variance with confidence 1-\alpha, one needs to calculate the so-called tolerance limit factors k. These values are computed for a given sample size n.

Usage

tol.lim.fac(n,p,a,mode="WW",m=30)

Arguments

Details

tol.lim.fac determines tolerance limits factors k by means of the fast and simple approximation due to Wald/Wolfowitz (1946) and of Gauss-Legendre quadrature like Odeh/Owen (1980), respectively, who used in fact the Simpson Rule. Then, by \bar x \pm k \cdot s one can build the tolerance intervals which cover at least proportion p of a normal distribution for given confidence level of 1-\alpha. \bar x and s stand for the sample mean and the sample standard deviation, respectively.

Value

Returns a single value which resembles the tolerance limit factor.

Author(s)

Sven Knoth

References

A. Wald, J. Wolfowitz (1946), Tolerance limits for a normal distribution, Annals of Mathematical Statistics 17, 208-215.

R. E. Odeh, D. B. Owen (1980), Tables for Normal Tolerance Limits, Sampling Plans, and Screening, Marcel Dekker, New York.

See Also

qnorm for the ”asymptotic” case – cf. second example.

Examples

n <- 2:10
p <- .95
a <- .05
kWW <- sapply(n,p=p,a=a,tol.lim.fac)
kEX <- sapply(n,p=p,a=a,mode="exact",tol.lim.fac)
print(cbind(n,kWW,kEX),digits=4)
## Odeh/Owen (1980), page 98, in Table 3.4.1
##  n  factor k
##  2  36.519
##  3   9.789
##  4   6.341
##  5   5.077
##  6   4.422
##  7   4.020
##  8   3.746
##  9   3.546
## 10   3.393

## n -> infty
n <- 10^{1:7}
p <- .95
a <- .05
kEX <- round(sapply(n,p=p,a=a,mode="exact",tol.lim.fac),digits=4)
kEXinf <- round(qnorm(1-a/2),digits=4)
print(rbind(cbind(n,kEX),c("infinity",kEXinf)),quote=FALSE)

Compute ARLs of EWMA residual control charts

Description

Computation of the (zero-state) Average Run Length (ARL) for EWMA residual control charts monitoring normal mean, variance, or mean and variance simultaneously. Additionally, the probability of misleading signals (PMS) is calculated.

Usage

x.res.ewma.arl(l, c, mu, alpha=0, n=5, hs=0, r=40)

s.res.ewma.arl(l, cu, sigma, mu=0, alpha=0, n=5, hs=1, r=40, qm=30)

xs.res.ewma.arl(lx, cx, ls, csu, mu, sigma, alpha=0,
n=5, hsx=0, rx=40, hss=1, rs=40, qm=30)

xs.res.ewma.pms(lx, cx, ls, csu, mu, sigma, type="3",
alpha=0, n=5, hsx=0, rx=40, hss=1, rs=40, qm=30)

Arguments

Details

The above list of functions provides the application of algorithms developed for iid data to the residual case. To be more precise, the underlying model is a sequence of normally distributed batches with size n with autocorrelation within the batch and independence between the batches (see also the references below). It is restricted to the classical EWMA chart types, that is two-sided for the mean, upper charts for the variance, and all equipped with fixed limits. The autocorrelation is modeled by an AR(1) process with parameter alpha. Additionally, with xs.res.ewma.pms the probability of misleading signals (PMS) of type is calculated. This is offered exclusively in this small collection so that for iid data this function has to be used too (with alpha=0).

Value

Return single values which resemble the ARL and the PMS, respectively.

Author(s)

Sven Knoth

References

S. Knoth, M. C. Morais, A. Pacheco, W. Schmid (2009), Misleading Signals in Simultaneous Residual Schemes for the Mean and Variance of a Stationary Process, Commun. Stat., Theory Methods 38, 2923-2943.

S. Knoth, W. Schmid, A. Schoene (2001), Simultaneous Shewhart-Type Charts for the Mean and the Variance of a Time Series, Frontiers of Statistical Quality Control 6, A. Lenz, H.-J. & Wilrich, P.-T. (Eds.), 6, 61-79.

S. Knoth, W. Schmid (2002) Monitoring the mean and the variance of a stationary process, Statistica Neerlandica 56, 77-100.

See Also

xewma.arl, sewma.arl, and xsewma.arl as more elaborated functions in the iid case.

Examples

## Not run: 
## S. Knoth, W. Schmid (2002)

cat("\nFragments of Table 2 (n=5, lambda.1=lambda.2)\n")

lambdas <- c(.5, .25, .1, .05)
L0 <- 500
n <- 5

crit <- NULL
for ( lambda in lambdas ) {
  cs <- xsewma.crit(lambda, lambda, L0, n-1) 
  x.e <- round(cs[1], digits=4)
  names(x.e) <- NULL
  s.e <- round((cs[3]-1) * sqrt((2-lambda)/lambda)*sqrt((n-1)/2), digits=4)
  names(s.e) <- NULL
  crit <- rbind(crit, data.frame(lambda, x.e, s.e))
}


## orinal values are (Markov chain approximation with 50 states)
# lambda x.e    s.e
#   0.50 3.2765 4.6439
#   0.25 3.2168 4.0149
#   0.10 3.0578 3.3376
#   0.05 2.8817 2.9103

print(crit)


cat("\nFragments of Table 4 (n=5, lambda.1=lambda.2=0.1)\n\n")

lambda <- .1
# the algorithm used in Knoth/Schmid is less accurate -- proceed with their values
cx <- x.e <- 3.0578
s.e <- 3.3376
csu <- 1 + s.e * sqrt(lambda/(2-lambda))*sqrt(2/(n-1))

alpha <- .3

a.values <- c((0:6)/4, 2)
d.values <- c(1 + (0:5)/10, 1.75 , 2)

arls <- NULL
for ( delta in d.values ) {
  row <- NULL
  for ( mu in a.values ) {
    arl <- round(xs.res.ewma.arl(lambda, cx, lambda, csu, mu*sqrt(n), delta, alpha=alpha, n=n),
                 digits=2)
    names(arl) <- NULL
    row <- c(row, arl)   
  }
  arls <- rbind(arls, data.frame(t(row)))
}
names(arls) <- a.values
rownames(arls) <- d.values

## orinal values are (now Monte-Carlo with 10^6 replicates)
#          0  0.25   0.5 0.75    1 1.25  1.5    2
#1    502.44 49.50 14.21 7.93 5.53 4.28 3.53 2.65
#1.1   73.19 32.91 13.33 7.82 5.52 4.29 3.54 2.66
#1.2   24.42 18.88 11.37 7.44 5.42 4.27 3.54 2.67
#1.3   13.11 11.83  9.09 6.74 5.18 4.17 3.50 2.66
#1.4    8.74  8.31  7.19 5.89 4.81 4.00 3.41 2.64
#1.5    6.50  6.31  5.80 5.08 4.37 3.76 3.28 2.59
#1.75   3.94  3.90  3.78 3.59 3.35 3.09 2.83 2.40
#2      2.85  2.84  2.80 2.73 2.63 2.51 2.39 2.14

print(arls)


## S. Knoth, M. C. Morais, A. Pacheco, W. Schmid (2009)

cat("\nFragments of Table 5 (n=5, lambda=0.1)\n\n")

d.values <- c(1.02, 1 + (1:5)/10, 1.75 , 2)

arl.x <- arl.s <- arl.xs <- PMS.3 <- NULL
for ( delta in d.values ) {
  arl.x  <- c(arl.x,  round(x.res.ewma.arl(lambda, cx/delta, 0, n=n),
                            digits=3))
  arl.s  <- c(arl.s,  round(s.res.ewma.arl(lambda, csu, delta, n=n),
                            digits=3))
  arl.xs <- c(arl.xs, round(xs.res.ewma.arl(lambda, cx, lambda, csu, 0, delta, n=n),
                            digits=3))
  PMS.3  <- c(PMS.3,  round(xs.res.ewma.pms(lambda, cx, lambda, csu, 0, delta, n=n),
                            digits=6))
}

## orinal values are (Markov chain approximation)
# delta   arl.x   arl.s  arl.xs PMS.3
#  1.02 833.086 518.935 323.324 0.381118
#  1.10 454.101  84.208  73.029 0.145005
#  1.20 250.665  25.871  24.432 0.071024
#  1.30 157.343  13.567  13.125 0.047193
#  1.40 108.112   8.941   8.734 0.035945
#  1.50  79.308   6.614   6.493 0.029499
#  1.75  44.128   3.995   3.942 0.021579
#  2.00  28.974   2.887   2.853 0.018220

print(cbind(delta=d.values, arl.x, arl.s, arl.xs, PMS.3))


cat("\nFragments of Table 6 (n=5, lambda=0.1)\n\n")

alphas <- c(-0.9, -0.5, -0.3, 0, 0.3, 0.5, 0.9)
deltas <- c(0.05, 0.25, 0.5, 0.75, 1, 1.25, 1.5, 2)

PMS.4 <- NULL
for ( ir in 1:length(deltas) ) {
  mu <- deltas[ir]*sqrt(n)
  pms <- NULL
  for ( alpha in alphas ) {
    pms <- c(pms, round(xs.res.ewma.pms(lambda, cx, lambda, csu, mu, 1, type="4", alpha=alpha, n=n),
                        digits=6))
  }
  PMS.4 <- rbind(PMS.4, data.frame(delta=deltas[ir], t(pms)))
}
names(PMS.4) <- c("delta", alphas)
rownames(PMS.4) <- NULL

## orinal values are (Markov chain approximation)
#  delta     -0.9     -0.5     -0.3        0      0.3      0.5      0.9
#   0.05 0.055789 0.224521 0.279842 0.342805 0.391299 0.418915 0.471386
#   0.25 0.003566 0.009522 0.014580 0.025786 0.044892 0.066584 0.192023
#   0.50 0.002994 0.001816 0.002596 0.004774 0.009259 0.015303 0.072945
#   0.75 0.006967 0.000703 0.000837 0.001529 0.003400 0.006424 0.046602
#   1.00 0.005098 0.000402 0.000370 0.000625 0.001589 0.003490 0.039978
#   1.25 0.000084 0.000266 0.000202 0.000300 0.000867 0.002220 0.039773
#   1.50 0.000000 0.000256 0.000120 0.000163 0.000531 0.001584 0.042734
#   2.00 0.000000 0.000311 0.000091 0.000056 0.000259 0.001029 0.054543

print(PMS.4)

## End(Not run)

Compute ARLs of CUSUM control charts under drift

Description

Computation of the (zero-state and other) Average Run Length (ARL) under drift for one-sided CUSUM control charts monitoring normal mean.

Usage

xDcusum.arl(k, h, delta, hs = 0, sided = "one",
    mode = "Gan", m = NULL, q = 1, r = 30, with0 = FALSE)

Arguments

Details

Based on Gan (1991) or Knoth (2003), the ARL is calculated for one-sided CUSUM control charts under drift. In case of Gan's framework, the usual ARL function with mu=m*delta is determined and recursively via m-1, m-2, ... 1 (or 0) the drift ARL determined. The framework of Knoth allows to calculate ARLs for varying parameters, such as control limits and distributional parameters. For details see the cited papers. Note that two-sided CUSUM charts under drift are difficult to treat.

Value

Returns a single value which resembles the ARL.

Author(s)

Sven Knoth

References

F. F. Gan (1992), CUSUM control charts under linear drift, Statistician 41, 71-84.

F. F. Gan (1996), Average Run Lengths for Cumulative Sum control chart under linear trend, Applied Statistics 45, 505-512.

S. Knoth (2003), EWMA schemes with non-homogeneous transition kernels, Sequential Analysis 22, 241-255.

S. Knoth (2012), More on Control Charting under Drift, in: Frontiers in Statistical Quality Control 10, H.-J. Lenz, W. Schmid and P.-T. Wilrich (Eds.), Physica Verlag, Heidelberg, Germany, 53-68.

C. Zou, Y. Liu and Z. Wang (2009), Comparisons of control schemes for monitoring the means of processes subject to drifts, Metrika 70, 141-163.

See Also

xcusum.arl and xcusum.ad for zero-state and steady-state ARL computation of CUSUM control charts for the classical step change model.

Examples

## Gan (1992)
## Table 1
## original values are
#   deltas  arl
#   0.0001  475
#   0.0005  261
#   0.0010  187
#   0.0020  129
#   0.0050  76.3
#   0.0100  52.0
#   0.0200  35.2
#   0.0500  21.4
#   0.1000  15.0
#   0.5000  6.95
#   1.0000  5.16
#   3.0000  3.30
## Not run: k <- .25
h <- 8
r <- 50
DxDcusum.arl <- Vectorize(xDcusum.arl, "delta")
deltas <- c(0.0001, 0.0005, 0.001, 0.002, 0.005, 0.01, 0.02, 0.05, 0.1, 0.5, 1, 3)
arl.like.Gan   <-
  round(DxDcusum.arl(k, h, deltas, r=r, with0=TRUE), digits=2)
arl.like.Knoth <-
  round(DxDcusum.arl(k, h, deltas, r=r, mode="Knoth", with0=TRUE), digits=2)
data.frame(deltas, arl.like.Gan, arl.like.Knoth)
## End(Not run)

## Zou et al. (2009)
## Table 1
## original values are
#  delta   arl1  arl2  arl3
#  0           ~ 1730
#  0.0005  345   412   470
#  0.001   231   275   317
#  0.005   86.6  98.6  112
#  0.01    56.9  61.8  69.3
#  0.05    22.6  21.6  22.7
#  0.1     15.4  14.7  14.2
#  0.5     6.60  5.54  5.17
#  1.0     4.63  3.80  3.45
#  2.0     3.17  2.67  2.32
#  3.0     2.79  2.04  1.96
#  4.0     2.10  1.98  1.74
## Not run: 
k1 <- 0.25
k2 <- 0.5
k3 <- 0.75
h1 <- 9.660
h2 <- 5.620
h3 <- 3.904
deltas <- c(0.0005, 0.001, 0.005, 0.01, 0.05, 0.1, 0.5, 1:4)
arl1 <- c(round(xcusum.arl(k1, h1, 0, r=r), digits=2),
          round(DxDcusum.arl(k1, h1, deltas, r=r), digits=2))
arl2 <- c(round(xcusum.arl(k2, h2, 0), digits=2),
          round(DxDcusum.arl(k2, h2, deltas, r=r), digits=2))
arl3 <- c(round(xcusum.arl(k3, h3, 0, r=r), digits=2),
          round(DxDcusum.arl(k3, h3, deltas, r=r), digits=2))
data.frame(delta=c(0, deltas), arl1, arl2, arl3)
## End(Not run)

Compute ARLs of EWMA control charts under drift

Description

Computation of the (zero-state and other) Average Run Length (ARL) under drift for different types of EWMA control charts monitoring normal mean.

Usage

xDewma.arl(l, c, delta, zr = 0, hs = 0, sided = "one", limits = "fix",
    mode = "Gan", m = NULL, q = 1, r = 40, with0 = FALSE)

Arguments

Details

Based on Gan (1991) or Knoth (2003), the ARL is calculated for EWMA control charts under drift. In case of Gan's framework, the usual ARL function with mu=m*delta is determined and recursively via m-1, m-2, ... 1 (or 0) the drift ARL determined. The framework of Knoth allows to calculate ARLs for varying parameters, such as control limits and distributional parameters. For details see the cited papers.

Value

Returns a single value which resembles the ARL.

Author(s)

Sven Knoth

References

F. F. Gan (1991), EWMA control chart under linear drift, J. Stat. Comput. Simulation 38, 181-200.

L. A. Aerne, C. W. Champ and S. E. Rigdon (1991), Evaluation of control charts under linear trend, Commun. Stat., Theory Methods 20, 3341-3349.

S. Knoth (2003), EWMA schemes with non-homogeneous transition kernels, Sequential Analysis 22, 241-255.

H. M. Fahmy and E. A. Elsayed (2006), Detection of linear trends in process mean, International Journal of Production Research 44, 487-504.

S. Knoth (2012), More on Control Charting under Drift, in: Frontiers in Statistical Quality Control 10, H.-J. Lenz, W. Schmid and P.-T. Wilrich (Eds.), Physica Verlag, Heidelberg, Germany, 53-68.

C. Zou, Y. Liu and Z. Wang (2009), Comparisons of control schemes for monitoring the means of processes subject to drifts, Metrika 70, 141-163.

See Also

xewma.arl and xewma.ad for zero-state and steady-state ARL computation of EWMA control charts for the classical step change model.

Examples

## Not run: 
DxDewma.arl <- Vectorize(xDewma.arl, "delta")
## Gan (1991)
## Table 1
## original values are
#  delta   arlE1  arlE2  arlE3
#  0       500    500    500
#  0.0001  482    460    424
#  0.0010  289    231    185
#  0.0020  210    162    129
#  0.0050  126     94.6   77.9
#  0.0100   81.7   61.3   52.7
#  0.0500   27.5   21.8   21.9
#  0.1000   17.0   14.2   15.3
#  1.0000    4.09   4.28   5.25
#  3.0000    2.60   2.90   3.43
#
lambda1 <- 0.676
lambda2 <- 0.242
lambda3 <- 0.047
h1 <- 2.204/sqrt(lambda1/(2-lambda1))
h2 <- 1.111/sqrt(lambda2/(2-lambda2))
h3 <- 0.403/sqrt(lambda3/(2-lambda3))
deltas <- c(.0001, .001, .002, .005, .01, .05, .1, 1, 3)
arlE10 <- round(xewma.arl(lambda1, h1, 0, sided="two"), digits=2)
arlE1 <- c(arlE10, round(DxDewma.arl(lambda1, h1, deltas, sided="two", with0=TRUE),
                         digits=2))
arlE20 <- round(xewma.arl(lambda2, h2, 0, sided="two"), digits=2)
arlE2 <- c(arlE20, round(DxDewma.arl(lambda2, h2, deltas, sided="two", with0=TRUE),
                         digits=2))
arlE30 <- round(xewma.arl(lambda3, h3, 0, sided="two"), digits=2)
arlE3 <- c(arlE30, round(DxDewma.arl(lambda3, h3, deltas, sided="two", with0=TRUE),
                         digits=2))
data.frame(delta=c(0, deltas), arlE1, arlE2, arlE3)

## do the same with more digits for the alarm threshold
L0 <- 500
h1 <- xewma.crit(lambda1, L0, sided="two")
h2 <- xewma.crit(lambda2, L0, sided="two")
h3 <- xewma.crit(lambda3, L0, sided="two")
lambdas <- c(lambda1, lambda2, lambda3)
hs <- c(h1, h2, h3) * sqrt(lambdas/(2-lambdas))
hs
# compare with Gan's values 2.204, 1.111, 0.403
round(hs, digits=3)

arlE10 <- round(xewma.arl(lambda1, h1, 0, sided="two"), digits=2)
arlE1 <- c(arlE10, round(DxDewma.arl(lambda1, h1, deltas, sided="two", with0=TRUE),
                         digits=2))
arlE20 <- round(xewma.arl(lambda2, h2, 0, sided="two"), digits=2)
arlE2 <- c(arlE20, round(DxDewma.arl(lambda2, h2, deltas, sided="two", with0=TRUE),
                         digits=2))
arlE30 <- round(xewma.arl(lambda3, h3, 0, sided="two"), digits=2)
arlE3 <- c(arlE30, round(DxDewma.arl(lambda3, h3, deltas, sided="two", with0=TRUE),
                         digits=2))
data.frame(delta=c(0, deltas), arlE1, arlE2, arlE3)

## Aerne et al. (1991) -- two-sided EWMA
## Table I (continued)
## original numbers are
#     delta  arlE1  arlE2  arlE3
#  0.000000  465.0  465.0  465.0
#  0.005623  77.01  85.93  102.68
#  0.007499  64.61  71.78  85.74
#  0.010000  54.20  59.74  71.22
#  0.013335  45.20  49.58  58.90
#  0.017783  37.76  41.06  48.54
#  0.023714  31.54  33.95  39.87
#  0.031623  26.36  28.06  32.68
#  0.042170  22.06  23.19  26.73
#  0.056234  18.49  19.17  21.84
#  0.074989  15.53  15.87  17.83
#  0.100000  13.07  13.16  14.55
#  0.133352  11.03  10.94  11.88
#  0.177828   9.33   9.12   9.71
#  0.237137   7.91   7.62   7.95
#  0.316228   6.72   6.39   6.52
#  0.421697   5.72   5.38   5.37
#  0.562341   4.88   4.54   4.44
#  0.749894   4.18   3.84   3.68
#  1.000000   3.58   3.27   3.07
#
lambda1 <- .133
lambda2 <- .25
lambda3 <- .5
cE1 <- 2.856
cE2 <- 2.974
cE3 <- 3.049
deltas <- 10^(-(18:0)/8)
arlE10 <- round(xewma.arl(lambda1, cE1, 0, sided="two"), digits=2)
arlE1 <- c(arlE10, round(DxDewma.arl(lambda1, cE1, deltas, sided="two"), digits=2))
arlE20 <- round(xewma.arl(lambda2, cE2, 0, sided="two"), digits=2)
arlE2 <- c(arlE20, round(DxDewma.arl(lambda2, cE2, deltas, sided="two"), digits=2))
arlE30 <- round(xewma.arl(lambda3, cE3, 0, sided="two"), digits=2)
arlE3 <- c(arlE30, round(DxDewma.arl(lambda3, cE3, deltas, sided="two"), digits=2))
data.frame(delta=c(0, round(deltas, digits=6)), arlE1, arlE2, arlE3)


## Fahmy/Elsayed (2006) -- two-sided EWMA
## Table 4 (Monte Carlo results, 10^4 replicates, change point at t=51!)
## original numbers are
#   delta     arl  s.e.
#   0.00  365.749  3.598
#   0.10   12.971  0.029
#   0.25    7.738  0.015
#   0.50    5.312  0.009
#   0.75    4.279  0.007
#   1.00    3.680  0.006
#   1.25    3.271  0.006
#   1.50    2.979  0.005
#   1.75    2.782  0.004
#   2.00    2.598  0.005
#
lambda <- 0.1
cE <- 2.7
deltas <- c(.1, (1:8)/4)
arlE1 <- c(round(xewma.arl(lambda, cE, 0, sided="two"), digits=3),
           round(DxDewma.arl(lambda, cE, deltas, sided="two"), digits=3))
arlE51 <- c(round(xewma.arl(lambda, cE, 0, sided="two", q=51)[51], digits=3),
     round(DxDewma.arl(lambda, cE, deltas, sided="two", mode="Knoth", q=51),
           digits=3))
data.frame(delta=c(0, deltas), arlE1, arlE51)

## additional Monte Carlo results with 10^8 replicates
#   delta   arl.q=1   s.e.    arl.q=51  s.e.
#   0.00    368.910   0.036   361.714   0.038
#   0.10     12.986   0.000    12.781   0.000
#   0.25      7.758   0.000     7.637   0.000
#   0.50      5.318   0.000     5.235   0.000
#   0.75      4.285   0.000     4.218   0.000
#   1.00      3.688   0.000     3.628   0.000
#   1.25      3.274   0.000     3.233   0.000
#   1.50      2.993   0.000     2.942   0.000
#   1.75      2.808   0.000     2.723   0.000
#   2.00      2.616   0.000     2.554   0.000

## Zou et al. (2009) -- one-sided EWMA
## Table 1
## original values are
#  delta   arl1  arl2  arl3
#  0           ~ 1730
#  0.0005  317   377   440
#  0.001   215   253   297
#  0.005   83.6  92.6  106
#  0.01    55.6  58.8  66.1
#  0.05    22.6  21.1  22.0
#  0.1     15.5  13.9  13.8
#  0.5     6.65  5.56  5.09
#  1.0     4.67  3.83  3.43
#  2.0     3.21  2.74  2.32
#  3.0     2.86  2.06  1.98
#  4.0     2.14  2.00  1.83
l1 <- 0.03479
l2 <- 0.11125
l3 <- 0.23052
c1 <- 2.711
c2 <- 3.033
c3 <- 3.161
zr <- -6
r  <- 50
deltas <- c(0.0005, 0.001, 0.005, 0.01, 0.05, 0.1, 0.5, 1:4)
arl1 <- c(round(xewma.arl(l1, c1, 0, zr=zr, r=r), digits=2),
          round(DxDewma.arl(l1, c1, deltas, zr=zr, r=r), digits=2))
arl2 <- c(round(xewma.arl(l2, c2, 0, zr=zr), digits=2),
          round(DxDewma.arl(l2, c2, deltas, zr=zr, r=r), digits=2))
arl3 <- c(round(xewma.arl(l3, c3, 0, zr=zr, r=r), digits=2),
          round(DxDewma.arl(l3, c3, deltas, zr=zr, r=r), digits=2))
data.frame(delta=c(0, deltas), arl1, arl2, arl3)

## End(Not run)

Compute ARLs of Shiryaev-Roberts schemes under drift

Description

Computation of the (zero-state and other) Average Run Length (ARL) under drift for Shiryaev-Roberts schemes monitoring normal mean.

Usage

xDgrsr.arl(k, g, delta, zr = 0, hs = NULL, sided = "one", m = NULL,
mode = "Gan", q = 1, r = 30, with0 = FALSE)

Arguments

Details

Based on Gan (1991) or Knoth (2003), the ARL is calculated for Shiryaev-Roberts schemes under drift. In case of Gan's framework, the usual ARL function with mu=m*delta is determined and recursively via m-1, m-2, ... 1 (or 0) the drift ARL determined. The framework of Knoth allows to calculate ARLs for varying parameters, such as control limits and distributional parameters. For details see the cited papers.

Value

Returns a single value which resembles the ARL.

Author(s)

Sven Knoth

References

F. F. Gan (1991), EWMA control chart under linear drift, J. Stat. Comput. Simulation 38, 181-200.

S. Knoth (2003), EWMA schemes with non-homogeneous transition kernels, Sequential Analysis 22, 241-255.

S. Knoth (2012), More on Control Charting under Drift, in: Frontiers in Statistical Quality Control 10, H.-J. Lenz, W. Schmid and P.-T. Wilrich (Eds.), Physica Verlag, Heidelberg, Germany, 53-68.

C. Zou, Y. Liu and Z. Wang (2009), Comparisons of control schemes for monitoring the means of processes subject to drifts, Metrika 70, 141-163.

See Also

xewma.arl and xewma.ad for zero-state and steady-state ARL computation of EWMA control charts for the classical step change model.

Examples

## Not run: 
## Monte Carlo example with 10^8 replicates
#   delta      arl    s.e.
#   0.0001 381.8240   0.0304
#   0.0005 238.4630   0.0148
#   0.001  177.4061   0.0097
#   0.002  125.9055   0.0061
#   0.005   75.7574   0.0031
#   0.01    50.2203   0.0018
#   0.02    32.9458   0.0011
#   0.05    18.9213   0.0005
#   0.1     12.6054   0.0003
#   0.5      5.2157   0.0001
#   1        3.6537   0.0001
#   3        2.0289   0.0000
k <- .5
L0 <- 500
zr <- -7
r <- 50
g <- xgrsr.crit(k, L0, zr=zr, r=r)
DxDgrsr.arl <- Vectorize(xDgrsr.arl, "delta")
deltas <- c(0.0001, 0.0005, 0.001, 0.002, 0.005, 0.01, 0.02, 0.05, 0.1, 0.5, 1, 3)
arls <- round(DxDgrsr.arl(k, g, deltas, zr=zr, r=r), digits=4)
data.frame(deltas, arls)

## End(Not run)

Compute ARLs of Shewhart control charts with and without runs rules under drift

Description

Computation of the zero-state Average Run Length (ARL) under drift for Shewhart control charts with and without runs rules monitoring normal mean.

Usage

xDshewhartrunsrules.arl(delta, c = 1, m = NULL, type = "12")

xDshewhartrunsrulesFixedm.arl(delta, c = 1, m = 100, type = "12")

Arguments

Details

Based on Gan (1991), the ARL is calculated for Shewhart control charts with and without runs rules under drift. The usual ARL function with mu=m*delta is determined and recursively via m-1, m-2, ... 1 (or 0) the drift ARL determined. xDshewhartrunsrulesFixedm.arl is the actual work horse, while xDshewhartrunsrules.arl provides a convenience wrapper. Note that Aerne et al. (1991) deployed a method that is quite similar to Gan's algorithm. For type see the help page of xshewhartrunsrules.arl.

Value

Returns a single value which resembles the ARL.

Author(s)

Sven Knoth

References

F. F. Gan (1991), EWMA control chart under linear drift, J. Stat. Comput. Simulation 38, 181-200.

L. A. Aerne, C. W. Champ and S. E. Rigdon (1991), Evaluation of control charts under linear trend, Commun. Stat., Theory Methods 20, 3341-3349.

See Also

xshewhartrunsrules.arl for zero-state ARL computation of Shewhart control charts with and without runs rules for the classical step change model.

Examples

## Aerne et al. (1991)
## Table I (continued)
## original numbers are
#     delta arl1of1 arl2of3 arl4of5  arl10
#  0.005623  136.67  120.90  105.34 107.08
#  0.007499  114.98  101.23   88.09  89.94
#  0.010000   96.03   84.22   73.31  75.23
#  0.013335   79.69   69.68   60.75  62.73
#  0.017783   65.75   57.38   50.18  52.18
#  0.023714   53.99   47.06   41.33  43.35
#  0.031623   44.15   38.47   33.99  36.00
#  0.042170   35.97   31.36   27.91  29.90
#  0.056234   29.21   25.51   22.91  24.86
#  0.074989   23.65   20.71   18.81  20.70
#  0.100000   19.11   16.79   15.45  17.29
#  0.133352   15.41   13.61   12.72  14.47
#  0.177828   12.41   11.03   10.50  12.14
#  0.237137    9.98    8.94    8.71  10.18
#  0.316228    8.02    7.25    7.26   8.45
#  0.421697    6.44    5.89    6.09   6.84
#  0.562341    5.17    4.80    5.15   5.48
#  0.749894    4.16    3.92    4.36   4.39
#  1.000000    3.35    3.22    3.63   3.52
c1of1 <- 3.069/3
c2of3 <- 2.1494/2
c4of5 <- 1.14
c10   <- 3.2425/3
DxDshewhartrunsrules.arl <- Vectorize(xDshewhartrunsrules.arl, "delta")
deltas <- 10^(-(18:0)/8)
arl1of1 <-
round(DxDshewhartrunsrules.arl(deltas, c=c1of1, type="1"), digits=2)
arl2of3 <-
round(DxDshewhartrunsrules.arl(deltas, c=c2of3, type="12"), digits=2)
arl4of5 <-
round(DxDshewhartrunsrules.arl(deltas, c=c4of5, type="13"), digits=2)
arl10 <- 
round(DxDshewhartrunsrules.arl(deltas, c=c10, type="SameSide10"), digits=2)
data.frame(delta=round(deltas, digits=6), arl1of1, arl2of3, arl4of5, arl10)

Compute steady-state ARLs of CUSUM control charts

Description

Computation of the steady-state Average Run Length (ARL) for different types of CUSUM control charts monitoring normal mean.

Usage

xcusum.ad(k, h, mu1, mu0 = 0, sided = "one", r = 30)

Arguments

Details

xcusum.ad determines the steady-state Average Run Length (ARL) by numerically solving the related ARL integral equation by means of the Nystroem method based on Gauss-Legendre quadrature and using the power method for deriving the largest in magnitude eigenvalue and the related left eigenfunction.

Value

Returns a single value which resembles the steady-state ARL.

Note

Be cautious in increasing the dimension parameter r for two-sided CUSUM schemes. The resulting matrix dimension is r^2 times r^2. Thus, go beyond 30 only on fast machines. This is the only case, were the package routines are based on the Markov chain approach. Moreover, the two-sided CUSUM scheme needs a two-dimensional Markov chain.

Author(s)

Sven Knoth

References

R. B. Crosier (1986), A new two-sided cumulative quality control scheme, Technometrics 28, 187-194.

See Also

xcusum.arl for zero-state ARL computation and xewma.ad for the steady-state ARL of EWMA control charts.

Examples

## comparison of zero-state (= worst case ) and steady-state performance
## for one-sided CUSUM control charts

k <- .5
h <- xcusum.crit(k,500)
mu <- c(0,.5,1,1.5,2)
arl <- sapply(mu,k=k,h=h,xcusum.arl)
ad <- sapply(mu,k=k,h=h,xcusum.ad)
round(cbind(mu,arl,ad),digits=2)

## Crosier (1986), Crosier's modified two-sided CUSUM
## He introduced the modification and evaluated it by means of
## Markov chain approximation

k <- .5
h2 <- 4
hC <- 3.73
mu <- c(0,.25,.5,.75,1,1.5,2,2.5,3,4,5)
ad2 <- sapply(mu,k=k,h=h2,sided="two",r=20,xcusum.ad)
adC <- sapply(mu,k=k,h=hC,sided="Crosier",xcusum.ad)
round(cbind(mu,ad2,adC),digits=2)

## results in the original paper are (in Table 5)
## 0.00 163.   164.
## 0.25  71.6   69.0
## 0.50  25.2   24.3
## 0.75  12.3   12.1
## 1.00   7.68   7.69
## 1.50   4.31   4.39
## 2.00   3.03   3.12
## 2.50   2.38   2.46
## 3.00   2.00   2.07
## 4.00   1.55   1.60
## 5.00   1.22   1.29

Compute ARLs of CUSUM control charts

Description

Computation of the (zero-state) Average Run Length (ARL) for different types of CUSUM control charts monitoring normal mean.

Usage

xcusum.arl(k, h, mu, hs = 0, sided = "one", method = "igl", q = 1, r = 30)

Arguments

Details

xcusum.arl determines the Average Run Length (ARL) by numerically solving the related ARL integral equation by means of the Nystroem method based on Gauss-Legendre quadrature.

Value

Returns a vector of length q which resembles the ARL and the sequence of conditional expected delays for q=1 and q>1, respectively.

Author(s)

Sven Knoth

References

A. L. Goel, S. M. Wu (1971), Determination of A.R.L. and a contour nomogram for CUSUM charts to control normal mean, Technometrics 13, 221-230.

D. Brook, D. A. Evans (1972), An approach to the probability distribution of cusum run length, Biometrika 59, 539-548.

J. M. Lucas, R. B. Crosier (1982), Fast initial response for cusum quality-control schemes: Give your cusum a headstart, Technometrics 24, 199-205.

L. C. Vance (1986), Average run lengths of cumulative sum control charts for controlling normal means, Journal of Quality Technology 18, 189-193.

K.-H. Waldmann (1986), Bounds for the distribution of the run length of one-sided and two-sided CUSUM quality control schemes, Technometrics 28, 61-67.

R. B. Crosier (1986), A new two-sided cumulative quality control scheme, Technometrics 28, 187-194.

See Also

xewma.arl for zero-state ARL computation of EWMA control charts and xcusum.ad for the steady-state ARL.

Examples

## Brook/Evans (1972), one-sided CUSUM
## Their results are based on the less accurate Markov chain approach.

k <- .5
h <- 3
round(c( xcusum.arl(k,h,0), xcusum.arl(k,h,1.5) ),digits=2)

## results in the original paper are L0 = 117.59, L1 = 3.75 (in Subsection 4.3).

## Lucas, Crosier (1982)
## (one- and) two-sided CUSUM with possible headstarts

k <- .5
h <- 4
mu <- c(0,.25,.5,.75,1,1.5,2,2.5,3,4,5)
arl1 <- sapply(mu,k=k,h=h,sided="two",xcusum.arl)
arl2 <- sapply(mu,k=k,h=h,hs=h/2,sided="two",xcusum.arl)
round(cbind(mu,arl1,arl2),digits=2)

## results in the original paper are (in Table 1)
## 0.00 168.   149.
## 0.25  74.2   62.7
## 0.50  26.6   20.1
## 0.75  13.3    8.97
## 1.00   8.38   5.29
## 1.50   4.75   2.86
## 2.00   3.34   2.01
## 2.50   2.62   1.59
## 3.00   2.19   1.32
## 4.00   1.71   1.07
## 5.00   1.31   1.01

## Vance (1986), one-sided CUSUM
## The first paper on using Nystroem method and Gauss-Legendre quadrature
## for solving the ARL integral equation (see as well Goel/Wu, 1971)

k <- 0
h <- 10
mu <- c(-.25,-.125,0,.125,.25,.5,.75,1)
round(cbind(mu,sapply(mu,k=k,h=h,xcusum.arl)),digits=2)

## results in the original paper are (in Table 1 incl. Goel/Wu (1971) results)
##  -0.25  2071.51
##  -0.125  400.28
##   0.0    124.66
##   0.125   59.30
##   0.25    36.71
##   0.50    20.37
##   0.75    14.06
##   1.00    10.75

## Waldmann (1986),
## one- and two-sided CUSUM

## one-sided case

k <- .5
h <- 3
mu <- c(-.5,0,.5)
round(sapply(mu,k=k,h=h,xcusum.arl),digits=2)

## results in the original paper are 1963, 117.4, and 17.35, resp.
## (in Tables 3, 1, and 5, resp.).

## two-sided case

k <- .6
h <- 3
round(xcusum.arl(k,h,-.2,sided="two"),digits=1)  # fits to Waldmann's setup

## result in the original paper is 65.4 (in Table 6).

## Crosier (1986), Crosier's modified two-sided CUSUM
## He introduced the modification and evaluated it by means of
## Markov chain approximation

k <- .5
h <- 3.73
mu <- c(0,.25,.5,.75,1,1.5,2,2.5,3,4,5)
round(cbind(mu,sapply(mu,k=k,h=h,sided="Crosier",xcusum.arl)),digits=2)

## results in the original paper are (in Table 3)
## 0.00 168.
## 0.25  70.7
## 0.50  25.1
## 0.75  12.5
## 1.00   7.92
## 1.50   4.49
## 2.00   3.17
## 2.50   2.49
## 3.00   2.09
## 4.00   1.60
## 5.00   1.22

## SAS/QC manual 1999
## one- and two-sided CUSUM schemes

## one-sided

k <- .25
h <- 8
mu <- 2.5
print(xcusum.arl(k,h,mu),digits=12)
print(xcusum.arl(k,h,mu,hs=.1),digits=12)

## original results are 4.1500836225 and 4.1061588131.

## two-sided

print(xcusum.arl(k,h,mu,sided="two"),digits=12)

## original result is 4.1500826715.

Compute decision intervals of CUSUM control charts

Description

Computation of the decision intervals (alarm limits) for different types of CUSUM control charts monitoring normal mean.

Usage

xcusum.crit(k, L0, mu0 = 0, hs = 0, sided = "one", r = 30)

Arguments

Details

xcusum.crit determines the decision interval (alarm limit) for given in-control ARL L0 by applying secant rule and using xcusum.arl().

Value

Returns a single value which resembles the decision interval h.

Author(s)

Sven Knoth

See Also

xcusum.arl for zero-state ARL computation.

Examples

k <- .5
incontrolARL <- c(500,5000,50000)
sapply(incontrolARL,k=k,xcusum.crit,r=10) # accuracy with 10 nodes
sapply(incontrolARL,k=k,xcusum.crit,r=20) # accuracy with 20 nodes
sapply(incontrolARL,k=k,xcusum.crit)      # accuracy with 30 nodes

Compute the CUSUM k and h for given in-control ARL L0 and out-of-control L1

Description

Computation of the reference value k and the alarm threshold h for one-sided CUSUM control charts monitoring normal mean, if the in-control ARL L0 and the out-of-control L1 are given.

Usage

xcusum.crit.L0L1(L0, L1, hs=0, sided="one", r=30, L1.eps=1e-6, k.eps=1e-8)

Arguments

Details

xcusum.crit.L0L1 determines the reference value k and the alarm threshold h for given in-control ARL L0 and out-of-control ARL L1 by applying secant rule and using xcusum.arl() and xcusum.crit(). These CUSUM design rules were firstly (and quite rarely afterwards) used by Ewan and Kemp.

Value

Returns two values which resemble the reference value k and the threshold h.

Author(s)

Sven Knoth

References

W. D. Ewan and K. W. Kemp (1960), Sampling inspection of continuous processes with no autocorrelation between successive results, Biometrika 47, 363-380.

K. W. Kemp (1962), The Use of Cumulative Sums for Sampling Inspection Schemes, Journal of the Royal Statistical Sociecty C, Applied Statistics, 10, 16-31.

See Also

xcusum.arl for zero-state ARL and xcusum.crit for threshold h computation.

Examples

## Table 2 from Ewan/Kemp (1960) -- one-sided CUSUM
#
# A.R.L. at A.Q.L.   A.R.L. at A.Q.L.     k      h
#       1000                3           1.12   2.40
#       1000                7           0.65   4.06
#        500                3           1.04   2.26
#        500                7           0.60   3.80
#        250                3           0.94   2.11
#        250                7           0.54   3.51
#
L0.set <- c(1000, 500, 250)
L1.set <- c(3, 7)
cat("\nL0\tL1\tk\th\n")
for ( L0 in L0.set ) {
  for ( L1 in L1.set ) {
    result <- round(xcusum.crit.L0L1(L0, L1), digits=2)
    cat(paste(L0, L1, result[1], result[2], sep="\t"), "\n")
  }
}
#
# two confirmation runs
xcusum.arl(0.54, 3.51, 0) # Ewan/Kemp
xcusum.arl(result[1], result[2], 0) # here
xcusum.arl(0.54, 3.51, 2*0.54) # Ewan/Kemp
xcusum.arl(result[1], result[2], 2*result[1]) # here
#
## Table II from Kemp (1962) -- two-sided CUSUM
#
#    Lr                  k
#             La=250   La=500   La=1000
#    2.5       1.05     1.17     1.27
#    3.0       0.94     1.035    1.13
#    4.0       0.78     0.85     0.92
#    5.0       0.68     0.74     0.80
#    6.0       0.60     0.655    0.71
#    7.5       0.52     0.57     0.62
#   10.0       0.43     0.48     0.52
#
L0.set <- c(250, 500, 1000)
L1.set <- c(2.5, 3:6, 7.5, 10)
cat("\nL1\tL0=250\tL0=500\tL0=1000\n")
for ( L1 in L1.set ) {
  cat(L1)
  for ( L0 in L0.set ) {
    result <- round(xcusum.crit.L0L1(L0, L1, sided="two"), digits=2)
    cat("\t", result[1])
  }
  cat("\n")
}

Compute the CUSUM reference value k for given in-control ARL and threshold h

Description

Computation of the reference value k for one-sided CUSUM control charts monitoring normal mean, if the in-control ARL L0 and the alarm threshold h are given.

Usage

xcusum.crit.L0h(L0, h, hs=0, sided="one", r=30, L0.eps=1e-6, k.eps=1e-8)

Arguments

Details

xcusum.crit.L0h determines the reference value k for given in-control ARL L0 and alarm level h by applying secant rule and using xcusum.arl(). Note that not for any combination of L0 and h a solution exists – for given L0 there is a maximal value for h to get a valid result k.

Value

Returns a single value which resembles the reference value k.

Author(s)

Sven Knoth

See Also

xcusum.arl for zero-state ARL computation.

Examples

L0 <- 100
h.max <- xcusum.crit(0, L0, 0)
hs <- (300:1)/100
hs <- hs[hs < h.max]
ks <- NULL
for ( h in hs ) ks <- c(ks, xcusum.crit.L0h(L0, h))  
k.max <- qnorm( 1 - 1/L0 )
plot(hs, ks, type="l", ylim=c(0, max(k.max, ks)), xlab="h", ylab="k")
abline(h=c(0, k.max), col="red")

Compute RL quantiles of CUSUM control charts

Description

Computation of quantiles of the Run Length (RL)for CUSUM control charts monitoring normal mean.

Usage

xcusum.q(k, h, mu, alpha, hs=0, sided="one", r=40)

Arguments

Details

Instead of the popular ARL (Average Run Length) quantiles of the CUSUM stopping time (Run Length) are determined. The algorithm is based on Waldmann's survival function iteration procedure.

Value

Returns a single value which resembles the RL quantile of order q.

Author(s)

Sven Knoth

References

K.-H. Waldmann (1986), Bounds for the distribution of the run length of one-sided and two-sided CUSUM quality control schemes, Technometrics 28, 61-67.

See Also

xcusum.arl for zero-state ARL computation of CUSUM control charts.

Examples

## Waldmann (1986), one-sided CUSUM, Table 2
## original values are 345, 82, 9

XCUSUM.Q <- Vectorize("xcusum.q", "alpha")
k <- .5
h <- 3
mu <- 0 # corresponds to Waldmann's -0.5
a.list <- c(.95, .5, .05)
rl.quantiles <- ceiling(XCUSUM.Q(k, h, mu, a.list))
cbind(a.list, rl.quantiles)

Compute the survival function of CUSUM run length

Description

Computation of the survival function of the Run Length (RL) for CUSUM control charts monitoring normal mean.

Usage

xcusum.sf(k, h, mu, n, hs=0, sided="one", r=40)

Arguments

Details

The survival function P(L>n) and derived from it also the cdf P(L<=n) and the pmf P(L=n) illustrate the distribution of the CUSUM run length. For large n the geometric tail could be exploited. That is, with reasonable large n the complete distribution is characterized. The algorithm is based on Waldmann's survival function iteration procedure.

Value

Returns a vector which resembles the survival function up to a certain point.

Author(s)

Sven Knoth

References

K.-H. Waldmann (1986), Bounds for the distribution of the run length of one-sided and two-sided CUSUM quality control schemes, Technometrics 28, 61-67.

See Also

xcusum.q for computation of CUSUM run length quantiles.

Examples

## Waldmann (1986), one-sided CUSUM, Table 2

k <- .5
h <- 3
mu <- 0 # corresponds to Waldmann's -0.5
SF <- xcusum.sf(k, h, 0, 1000)
plot(1:length(SF), SF, type="l", xlab="n", ylab="P(L>n)", ylim=c(0,1))
#

Compute steady-state ARLs of EWMA control charts

Description

Computation of the steady-state Average Run Length (ARL) for different types of EWMA control charts monitoring normal mean.

Usage

xewma.ad(l, c, mu1, mu0=0, zr=0, z0=0, sided="one", limits="fix",
steady.state.mode="conditional", r=40)

Arguments

Details

xewma.ad determines the steady-state Average Run Length (ARL) by numerically solving the related ARL integral equation by means of the Nystroem method based on Gauss-Legendre quadrature and using the power method for deriving the largest in magnitude eigenvalue and the related left eigenfunction.

Value

Returns a single value which resembles the steady-state ARL.

Author(s)

Sven Knoth

References

R. B. Crosier (1986), A new two-sided cumulative quality control scheme, Technometrics 28, 187-194.

S. V. Crowder (1987), A simple method for studying run-length distributions of exponentially weighted moving average charts, Technometrics 29, 401-407.

J. M. Lucas and M. S. Saccucci (1990), Exponentially weighted moving average control schemes: Properties and enhancements, Technometrics 32, 1-12.

See Also

xewma.arl for zero-state ARL computation and xcusum.ad for the steady-state ARL of CUSUM control charts.

Examples

## comparison of zero-state (= worst case ) and steady-state performance
## for two-sided EWMA control charts

l <- .1
c <- xewma.crit(l,500,sided="two")
mu <- c(0,.5,1,1.5,2)
arl <- sapply(mu,l=l,c=c,sided="two",xewma.arl)
ad <- sapply(mu,l=l,c=c,sided="two",xewma.ad)
round(cbind(mu,arl,ad),digits=2)

## Lucas/Saccucci (1990)
## two-sided EWMA

## with fixed limits
l1 <- .5
l2 <- .03
c1 <- 3.071
c2 <- 2.437
mu <- c(0,.25,.5,.75,1,1.5,2,2.5,3,3.5,4,5)
ad1 <- sapply(mu,l=l1,c=c1,sided="two",xewma.ad)
ad2 <- sapply(mu,l=l2,c=c2,sided="two",xewma.ad)
round(cbind(mu,ad1,ad2),digits=2)

## original results are (in Table 3)
## 0.00 499.   480.  
## 0.25 254.    74.1
## 0.50  88.4   28.6
## 0.75  35.7   17.3
## 1.00  17.3   12.5
## 1.50   6.44   8.00
## 2.00   3.58   5.95
## 2.50   2.47   4.78
## 3.00   1.91   4.02
## 3.50   1.58   3.49
## 4.00   1.36   3.09
## 5.00   1.10   2.55

Compute ARLs of EWMA control charts

Description

Computation of the (zero-state) Average Run Length (ARL) for different types of EWMA control charts monitoring normal mean.

Usage

xewma.arl(l,cE,mu,zr=0,hs=0,sided="one",limits="fix",q=1,
steady.state.mode="conditional",r=40)

Arguments

Details

In case of the EWMA chart with fixed control limits, xewma.arl determines the Average Run Length (ARL) by numerically solving the related ARL integral equation by means of the Nystroem method based on Gauss-Legendre quadrature. If limits is not "fix", then the method presented in Knoth (2003) is utilized. Note that for one-sided EWMA charts (sided="one"), only "vacl" and "stat" are deployed, while for two-sided ones (sided="two") also "fir", "both" (combination of "fir" and "vacl"), "Steiner" and "cfar" are implemented. For details see Knoth (2004).

Value

Except for the fixed limits EWMA charts it returns a single value which resembles the ARL. For fixed limits charts, it returns a vector of length q which resembles the ARL and the sequence of conditional expected delays for q=1 and q>1, respectively.

Author(s)

Sven Knoth

References

K.-H. Waldmann (1986), Bounds for the distribution of the run length of geometric moving average charts, Appl. Statist. 35, 151-158.

S. V. Crowder (1987), A simple method for studying run-length distributions of exponentially weighted moving average charts, Technometrics 29, 401-407.

J. M. Lucas and M. S. Saccucci (1990), Exponentially weighted moving average control schemes: Properties and enhancements, Technometrics 32, 1-12.

S. Chandrasekaran, J. R. English and R. L. Disney (1995), Modeling and analysis of EWMA control schemes with variance-adjusted control limits, IIE Transactions 277, 282-290.

T. R. Rhoads, D. C. Montgomery and C. M. Mastrangelo (1996), Fast initial response scheme for exponentially weighted moving average control chart, Quality Engineering 9, 317-327.

S. H. Steiner (1999), EWMA control charts with time-varying control limits and fast initial response, Journal of Quality Technology 31, 75-86.

S. Knoth (2003), EWMA schemes with non-homogeneous transition kernels, Sequential Analysis 22, 241-255.

S. Knoth (2004), Fast initial response features for EWMA Control Charts, Statistical Papers 46, 47-64.

See Also

xcusum.arl for zero-state ARL computation of CUSUM control charts and xewma.ad for the steady-state ARL.

Examples

## Waldmann (1986), one-sided EWMA
l <- .75
round(xewma.arl(l,2*sqrt((2-l)/l),0,zr=-4*sqrt((2-l)/l)),digits=1)
l <- .5
round(xewma.arl(l,2*sqrt((2-l)/l),0,zr=-4*sqrt((2-l)/l)),digits=1)
## original values are 209.3 and 3907.5 (in Table 2).

## Waldmann (1986), two-sided EWMA with fixed control limits
l <- .75
round(xewma.arl(l,2*sqrt((2-l)/l),0,sided="two"),digits=1)
l <- .5
round(xewma.arl(l,2*sqrt((2-l)/l),0,sided="two"),digits=1)
## original values are 104.0 and 1952 (in Table 1).

## Crowder (1987), two-sided EWMA with fixed control limits
l1 <- .5
l2 <- .05
cE <- 2
mu <- (0:16)/4
arl1 <- sapply(mu,l=l1,cE=cE,sided="two",xewma.arl)
arl2 <- sapply(mu,l=l2,cE=cE,sided="two",xewma.arl)
round(cbind(mu,arl1,arl2),digits=2)

## original results are (in Table 1)
## 0.00 26.45 127.53
## 0.25 20.12  43.94
## 0.50 11.89  18.97
## 0.75  7.29  11.64
## 1.00  4.91   8.38
## 1.25  3.95*  6.56
## 1.50  2.80   5.41
## 1.75  2.29   4.62
## 2.00  1.94   4.04
## 2.25  1.70   3.61
## 2.50  1.51   3.26
## 2.75  1.37   2.99
## 3.00  1.26   2.76
## 3.25  1.18   2.56
## 3.50  1.12   2.39
## 3.75  1.08   2.26
## 4.00  1.05   2.15  (* -- in Crowder (1987) typo!?)

## Lucas/Saccucci (1990)
## two-sided EWMA

## with fixed limits
l1 <- .5
l2 <- .03
c1 <- 3.071
c2 <- 2.437
mu <- c(0,.25,.5,.75,1,1.5,2,2.5,3,3.5,4,5)
arl1 <- sapply(mu,l=l1,cE=c1,sided="two",xewma.arl)
arl2 <- sapply(mu,l=l2,cE=c2,sided="two",xewma.arl)
round(cbind(mu,arl1,arl2),digits=2)

## original results are (in Table 3)
## 0.00 500.   500.
## 0.25 255.    76.7
## 0.50  88.8   29.3
## 0.75  35.9   17.6
## 1.00  17.5   12.6
## 1.50   6.53   8.07
## 2.00   3.63   5.99
## 2.50   2.50   4.80
## 3.00   1.93   4.03
## 3.50   1.58   3.49
## 4.00   1.34   3.11
## 5.00   1.07   2.55

## Not run: 
## with fir feature
l1 <- .5
l2 <- .03
c1 <- 3.071
c2 <- 2.437
hs1 <- c1/2
hs2 <- c2/2
mu <- c(0,.5,1,2,3,5)
arl1 <- sapply(mu,l=l1,cE=c1,hs=hs1,sided="two",limits="fir",xewma.arl)
arl2 <- sapply(mu,l=l2,cE=c2,hs=hs2,sided="two",limits="fir",xewma.arl)
round(cbind(mu,arl1,arl2),digits=2)

## original results are (in Table 5)
## 0.0 487.   406.
## 0.5  86.1   18.4
## 1.0  15.9    7.36
## 2.0   2.87   3.43
## 3.0   1.45   2.34
## 5.0   1.01   1.57

## Chandrasekaran, English, Disney (1995)
## two-sided EWMA with fixed and variance adjusted limits (vacl)

l1 <- .25
l2 <- .1
c1s <- 2.9985
c1n <- 3.0042
c2s <- 2.8159
c2n <- 2.8452
mu <- c(0,.25,.5,.75,1,2)
arl1s <- sapply(mu,l=l1,cE=c1s,sided="two",xewma.arl)
arl1n <- sapply(mu,l=l1,cE=c1n,sided="two",limits="vacl",xewma.arl)
arl2s <- sapply(mu,l=l2,cE=c2s,sided="two",xewma.arl)
arl2n <- sapply(mu,l=l2,cE=c2n,sided="two",limits="vacl",xewma.arl)
round(cbind(mu,arl1s,arl1n,arl2s,arl2n),digits=2)

## original results are (in Table 2)
## 0.00 500.   500.   500.   500.
## 0.25 170.09 167.54 105.90  96.6
## 0.50  48.14  45.65  31.08  24.35
## 0.75  20.02  19.72  15.71  10.74
## 1.00  11.07   9.37  10.23   6.35
## 2.00   3.59   2.64   4.32   2.73

## The results in Chandrasekaran, English, Disney (1995) are not
## that accurate. Let us consider the more appropriate comparison

c1s <- xewma.crit(l1,500,sided="two")
c1n <- xewma.crit(l1,500,sided="two",limits="vacl")
c2s <- xewma.crit(l2,500,sided="two")
c2n <- xewma.crit(l2,500,sided="two",limits="vacl")
mu <- c(0,.25,.5,.75,1,2)
arl1s <- sapply(mu,l=l1,cE=c1s,sided="two",xewma.arl)
arl1n <- sapply(mu,l=l1,cE=c1n,sided="two",limits="vacl",xewma.arl)
arl2s <- sapply(mu,l=l2,cE=c2s,sided="two",xewma.arl)
arl2n <- sapply(mu,l=l2,cE=c2n,sided="two",limits="vacl",xewma.arl)
round(cbind(mu,arl1s,arl1n,arl2s,arl2n),digits=2)

## which demonstrate the abilities of the variance-adjusted limits
## scheme more explicitely.

## Rhoads, Montgomery, Mastrangelo (1996)
## two-sided EWMA with fixed and variance adjusted limits (vacl),
## with fir and both features

l <- .03
cE <- 2.437
mu <- c(0,.5,1,1.5,2,3,4)
sl <- sqrt(l*(2-l))
arlfix  <- sapply(mu,l=l,cE=cE,sided="two",xewma.arl)
arlvacl <- sapply(mu,l=l,cE=cE,sided="two",limits="vacl",xewma.arl)
arlfir  <- sapply(mu,l=l,cE=cE,hs=c/2,sided="two",limits="fir",xewma.arl)
arlboth <- sapply(mu,l=l,cE=cE,hs=c/2*sl,sided="two",limits="both",xewma.arl)
round(cbind(mu,arlfix,arlvacl,arlfir,arlboth),digits=1)

## original results are (in Table 1)
## 0.0 477.3* 427.9* 383.4* 286.2*
## 0.5  29.7   20.0   18.6   12.8
## 1.0  12.5    6.5    7.4    3.6
## 1.5   8.1    3.3    4.6    1.9
## 2.0   6.0    2.2    3.4    1.4
## 3.0   4.0    1.3    2.4    1.0
## 4.0   3.1    1.1    1.9    1.0
## * -- the in-control values differ sustainably from the true values!

## Steiner (1999)
## two-sided EWMA control charts with various modifications

## fixed vs. variance adjusted limits

l <- .05
cE <- 3
mu <- c(0,.25,.5,.75,1,1.5,2,2.5,3,3.5,4)
arlfix <- sapply(mu,l=l,cE=cE,sided="two",xewma.arl)
arlvacl <- sapply(mu,l=l,cE=cE,sided="two",limits="vacl",xewma.arl)
round(cbind(mu,arlfix,arlvacl),digits=1)

## original results are (in Table 2)
## 0.00 1379.0   1353.0
## 0.25  135.0    127.0
## 0.50   37.4     32.5 
## 0.75   20.0     15.6
## 1.00   13.5      9.0
## 1.50    8.3      4.5
## 2.00    6.0      2.8
## 2.50    4.8      2.0
## 3.00    4.0      1.6
## 3.50    3.4      1.3
## 4.00    3.0      1.1

## fir, both, and Steiner's modification

l <- .03
cfir <- 2.44
cboth <- 2.54
cstein <- 2.55
hsfir <- cfir/2
hsboth <- cboth/2*sqrt(l*(2-l))
mu <- c(0,.5,1,1.5,2,3,4)
arlfir <- sapply(mu,l=l,cE=cfir,hs=hsfir,sided="two",limits="fir",xewma.arl)
arlboth <- sapply(mu,l=l,cE=cboth,hs=hsboth,sided="two",limits="both",xewma.arl)
arlstein <- sapply(mu,l=l,cE=cstein,sided="two",limits="Steiner",xewma.arl)
round(cbind(mu,arlfir,arlboth,arlstein),digits=1)

## original values are (in Table 5)
## 0.0 383.0   384.0   391.0
## 0.5  18.6    14.9    13.8
## 1.0   7.4     3.9     3.6
## 1.5   4.6     2.0     1.8
## 2.0   3.4     1.4     1.3
## 3.0   2.4     1.1     1.0
## 4.0   1.9     1.0     1.0

## SAS/QC manual 1999
## two-sided EWMA control charts with fixed limits

l <- .25
cE <- 3
mu <- 1
print(xewma.arl(l,cE,mu,sided="two"),digits=11)

# original value is 11.154267016.

## Some recent examples for one-sided EWMA charts
## with varying limits and in the so-called stationary mode

# 1. varying limits = "vacl"

lambda <- .1
L0 <- 500

## Monte Carlo results (10^9 replicates)
# mu    ARL      s.e.
# 0     500.00   0.0160
# 0.5   21.637   0.0006
# 1     6.7596   0.0001
# 1.5   3.5398   0.0001
# 2     2.3038   0.0000
# 2.5   1.7004   0.0000
# 3     1.3675   0.0000

zr <- -6
r <- 50
cE <- xewma.crit(lambda, L0, zr=zr, limits="vacl", r=r)
Mxewma.arl <- Vectorize(xewma.arl, "mu")
mus <- (0:6)/2
arls <- round(Mxewma.arl(lambda, cE, mus, zr=zr, limits="vacl", r=r), digits=4)
data.frame(mus, arls)

# 2. stationary mode, i. e. limits = "stat"

## Monte Carlo results (10^9 replicates)
# mu    ARL      s.e.
# 0     500.00   0.0159
# 0.5   22.313   0.0006
# 1     7.2920   0.0001
# 1.5   3.9064   0.0001
# 2     2.5131   0.0000
# 2.5   1.7983   0.0000
# 3     1.4029   0.0000

cE <- xewma.crit(lambda, L0, zr=zr, limits="stat", r=r)
arls <- round(Mxewma.arl(lambda, cE, mus, zr=zr, limits="stat", r=r), digits=4)
data.frame(mus, arls)

## End(Not run)

Compute ARL function of EWMA control charts

Description

Computation of the (zero-state) Average Run Length (ARL) function for different types of EWMA control charts monitoring normal mean.

Usage

xewma.arl.f(l,c,mu,zr=0,sided="one",limits="fix",r=40)

Arguments

Details

It is a convenience function to yield the ARL as function of the head start hs. For more details see xewma.arl.

Value

It returns a function of a single argument, hs=x which maps the head-start value hs to the ARL.

Author(s)

Sven Knoth

References

S. V. Crowder (1987), A simple method for studying run-length distributions of exponentially weighted moving average charts, Technometrics 29, 401-407.

See Also

xewma.arl for zero-state ARL for one specific head-start hs.

Examples

# will follow

Compute ARLs of EWMA control charts in case of estimated parameters

Description

Computation of the (zero-state) Average Run Length (ARL) for different types of EWMA control charts monitoring normal mean if the in-control mean, standard deviation, or both are estimated by a pre run.

Usage

xewma.arl.prerun(l, c, mu, zr=0, hs=0, sided="two", limits="fix", q=1,
size=100, df=NULL, estimated="mu", qm.mu=30, qm.sigma=30, truncate=1e-10)

xewma.crit.prerun(l, L0, mu, zr=0, hs=0, sided="two", limits="fix", size=100,
df=NULL, estimated="mu", qm.mu=30, qm.sigma=30, truncate=1e-10,
c.error=1e-12, L.error=1e-9, OUTPUT=FALSE)

Arguments

Details

Essentially, the ARL function xewma.arl is convoluted with the distribution of the sample mean, standard deviation or both. For details see Jones/Champ/Rigdon (2001) and Knoth (2014?).

Value

Returns a single value which resembles the ARL.

Author(s)

Sven Knoth

References

L. A. Jones, C. W. Champ, S. E. Rigdon (2001), The performance of exponentially weighted moving average charts with estimated parameters, Technometrics 43, 156-167.

S. Knoth (2003), EWMA schemes with non-homogeneous transition kernels, Sequential Analysis 22, 241-255.

S. Knoth (2004), Fast initial response features for EWMA Control Charts, Statistical Papers 46, 47-64.

S. Knoth (2014?), tbd, tbd, tbd-tbd.

See Also

xewma.arl for the usual zero-state ARL computation.

Examples

## Jones/Champ/Rigdon (2001)

c4m <- function(m, n) sqrt(2)*gamma( (m*(n-1)+1)/2 )/sqrt( m*(n-1) )/gamma( m*(n-1)/2 )

n <- 5 # sample size
m <- 20 # pre run with 20 samples of size n = 5
C4m <- c4m(m, n) # needed for bias correction

# Table 1, 3rd column
lambda <- 0.2
L <- 2.636

xewma.ARL <- Vectorize("xewma.arl", "mu")
xewma.ARL.prerun <- Vectorize("xewma.arl.prerun", "mu")

mu <- c(0, .25, .5, 1, 1.5, 2)
ARL <- round(xewma.ARL(lambda, L, mu, sided="two"), digits=2)
p.ARL <- round(xewma.ARL.prerun(lambda, L/C4m, mu, sided="two",
size=m, df=m*(n-1), estimated="both", qm.mu=70), digits=2)

# Monte-Carlo with 10^8 repetitions: 200.325 (0.020) and 144.458 (0.022)
cbind(mu, ARL, p.ARL)

## Not run: 
# Figure 5, subfigure r = 0.2
mu_ <- (0:85)/40
ARL_ <- round(xewma.ARL(lambda, L, mu_, sided="two"), digits=2)
p.ARL_ <- round(xewma.ARL.prerun(lambda, L/C4m, mu_, sided="two",
size=m, df=m*(n-1), estimated="both"), digits=2)

plot(mu_, ARL_, type="l", xlab=expression(delta), ylab="ARL", xlim=c(0,2))
abline(v=0, h=0, col="grey", lwd=.7)
points(mu, ARL, pch=5)
lines(mu_, p.ARL_, col="blue")
points(mu, p.ARL, pch=18, col ="blue")
legend("topright", c("Known", "Estimated"), col=c("black", "blue"),
lty=1, pch=c(5, 18))

## End(Not run)

Compute critical values of EWMA control charts

Description

Computation of the critical values (similar to alarm limits) for different types of EWMA control charts monitoring normal mean.

Usage

xewma.crit(l,L0,mu0=0,zr=0,hs=0,sided="one",limits="fix",r=40,c0=NULL,nmax=10000)

Arguments

Details

xewma.crit determines the critical values (similar to alarm limits) for given in-control ARL L0 by applying secant rule and using xewma.arl().

Value

Returns a single value which resembles the critical value c.

Author(s)

Sven Knoth

References

S. V. Crowder (1989), Design of exponentially weighted moving average schemes, Journal of Quality Technology 21, 155-162.

See Also

xewma.arl for zero-state ARL computation.

Examples

l <- .1
incontrolARL <- c(500,5000,50000)
sapply(incontrolARL,l=l,sided="two",xewma.crit,r=35) # accuracy with 35 nodes
sapply(incontrolARL,l=l,sided="two",xewma.crit)      # accuracy with 40 nodes
sapply(incontrolARL,l=l,sided="two",xewma.crit,r=50) # accuracy with 50 nodes

## Crowder (1989)
## two-sided EWMA control charts with fixed limits

l <- c(.05,.1,.15,.2,.25)
L0 <- 250
round(sapply(l,L0=L0,sided="two",xewma.crit),digits=2)

## original values are 2.32, 2.55, 2.65, 2.72, and 2.76.

Compute RL quantiles of EWMA control charts

Description

Computation of quantiles of the Run Length (RL) for EWMA control charts monitoring normal mean.

Usage

xewma.q(l, c, mu, alpha, zr=0, hs=0, sided="two", limits="fix", q=1, r=40)

xewma.q.crit(l, L0, mu, alpha, zr=0, hs=0, sided="two", limits="fix", r=40,
c.error=1e-12, a.error=1e-9, OUTPUT=FALSE)

Arguments

Details

Instead of the popular ARL (Average Run Length) quantiles of the EWMA stopping time (Run Length) are determined. The algorithm is based on Waldmann's survival function iteration procedure. If limits is not "fix", then the method presented in Knoth (2003) is utilized. Note that for one-sided EWMA charts (sided="one"), only "vacl" and "stat" are deployed, while for two-sided ones (sided="two") also "fir", "both" (combination of "fir" and "vacl"), and "Steiner" are implemented. For details see Knoth (2004).

Value

Returns a single value which resembles the RL quantile of order q.

Author(s)

Sven Knoth

References

F. F. Gan (1993), An optimal design of EWMA control charts based on the median run length, J. Stat. Comput. Simulation 45, 169-184.

S. Knoth (2003), EWMA schemes with non-homogeneous transition kernels, Sequential Analysis 22, 241-255.

S. Knoth (2004), Fast initial response features for EWMA Control Charts, Statistical Papers 46, 47-64.

S. Knoth (2015), Run length quantiles of EWMA control charts monitoring normal mean or/and variance, International Journal of Production Research 53, 4629-4647.

K.-H. Waldmann (1986), Bounds for the distribution of the run length of geometric moving average charts, Appl. Statist. 35, 151-158.

See Also

xewma.arl for zero-state ARL computation of EWMA control charts.

Examples

## Gan (1993), two-sided EWMA with fixed control limits
## some values of his Table 1 -- any median RL should be 500
XEWMA.Q <- Vectorize("xewma.q", c("l", "c"))
G.lambda <- c(.05, .1, .15, .2, .25)
G.h <- c(.441, .675, .863, 1.027, 1.177)
MEDIAN <- ceiling(XEWMA.Q(G.lambda, G.h/sqrt(G.lambda/(2-G.lambda)),
0, .5, sided="two"))
print(cbind(G.lambda, MEDIAN))

## increase accuracy of thresholds

# (i) calculate threshold for given in-control median value by
#     deplyoing secant rule
XEWMA.q.crit <- Vectorize("xewma.q.crit", "l")

# (ii) re-calculate the thresholds and remove the standardization step
L0 <- 500
G.h.new <- XEWMA.q.crit(G.lambda, L0, 0, .5, sided="two")
G.h.new <- round(G.h.new * sqrt(G.lambda/(2-G.lambda)), digits=5)

# (iii) compare Gan's original values and the new ones with 5 digits
print(cbind(G.lambda, G.h.new, G.h))

# (iv) calculate the new medians
MEDIAN <- ceiling(XEWMA.Q(G.lambda, G.h.new/sqrt(G.lambda/(2-G.lambda)),
0, .5, sided="two"))
print(cbind(G.lambda, MEDIAN))

Compute RL quantiles of EWMA control charts in case of estimated parameters

Description

Computation of quantiles of the Run Length (RL) for EWMA control charts monitoring normal mean if the in-control mean, standard deviation, or both are estimated by a pre run.

Usage

xewma.q.prerun(l, c, mu, p, zr=0, hs=0, sided="two", limits="fix", q=1, size=100,
df=NULL, estimated="mu", qm.mu=30, qm.sigma=30, truncate=1e-10, bound=1e-10)

xewma.q.crit.prerun(l, L0, mu, p, zr=0, hs=0, sided="two", limits="fix", size=100,
df=NULL, estimated="mu", qm.mu=30, qm.sigma=30, truncate=1e-10, bound=1e-10,
c.error=1e-10, p.error=1e-9, OUTPUT=FALSE)

Arguments

Details

Essentially, the ARL function xewma.q is convoluted with the distribution of the sample mean, standard deviation or both. For details see Jones/Champ/Rigdon (2001) and Knoth (2014?).

Value

Returns a single value which resembles the RL quantile of order q.

Author(s)

Sven Knoth

References

L. A. Jones, C. W. Champ, S. E. Rigdon (2001), The performance of exponentially weighted moving average charts with estimated parameters, Technometrics 43, 156-167.

S. Knoth (2003), EWMA schemes with non-homogeneous transition kernels, Sequential Analysis 22, 241-255.

S. Knoth (2004), Fast initial response features for EWMA Control Charts, Statistical Papers 46, 47-64.

S. Knoth (2014?), tbd, tbd, tbd-tbd.

K.-H. Waldmann (1986), Bounds for the distribution of the run length of geometric moving average charts, Appl. Statist. 35, 151-158.

See Also

xewma.q for the usual RL quantiles computation of EWMA control charts.

Examples

## Jones/Champ/Rigdon (2001)

c4m <- function(m, n) sqrt(2)*gamma( (m*(n-1)+1)/2 )/sqrt( m*(n-1) )/gamma( m*(n-1)/2 )

n <- 5 # sample size
m <- 20 # pre run with 20 samples of size n = 5
C4m <- c4m(m, n) # needed for bias correction

# Table 1, 3rd column
lambda <- 0.2
L <- 2.636

xewma.Q <- Vectorize("xewma.q", "mu")
xewma.Q.prerun <- Vectorize("xewma.q.prerun", "mu")

mu <- c(0, .25, .5, 1, 1.5, 2)
Q1  <- ceiling(xewma.Q(lambda, L, mu, 0.1, sided="two"))
Q2  <- ceiling(xewma.Q(lambda, L, mu, 0.5, sided="two"))
Q3  <- ceiling(xewma.Q(lambda, L, mu, 0.9, sided="two"))

cbind(mu, Q1, Q2, Q3)

## Not run: 
p.Q1 <- xewma.Q.prerun(lambda, L/C4m, mu, 0.1, sided="two", 
size=m, df=m*(n-1), estimated="both")
p.Q2 <- xewma.Q.prerun(lambda, L/C4m, mu, 0.5, sided="two",
size=m, df=m*(n-1), estimated="both")
p.Q3 <- xewma.Q.prerun(lambda, L/C4m, mu, 0.9, sided="two",
size=m, df=m*(n-1), estimated="both")

cbind(mu, p.Q1, p.Q2, p.Q3)

## End(Not run)

## original values are
#    mu Q1  Q2  Q3 p.Q1 p.Q2 p.Q3
#  0.00 25 140 456   13   73  345
#  0.25 12  56 174    9   46  253
#  0.50  7  20  56    6   20  101
#  1.00  4   7  15    3    7   18
#  1.50  3   4   7    2    4    8
#  2.00  2   3   5    2    3    5

Compute the survival function of EWMA run length

Description

Computation of the survival function of the Run Length (RL) for EWMA control charts monitoring normal mean.

Usage

xewma.sf(l, c, mu, n, zr=0, hs=0, sided="one", limits="fix", q=1, r=40)

Arguments

Details

The survival function P(L>n) and derived from it also the cdf P(L<=n) and the pmf P(L=n) illustrate the distribution of the EWMA run length. For large n the geometric tail could be exploited. That is, with reasonable large n the complete distribution is characterized. The algorithm is based on Waldmann's survival function iteration procedure. For varying limits and for change points after 1 the algorithm from Knoth (2004) is applied. Note that for one-sided EWMA charts (sided="one"), only "vacl" and "stat" are deployed, while for two-sided ones (sided="two") also "fir", "both" (combination of "fir" and "vacl"), and "Steiner" are implemented. For details see Knoth (2004).

Value

Returns a vector which resembles the survival function up to a certain point.

Author(s)

Sven Knoth

References

F. F. Gan (1993), An optimal design of EWMA control charts based on the median run length, J. Stat. Comput. Simulation 45, 169-184.

S. Knoth (2003), EWMA schemes with non-homogeneous transition kernels, Sequential Analysis 22, 241-255.

S. Knoth (2004), Fast initial response features for EWMA Control Charts, Statistical Papers 46, 47-64.

K.-H. Waldmann (1986), Bounds for the distribution of the run length of geometric moving average charts, Appl. Statist. 35, 151-158.

See Also

xewma.arl for zero-state ARL computation of EWMA control charts.

Examples

## Gan (1993), two-sided EWMA with fixed control limits
## some values of his Table 1 -- any median RL should be 500

G.lambda <- c(.05, .1, .15, .2, .25)
G.h <- c(.441, .675, .863, 1.027, 1.177)/sqrt(G.lambda/(2-G.lambda))

for ( i in 1:length(G.lambda) ) {
  SF <- xewma.sf(G.lambda[i], G.h[i], 0, 1000)
  if (i==1) plot(1:length(SF), SF, type="l", xlab="n", ylab="P(L>n)")
  else lines(1:length(SF), SF, col=i)
}

Compute the survival function of EWMA run length in case of estimated parameters

Description

Computation of the survival function of the Run Length (RL) for EWMA control charts monitoring normal mean if the in-control mean, standard deviation, or both are estimated by a pre run.

Usage

xewma.sf.prerun(l, c, mu, n, zr=0, hs=0, sided="one", limits="fix", q=1,
size=100, df=NULL, estimated="mu", qm.mu=30, qm.sigma=30,
truncate=1e-10, tail_approx=TRUE, bound=1e-10)

Arguments

Details

The survival function P(L>n) and derived from it also the cdf P(L<=n) and the pmf P(L=n) illustrate the distribution of the EWMA run length...

Value

Returns a vector which resembles the survival function up to a certain point.

Author(s)

Sven Knoth

References

F. F. Gan (1993), An optimal design of EWMA control charts based on the median run length, J. Stat. Comput. Simulation 45, 169-184.

S. Knoth (2003), EWMA schemes with non-homogeneous transition kernels, Sequential Analysis 22, 241-255.

S. Knoth (2004), Fast initial response features for EWMA Control Charts, Statistical Papers 46, 47-64.

L. A. Jones, C. W. Champ, S. E. Rigdon (2001), The performance of exponentially weighted moving average charts with estimated parameters, Technometrics 43, 156-167.

K.-H. Waldmann (1986), Bounds for the distribution of the run length of geometric moving average charts, Appl. Statist. 35, 151-158.

See Also

xewma.sf for the RL survival function of EWMA control charts w/o pre run uncertainty.

Examples

## Jones/Champ/Rigdon (2001)

c4m <- function(m, n) sqrt(2)*gamma( (m*(n-1)+1)/2 )/sqrt( m*(n-1) )/gamma( m*(n-1)/2 )

n <- 5 # sample size

# Figure 6, subfigure r=0.1
lambda <- 0.1
L <- 2.454

CDF0 <- 1 - xewma.sf(lambda, L, 0, 600, sided="two")
m <- 10 # pre run size
CDF1 <- 1 - xewma.sf.prerun(lambda, L/c4m(m,n), 0, 600, sided="two",
size=m, df=m*(n-1), estimated="both")
m <- 20
CDF2 <- 1 - xewma.sf.prerun(lambda, L/c4m(m,n), 0, 600, sided="two",
size=m, df=m*(n-1), estimated="both")
m <- 50
CDF3 <- 1 - xewma.sf.prerun(lambda, L/c4m(m,n), 0, 600, sided="two",
size=m, df=m*(n-1), estimated="both")

plot(CDF0, type="l", xlab="t", ylab=expression(P(T<=t)), xlim=c(0,500), ylim=c(0,1))
abline(v=0, h=c(0,1), col="grey", lwd=.7)
points((1:5)*100, CDF0[(1:5)*100], pch=18)
lines(CDF1, col="blue")
points((1:5)*100, CDF1[(1:5)*100], pch=2, col="blue")
lines(CDF2, col="red")
points((1:5)*100, CDF2[(1:5)*100], pch=16, col="red")
lines(CDF3, col="green")
points((1:5)*100, CDF3[(1:5)*100], pch=5, col="green")

legend("bottomright", c("Known", "m=10, n=5", "m=20, n=5", "m=50, n=5"),
       col=c("black", "blue", "red", "green"), pch=c(18, 2, 16, 5), lty=1)

Compute steady-state ARLs of Shiryaev-Roberts schemes

Description

Computation of the steady-state Average Run Length (ARL) for Shiryaev-Roberts schemes monitoring normal mean.

Usage

xgrsr.ad(k, g, mu1, mu0 = 0, zr = 0, sided = "one", MPT = FALSE, r = 30)

Arguments

Details

xgrsr.ad determines the steady-state Average Run Length (ARL) by numerically solving the related ARL integral equation by means of the Nystroem method based on Gauss-Legendre quadrature.

Value

Returns a single value which resembles the steady-state ARL.

Author(s)

Sven Knoth

References

S. Knoth (2006), The art of evaluating monitoring schemes – how to measure the performance of control charts? S. Lenz, H. & Wilrich, P. (ed.), Frontiers in Statistical Quality Control 8, Physica Verlag, Heidelberg, Germany, 74-99.

G. Moustakides, A. Polunchenko, A. Tartakovsky (2009), Numerical comparison of CUSUM and Shiryaev-Roberts procedures for detectin changes in distributions, Communications in Statistics: Theory and Methods 38, 3225-3239.

See Also

xewma.arl and xcusum-arl for zero-state ARL computation of EWMA and CUSUM control charts, respectively, and xgrsr.arl for the zero-state ARL.

Examples

## Small study to identify appropriate reflection border to mimic unreflected schemes
k <- .5
g <- log(390)
zrs <- -(0:10)
ZRxgrsr.ad <- Vectorize(xgrsr.ad, "zr")
ads <- ZRxgrsr.ad(k, g, 0, zr=zrs)
data.frame(zrs, ads)

## Table 2 from Knoth (2006)
## original values are
#  mu   arl
#  0    689
#  0.5  30
#  1    8.9
#  1.5  5.1
#  2    3.6
#  2.5  2.8
#  3    2.4
#
k <- .5
g <- log(390)
zr <- -5 # see first example
mus <- (0:6)/2
Mxgrsr.ad <- Vectorize(xgrsr.ad, "mu1")
ads <- round(Mxgrsr.ad(k, g, mus, zr=zr), digits=1)
data.frame(mus, ads)

## Table 4 from Moustakides et al. (2009)
## original values are
# gamma  A        STADD/steady-state ARL
# 50     28.02    4.37
# 100    56.04    5.46
# 500    280.19   8.33
# 1000   560.37   9.64
# 5000   2801.75  12.79
# 10000  5603.7   14.17
Gxgrsr.ad  <- Vectorize("xgrsr.ad", "g")
As <- c(28.02, 56.04, 280.19, 560.37, 2801.75, 5603.7)
gs <- log(As)
theta <- 1
zr <- -6
ads <- round(Gxgrsr.ad(theta/2, gs, theta, zr=zr, r=100), digits=2)
data.frame(As, ads)

Compute (zero-state) ARLs of Shiryaev-Roberts schemes

Description

Computation of the (zero-state) Average Run Length (ARL) for Shiryaev-Roberts schemes monitoring normal mean.

Usage

xgrsr.arl(k, g, mu, zr = 0, hs=NULL, sided = "one", q = 1, MPT = FALSE, r = 30)

Arguments

Details

xgrsr.arl determines the Average Run Length (ARL) by numerically solving the related ARL integral equation by means of the Nystroem method based on Gauss-Legendre quadrature.

Value

Returns a vector of length q which resembles the ARL and the sequence of conditional expected delays for q=1 and q>1, respectively.

Author(s)

Sven Knoth

References

S. Knoth (2006), The art of evaluating monitoring schemes – how to measure the performance of control charts? S. Lenz, H. & Wilrich, P. (ed.), Frontiers in Statistical Quality Control 8, Physica Verlag, Heidelberg, Germany, 74-99.

G. Moustakides, A. Polunchenko, A. Tartakovsky (2009), Numerical comparison of CUSUM and Shiryaev-Roberts procedures for detecting changes in distributions, Communications in Statistics: Theory and Methods 38, 3225-3239.

See Also

xewma.arl and xcusum-arl for zero-state ARL computation of EWMA and CUSUM control charts, respectively, and xgrsr.ad for the steady-state ARL.

Examples

## Small study to identify appropriate reflection border to mimic unreflected schemes
k <- .5
g <- log(390)
zrs <- -(0:10)
ZRxgrsr.arl <- Vectorize(xgrsr.arl, "zr")
arls <- ZRxgrsr.arl(k, g, 0, zr=zrs)
data.frame(zrs, arls)

## Table 2 from Knoth (2006)
## original values are
#  mu   arl
#  0    697
#  0.5  33
#  1    10.4
#  1.5  6.2
#  2    4.4
#  2.5  3.5
#  3    2.9
#
k <- .5
g <- log(390)
zr <- -5 # see first example
mus <- (0:6)/2
Mxgrsr.arl <- Vectorize(xgrsr.arl, "mu")
arls <- round(Mxgrsr.arl(k, g, mus, zr=zr), digits=1)
data.frame(mus, arls)

XGRSR.arl  <- Vectorize("xgrsr.arl", "g")
zr <- -6

## Table 2 from Moustakides et al. (2009)
## original values are
# gamma   A     ARL/E_infty(L) SADD/E_1(L)
#   50   47.17      50.29        41.40
#  100   94.34     100.28        72.32
#  500  471.70     500.28       209.44
# 1000  943.41    1000.28       298.50
# 5000 4717.04    5000.24       557.87
#10000 9434.08   10000.17       684.17

theta <- .1
As2 <- c(47.17, 94.34, 471.7, 943.41, 4717.04, 9434.08)
gs2 <- log(As2)
arls0 <- round(XGRSR.arl(theta/2, gs2, 0, zr=-5, r=300, MPT=TRUE), digits=2)
arls1 <- round(XGRSR.arl(theta/2, gs2, theta, zr=-5, r=300, MPT=TRUE), digits=2)
data.frame(As2, arls0, arls1)

## Table 3 from Moustakides et al. (2009)
## original values are
# gamma   A     ARL/E_infty(L) SADD/E_1(L)
#   50   37.38      49.45        12.30
#  100   74.76      99.45        16.60
#  500  373.81     499.45        28.05
# 1000  747.62     999.45        33.33
# 5000 3738.08    4999.45        45.96
#10000 7476.15    9999.24        51.49

theta <- .5
As3 <- c(37.38, 74.76, 373.81, 747.62, 3738.08, 7476.15)
gs3 <- log(As3)
arls0 <- round(XGRSR.arl(theta/2, gs3, 0, zr=-5, r=70, MPT=TRUE), digits=2)
arls1 <- round(XGRSR.arl(theta/2, gs3, theta, zr=-5, r=70, MPT=TRUE), digits=2)
data.frame(As3, arls0, arls1)

## Table 4 from Moustakides et al. (2009)
## original values are
# gamma   A     ARL/E_infty(L) SADD/E_1(L)
#   50   28.02      49.78         4.98
#  100   56.04      99.79         6.22
#  500  280.19     499.79         9.30
# 1000  560.37     999.79        10.66
# 5000 2801.85    5000.93        13.86
#10000 5603.70    9999.87        15.24

theta <- 1
As4 <- c(28.02, 56.04, 280.19, 560.37, 2801.85, 5603.7)
gs4 <- log(As4)
arls0 <- round(XGRSR.arl(theta/2, gs4, 0, zr=-6, r=40, MPT=TRUE), digits=2)
arls1 <- round(XGRSR.arl(theta/2, gs4, theta, zr=-6, r=40, MPT=TRUE), digits=2)
data.frame(As4, arls0, arls1)

Compute alarm thresholds for Shiryaev-Roberts schemes

Description

Computation of the alarm thresholds (alarm limits) for Shiryaev-Roberts schemes monitoring normal mean.

Usage

xgrsr.crit(k, L0, mu0 = 0, zr = 0, hs = NULL, sided = "one", MPT = FALSE, r = 30)

Arguments

Details

xgrsr.crit determines the alarm threshold (alarm limit) for given in-control ARL L0 by applying secant rule and using xgrsr.arl().

Value

Returns a single value which resembles the alarm limit g.

Author(s)

Sven Knoth

References

G. Moustakides, A. Polunchenko, A. Tartakovsky (2009), Numerical comparison of CUSUM and Shiryaev-Roberts procedures for detecting changes in distributions, Communications in Statistics: Theory and Methods 38, 3225-3239.r.

See Also

xgrsr.arl for zero-state ARL computation.

Examples

## Table 4 from Moustakides et al. (2009)
## original values are
# gamma/L0  A/exp(g)
# 50        28.02
# 100       56.04
# 500       280.19
# 1000      560.37
# 5000      2801.75
# 10000     5603.7
theta <- 1
zr <- -6
r <- 100
Lxgrsr.crit  <- Vectorize("xgrsr.crit", "L0")
L0s <- c(50, 100, 500, 1000, 5000, 10000)
gs <- Lxgrsr.crit(theta/2, L0s, zr=zr, r=r)
data.frame(L0s, gs, A=round(exp(gs), digits=2))

Compute ARLs of simultaneous EWMA control charts (mean and variance charts)

Description

Computation of the (zero-state) Average Run Length (ARL) for different types of simultaneous EWMA control charts (based on the sample mean and the sample variance S^2) monitoring normal mean and variance.

Usage

xsewma.arl(lx, cx, ls, csu, df, mu, sigma, hsx=0, Nx=40, csl=0,
hss=1, Ns=40, s2.on=TRUE, sided="upper", qm=30)

Arguments

Details

xsewma.arl determines the Average Run Length (ARL) by an extension of Gan's (derived from ideas already published by Waldmann) algorithm. The variance EWMA part is treated similarly to the ARL calculation method deployed for the single variance EWMA charts in Knoth (2005), that is, by means of collocation (Chebyshev polynomials). For more details see Knoth (2007).

Value

Returns a single value which resembles the ARL.

Author(s)

Sven Knoth

References

K. Waldmann (1986), Bounds for the distribution of the run length of geometric moving average charts, J. R. Stat. Soc., Ser. C, Appl. Stat. 35, 151-158.

F. F. Gan (1995), Joint monitoring of process mean and variance using exponentially weighted moving average control charts, Technometrics 37, 446-453.

S. Knoth (2005), Accurate ARL computation for EWMA-S^2 control charts, Statistics and Computing 15, 341-352.

S. Knoth (2007), Accurate ARL calculation for EWMA control charts monitoring simultaneously normal mean and variance, Sequential Analysis 26, 251-264.

See Also

xewma.arl and sewma.arl for zero-state ARL computation of single mean and variance EWMA control charts, respectively.

Examples

## Knoth (2007)
## collocation results in Table 1
## Monte Carlo with 10^9 replicates: 252.307 +/- 0.0078

# process parameters
mu <- 0
sigma <- 1
# subgroup size n=5, df=n-1
df  <- 4
# lambda of mean chart
lx  <- .134
# c_mu^* = .345476571 = cx/sqrt(n) * sqrt(lx/(2-lx)
cx  <- .345476571*sqrt(df+1)/sqrt(lx/(2-lx))
# lambda of variance chart
ls  <- .1
# c_sigma = .477977
csu <- 1 + .477977
# matrix dimensions for mean and variance part
Nx  <- 25
Ns  <- 25
# mode of variance chart
SIDED <- "upper"

arl <- xsewma.arl(lx, cx, ls, csu, df, mu, sigma, Nx=Nx, Ns=Ns, sided=SIDED)
arl

Compute critical values of simultaneous EWMA control charts (mean and variance charts)

Description

Computation of the critical values (similar to alarm limits) for different types of simultaneous EWMA control charts (based on the sample mean and the sample variance S^2) monitoring normal mean and variance.

Usage

xsewma.crit(lx, ls, L0, df, mu0=0, sigma0=1, cu=NULL, hsx=0,
hss=1, s2.on=TRUE, sided="upper", mode="fixed", Nx=30, Ns=40, qm=30)

Arguments

Details

xsewma.crit determines the critical values (similar to alarm limits) for given in-control ARL L0 by applying secant rule and using xsewma.arl(). In case of sided="two" and mode="unbiased" a two-dimensional secant rule is applied that also ensures that the maximum of the ARL function for given standard deviation is attained at sigma0. See Knoth (2007) for details and application.

Value

Returns the critical value of the two-sided mean EWMA chart and the lower and upper controls limit cl and cu of the variance EWMA chart.

Author(s)

Sven Knoth

References

S. Knoth (2007), Accurate ARL calculation for EWMA control charts monitoring simultaneously normal mean and variance, Sequential Analysis 26, 251-264.

See Also

xsewma.arl for calculation of ARL of simultaneous EWMA charts.

Examples

## Knoth (2007)
## results in Table 2

# subgroup size n=5, df=n-1
df  <- 4
# lambda of mean chart
lx  <- .134
# lambda of variance chart
ls  <- .1
# in-control ARL 
L0 <- 252.3
# matrix dimensions for mean and variance part
Nx  <- 25
Ns  <- 25
# mode of variance chart
SIDED <- "upper"

crit <- xsewma.crit(lx, ls, L0, df, sided=SIDED, Nx=Nx, Ns=Ns)
crit

## output as used in Knoth (2007)
crit["cx"]/sqrt(df+1)*sqrt(lx/(2-lx))
crit["cu"] - 1

Compute critical values of simultaneous EWMA control charts (mean and variance charts) for given RL quantile

Description

Computation of the critical values (similar to alarm limits) for different types of simultaneous EWMA control charts (based on the sample mean and the sample variance S^2) monitoring normal mean and variance.

Usage

xsewma.q(lx, cx, ls, csu, df, alpha, mu, sigma, hsx=0,
Nx=40, csl=0, hss=1, Ns=40, sided="upper", qm=30)

xsewma.q.crit(lx, ls, L0, alpha, df, mu0=0, sigma0=1, csu=NULL,
hsx=0, hss=1, sided="upper", mode="fixed", Nx=20, Ns=40, qm=30,
c.error=1e-12, a.error=1e-9)

Arguments

Details

Instead of the popular ARL (Average Run Length) quantiles of the EWMA stopping time (Run Length) are determined. The algorithm is based on Waldmann's survival function iteration procedure and on Knoth (2007). xsewma.q.crit determines the critical values (similar to alarm limits) for given in-control RL quantile L0 at level alpha by applying secant rule and using xsewma.sf(). In case of sided="two" and mode="unbiased" a two-dimensional secant rule is applied that also ensures that the maximum of the RL cdf for given standard deviation is attained at sigma0.

Value

Returns a single value which resembles the RL quantile of order alpha and the critical value of the two-sided mean EWMA chart and the lower and upper controls limit csl and csu of the variance EWMA chart, respectively.

Author(s)

Sven Knoth

References

S. Knoth (2007), Accurate ARL calculation for EWMA control charts monitoring simultaneously normal mean and variance, Sequential Analysis 26, 251-264.

See Also

xsewma.arl for calculation of ARL of simultaneous EWMA charts and xsewma.sf for the RL survival function.

Examples

## will follow

Compute the survival function of simultaneous EWMA control charts (mean and variance charts)

Description

Computation of the survival function of the Run Length (RL) for EWMA control charts monitoring simultaneously normal mean and variance.

Usage

xsewma.sf(n, lx, cx, ls, csu, df, mu, sigma, hsx=0, Nx=40,
csl=0, hss=1, Ns=40, sided="upper", qm=30) 

Arguments

Details

The survival function P(L>n) and derived from it also the cdf P(L<=n) and the pmf P(L=n) illustrate the distribution of the EWMA run length. For large n the geometric tail could be exploited. That is, with reasonable large n the complete distribution is characterized. The algorithm is based on Waldmann's survival function iteration procedure and on results in Knoth (2007).

Value

Returns a vector which resembles the survival function up to a certain point.

Author(s)

Sven Knoth

References

S. Knoth (2007), Accurate ARL calculation for EWMA control charts monitoring simultaneously normal mean and variance, Sequential Analysis 26, 251-264.

K.-H. Waldmann (1986), Bounds for the distribution of the run length of geometric moving average charts, Appl. Statist. 35, 151-158.

See Also

xsewma.arl for zero-state ARL computation of simultaneous EWMA control charts.

Examples

## will follow

Compute ARLs of modified Shewhart control charts for AR(1) data

Description

Computation of the (zero-state) Average Run Length (ARL) for modified Shewhart charts deployed to the original AR(1) data.

Usage

xshewhart.ar1.arl(alpha, cS, delta=0, N1=50, N2=30)

Arguments

Details

Following the idea of Schmid (1995), 1- alpha times the data turns out to be an EWMA smoothing of the underlying AR(1) residuals. Hence, by combining the solution of the EWMA ARL integral equation and the stationary distribution of the AR(1) data (normal distribution is assumed) one gets easily the overall ARL.

Value

It returns a single value resembling the zero-state ARL of a modified Shewhart chart.

Author(s)

Sven Knoth

References

S. Knoth, W. Schmid (2004). Control charts for time series: A review. In Frontiers in Statistical Quality Control 7, edited by H.-J. Lenz, P.-T. Wilrich, 210-236, Physica-Verlag.

H. Kramer, W. Schmid (2000). The influence of parameter estimation on the ARL of Shewhart type charts for time series. Statistical Papers 41(2), 173-196.

W. Schmid (1995). On the run length of a Shewhart chart for correlated data. Statistical Papers 36(1), 111-130.

See Also

xewma.arl for zero-state ARL computation of EWMA control charts.

Examples

## Table 1 in Kramer/Schmid (2000)

cS <- 3.09023
a  <- seq(0, 4, by=.5)
row1 <- row2 <- row3 <- NULL
for ( i in 1:length(a) ) {
  row1 <- c(row1, round(xshewhart.ar1.arl( 0.4, cS, delta=a[i]), digits=2))
  row2 <- c(row2, round(xshewhart.ar1.arl( 0.2, cS, delta=a[i]), digits=2))
  row3 <- c(row3, round(xshewhart.ar1.arl(-0.2, cS, delta=a[i]), digits=2))
}

results <- rbind(row1, row2, row3)
results

# original values are
# row1 515.44 215.48 61.85 21.63 9.19 4.58 2.61 1.71 1.29
# row2 502.56 204.97 56.72 19.13 7.95 3.97 2.33 1.59 1.25
# row3 502.56 201.41 54.05 17.42 6.89 3.37 2.03 1.46 1.20

Compute ARLs of Shewhart control charts with and without runs rules

Description

Computation of the (zero-state and steady-state) Average Run Length (ARL) for Shewhart control charts with and without runs rules monitoring normal mean.

Usage

xshewhartrunsrules.arl(mu, c = 1, type = "12")

xshewhartrunsrules.crit(L0, mu = 0, type = "12")

xshewhartrunsrules.ad(mu1, mu0 = 0, c = 1, type = "12")

xshewhartrunsrules.matrix(mu, c = 1, type = "12")

Arguments

Details

xshewhartrunsrules.arl determines the zero-state Average Run Length (ARL) based on the Markov chain approach given in Champ and Woodall (1987). xshewhartrunsrules.matrix provides the corresponding transition matrix that is also used in xDshewhartrunsrules.arl (ARL for control charting drift). xshewhartrunsrules.crit allows to find the normalization constant c to attain a fixed in-control ARL. Typically this is needed to calibrate the chart. With xshewhartrunsrules.ad the steady-state ARL is calculated. With the argument type certain runs rules could be set. The following list gives an overview.

"1"

The classical Shewhart chart with +/- 3 c sigma control limits (c is typically equal to 1 and can be changed by the argument c).

"12"

The classic and the following runs rule: 2 of 3 are beyond +/- 2 c sigma on the same side of the chart.

"13"

The classic and the following runs rule: 4 of 5 are beyond +/- 1 c sigma on the same side of the chart.

"14"

The classic and the following runs rule: 8 of 8 are on the same side of the chart (referring to the center line).

Value

Returns a single value which resembles the zero-state or steady-state ARL. xshewhartrunsrules.matrix returns a matrix.

Author(s)

Sven Knoth

References

C. W. Champ and W. H. Woodall (1987), Exact results for Shewhart control charts with supplementary runs rules, Technometrics 29, 393-399.

See Also

xDshewhartrunsrules.arl for zero-state ARL of Shewhart control charts with or without runs rules under drift.

Examples

## Champ/Woodall (1987)
## Table 1
mus <- (0:15)/5
Mxshewhartrunsrules.arl <- Vectorize(xshewhartrunsrules.arl, "mu")
# standard (1 of 1 beyond 3 sigma) Shewhart chart without runs rules
C1 <- round(Mxshewhartrunsrules.arl(mus, type="1"), digits=2)
# standard + runs rule: 2 of 3 beyond 2 sigma on the same side
C12 <- round(Mxshewhartrunsrules.arl(mus, type="12"), digits=2)
# standard + runs rule: 4 of 5 beyond 1 sigma on the same side
C13 <- round(Mxshewhartrunsrules.arl(mus, type="13"), digits=2)
# standard + runs rule: 8 of 8 on the same side of the center line
C14 <- round(Mxshewhartrunsrules.arl(mus, type="14"), digits=2)

## original results are
#  mus     C1    C12    C13    C14                                    
#  0.0 370.40 225.44 166.05 152.73                                    
#  0.2 308.43 177.56 120.70 110.52                                    
#  0.4 200.08 104.46  63.88  59.76                                    
#  0.6 119.67  57.92  33.99  33.64                                    
#  0.8  71.55  33.12  19.78  21.07                                    
#  1.0  43.89  20.01  12.66  14.58                                    
#  1.2  27.82  12.81   8.84  10.90                                    
#  1.4  18.25   8.69   6.62   8.60                                    
#  1.6  12.38   6.21   5.24   7.03                                    
#  1.8   8.69   4.66   4.33   5.85                                    
#  2.0   6.30   3.65   3.68   4.89                                    
#  2.2   4.72   2.96   3.18   4.08                                    
#  2.4   3.65   2.48   2.78   3.38                                    
#  2.6   2.90   2.13   2.43   2.81                                    
#  2.8   2.38   1.87   2.14   2.35                                    
#  3.0   2.00   1.68   1.89   1.99

data.frame(mus, C1, C12, C13, C14)

## plus calibration, i. e. L0=250 (the maximal value for "14" is 255!
L0 <- 250
c1  <- xshewhartrunsrules.crit(L0, type = "1")
c12 <- xshewhartrunsrules.crit(L0, type = "12")
c13 <- xshewhartrunsrules.crit(L0, type = "13")
c14 <- xshewhartrunsrules.crit(L0, type = "14")
C1  <- round(Mxshewhartrunsrules.arl(mus, c=c1,  type="1"), digits=2)
C12 <- round(Mxshewhartrunsrules.arl(mus, c=c12, type="12"), digits=2)
C13 <- round(Mxshewhartrunsrules.arl(mus, c=c13, type="13"), digits=2)
C14 <- round(Mxshewhartrunsrules.arl(mus, c=c14, type="14"), digits=2)
data.frame(mus, C1, C12, C13, C14)

## and the steady-state ARL
Mxshewhartrunsrules.ad <- Vectorize(xshewhartrunsrules.ad, "mu1")
C1  <- round(Mxshewhartrunsrules.ad(mus, c=c1,  type="1"), digits=2)
C12 <- round(Mxshewhartrunsrules.ad(mus, c=c12, type="12"), digits=2)
C13 <- round(Mxshewhartrunsrules.ad(mus, c=c13, type="13"), digits=2)
C14 <- round(Mxshewhartrunsrules.ad(mus, c=c14, type="14"), digits=2)
data.frame(mus, C1, C12, C13, C14)

Compute ARLs of CUSUM control charts

Description

Computation of the (zero-state) Average Run Length (ARL) for different types of CUSUM control charts monitoring normal mean.

Usage

xtcusum.arl(k, h, df, mu, hs = 0, sided="one", mode="tan", r=30)

Arguments

Details

xtcusum.arl determines the Average Run Length (ARL) by numerically solving the related ARL integral equation by means of the Nystroem method based on Gauss-Legendre quadrature.

Value

Returns a single value which resembles the ARL.

Author(s)

Sven Knoth

References

A. L. Goel, S. M. Wu (1971), Determination of A.R.L. and a contour nomogram for CUSUM charts to control normal mean, Technometrics 13, 221-230.

D. Brook, D. A. Evans (1972), An approach to the probability distribution of cusum run length, Biometrika 59, 539-548.

J. M. Lucas, R. B. Crosier (1982), Fast initial response for cusum quality-control schemes: Give your cusum a headstart, Technometrics 24, 199-205.

L. C. Vance (1986), Average run lengths of cumulative sum control charts for controlling normal means, Journal of Quality Technology 18, 189-193.

K.-H. Waldmann (1986), Bounds for the distribution of the run length of one-sided and two-sided CUSUM quality control schemes, Technometrics 28, 61-67.

R. B. Crosier (1986), A new two-sided cumulative quality control scheme, Technometrics 28, 187-194.

See Also

xtewma.arl for zero-state ARL computation of EWMA control charts and xtcusum.arl for the zero-state ARL of CUSUM for normal data.

Examples

## will follow

Compute steady-state ARLs of EWMA control charts, t distributed data

Description

Computation of the steady-state Average Run Length (ARL) for different types of EWMA control charts monitoring the mean of t distributed data.

Usage

xtewma.ad(l, c, df, mu1, mu0=0, zr=0, z0=0, sided="one", limits="fix",
steady.state.mode="conditional", mode="tan", r=40)

Arguments

Details

xtewma.ad determines the steady-state Average Run Length (ARL) by numerically solving the related ARL integral equation by means of the Nystroem method based on Gauss-Legendre quadrature and using the power method for deriving the largest in magnitude eigenvalue and the related left eigenfunction.

Value

Returns a single value which resembles the steady-state ARL.

Author(s)

Sven Knoth

References

R. B. Crosier (1986), A new two-sided cumulative quality control scheme, Technometrics 28, 187-194.

S. V. Crowder (1987), A simple method for studying run-length distributions of exponentially weighted moving average charts, Technometrics 29, 401-407.

J. M. Lucas and M. S. Saccucci (1990), Exponentially weighted moving average control schemes: Properties and enhancements, Technometrics 32, 1-12.

See Also

xtewma.arl for zero-state ARL computation and xewma.ad for the steady-state ARL for normal data.

Examples

## will follow

Compute ARLs of EWMA control charts, t distributed data

Description

Computation of the (zero-state) Average Run Length (ARL) for different types of EWMA control charts monitoring the mean of t distributed data.

Usage

xtewma.arl(l,c,df,mu,zr=0,hs=0,sided="two",limits="fix",mode="tan",q=1,r=40)

Arguments

Details

In case of the EWMA chart with fixed control limits, xtewma.arl determines the Average Run Length (ARL) by numerically solving the related ARL integral equation by means of the Nystroem method based on Gauss-Legendre quadrature. If limits is "vacl", then the method presented in Knoth (2003) is utilized. Other values (normal case) for limits are not yet supported.

Value

Except for the fixed limits EWMA charts it returns a single value which resembles the ARL. For fixed limits charts, it returns a vector of length q which resembles the ARL and the sequence of conditional expected delays for q=1 and q>1, respectively.

Author(s)

Sven Knoth

References

K.-H. Waldmann (1986), Bounds for the distribution of the run length of geometric moving average charts, Appl. Statist. 35, 151-158.

S. V. Crowder (1987), A simple method for studying run-length distributions of exponentially weighted moving average charts, Technometrics 29, 401-407.

J. M. Lucas and M. S. Saccucci (1990), Exponentially weighted moving average control schemes: Properties and enhancements, Technometrics 32, 1-12.

C. M. Borror, D. C. Montgomery, and G. C. Runger (1999), Robustness of the EWMA control chart to non-normality , Journal of Quality Technology 31, 309-316.

S. Knoth (2003), EWMA schemes with non-homogeneous transition kernels, Sequential Analysis 22, 241-255.

S. Knoth (2004), Fast initial response features for EWMA Control Charts, Statistical Papers 46, 47-64.

See Also

xewma.arl for zero-state ARL computation of EWMA control charts in the normal case.

Examples

##  Borror/Montgomery/Runger (1999), Table 3
lambda <- 0.1
cE <- 2.703
df <- c(4, 6, 8, 10, 15, 20, 30, 40, 50)
L0 <- rep(NA, length(df))
for ( i in 1:length(df) ) {
  L0[i] <- round(xtewma.arl(lambda, cE*sqrt(df[i]/(df[i]-2)), df[i], 0), digits=0)
}
data.frame(df, L0)

Compute RL quantiles of EWMA control charts

Description

Computation of quantiles of the Run Length (RL) for EWMA control charts monitoring normal mean.

Usage

xtewma.q(l, c, df, mu, alpha, zr=0, hs=0, sided="two", limits="fix", mode="tan",
q=1, r=40)

xtewma.q.crit(l, L0, df, mu, alpha, zr=0, hs=0, sided="two", limits="fix", mode="tan",
r=40, c.error=1e-12, a.error=1e-9, OUTPUT=FALSE)

Arguments

Details

Instead of the popular ARL (Average Run Length) quantiles of the EWMA stopping time (Run Length) are determined. The algorithm is based on Waldmann's survival function iteration procedure. If limits is "vacl", then the method presented in Knoth (2003) is utilized. For details see Knoth (2004).

Value

Returns a single value which resembles the RL quantile of order q.

Author(s)

Sven Knoth

References

F. F. Gan (1993), An optimal design of EWMA control charts based on the median run length, J. Stat. Comput. Simulation 45, 169-184.

S. Knoth (2003), EWMA schemes with non-homogeneous transition kernels, Sequential Analysis 22, 241-255.

S. Knoth (2004), Fast initial response features for EWMA Control Charts, Statistical Papers 46, 47-64.

K.-H. Waldmann (1986), Bounds for the distribution of the run length of geometric moving average charts, Appl. Statist. 35, 151-158.

See Also

xewma.q for RL quantile computation of EWMA control charts in the normal case.

Examples

## will follow

Compute the survival function of EWMA run length

Description

Computation of the survival function of the Run Length (RL) for EWMA control charts monitoring normal mean.

Usage

xtewma.sf(l, c, df, mu, n, zr=0, hs=0, sided="two", limits="fix", mode="tan", q=1, r=40)

Arguments

Details

The survival function P(L>n) and derived from it also the cdf P(L<=n) and the pmf P(L=n) illustrate the distribution of the EWMA run length. For large n the geometric tail could be exploited. That is, with reasonable large n the complete distribution is characterized. The algorithm is based on Waldmann's survival function iteration procedure. For varying limits and for change points after 1 the algorithm from Knoth (2004) is applied. For details see Knoth (2004).

Value

Returns a vector which resembles the survival function up to a certain point.

Author(s)

Sven Knoth

References

F. F. Gan (1993), An optimal design of EWMA control charts based on the median run length, J. Stat. Comput. Simulation 45, 169-184.

S. Knoth (2003), EWMA schemes with non-homogeneous transition kernels, Sequential Analysis 22, 241-255.

S. Knoth (2004), Fast initial response features for EWMA Control Charts, Statistical Papers 46, 47-64.

K.-H. Waldmann (1986), Bounds for the distribution of the run length of geometric moving average charts, Appl. Statist. 35, 151-158.

See Also

xewma.sf for survival function computation of EWMA control charts in the normal case.

Examples

## will follow

Compute ARLs of modified Shewhart control charts for AR(1) data with Student t residuals

Description

Computation of the (zero-state) Average Run Length (ARL) for modified Shewhart charts deployed to the original AR(1) data where the residuals follow a Student t distribution.

Usage

xtshewhart.ar1.arl(alpha, cS, df, delta=0, N1=50, N2=30, N3=2*N2, INFI=10, mode="tan")

Arguments

Details

Following the idea of Schmid (1995), 1-alpha times the data turns out to be an EWMA smoothing of the underlying AR(1) residuals. Hence, by combining the solution of the EWMA ARL integral equation and the stationary distribution of the AR(1) data (here Student t distribution is assumed) one gets easily the overall ARL.

Value

It returns a single value resembling the zero-state ARL of a modified Shewhart chart.

Author(s)

Sven Knoth

References

J. Andel, K. Hrach (2000). On calculation of stationary density of autoregressive processes. Kybernetika, Institute of Information Theory and Automation AS CR 36(3), 311-319.

H. Kramer, W. Schmid (2000). The influence of parameter estimation on the ARL of Shewhart type charts for time series. Statistical Papers 41(2), 173-196.

W. Schmid (1995). On the run length of a Shewhart chart for correlated data. Statistical Papers 36(1), 111-130.

See Also

xtewma.arl for zero-state ARL computation of EWMA control charts in case of Student t distributed data.

Examples

## will follow