spc (version 0.7.1)

mewma.arl: 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)

Value

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

Arguments

l

smoothing parameter lambda of the MEWMA control chart.

cE

alarm threshold of the MEWMA control chart.

p

dimension of multivariate normal distribution.

delta

magnitude of the potential change, delta=0 refers to the in-control state.

hs

so-called headstart (enables fast initial response) -- must be non-negative.

r

number of quadrature nodes -- dimension of the resulting linear equation system for delta = 0. For non-zero delta this dimension is mostly r^2 (Markov chain approximation leads to some larger values). Caution: If ntype is set to "co" (collocation), then values of r larger than 20 lead to large computing times. For the other selections this would happen for values larger than 40.

ntype

choose the numerical algorithm to solve the ARL integral equation. For delta=0: Possible values are "gl", "gl2" (gauss-legendre, classic and with variables change: square), "co" (collocation, for delta > 0 with sin transformation), "ra" (radau), "cc" (clenshaw-curtis), "mc" (markov chain), and "sr" (simpson rule). For delta larger than 0, some more values besides the others are possible: "gl3", "gl4", "gl5" (gauss-legendre with a further change in variables: sin, tan, sinh), "co2", "co3" (collocation with some trimming and tan as quadrature stabilizing transformations, respectively). If it is set to NULL (the default), then for delta=0 then "gl2" is chosen. If delta larger than 0, then for p equal 2 or 4 "gl3" and for all other values "gl5" is taken. "ra" denotes the method used in Rigdon (1995a). "mc" denotes the Markov chain approximation.

type

switch between "cond" and "cycl" for differentiating between the conditional (no false alarm) and the cyclical (after false alarm re-start in hs), respectively.

n

number of quadrature nodes for Calculating the steady-state ARL integral(s).

qm0,qm1

number of collocation quadrature nodes for the out-of-control case (qm0 for the inner integral, qm1 for the outer one), that is, for positive delta, and for the in-control case (now only qm0 is deployed) if via ntype the collocation procedure is requested.

Author

Sven Knoth

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.

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
if (FALSE) {
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"))
}

# Runger/Prabhu (1996), p. 1704, Tab. 1
if (FALSE) {
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
}

# Knoth (2017), p. 85, Tab. 3, rows with p=3
if (FALSE) {
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) )
}

# Prabhu/Runger (1997), p. 13, Tab. 3
if (FALSE) {
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"))
}