## 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.c.of.quantile <- function(l, L0, mu, p, zr=0, hs=0, sided="one", mode="integer", r=40) {
c2 <- 0
a2 <- 0
while ( a2 < L0 ) {
c2 <- c2 + .5
a2 <- xewma.q(l, c2, mu, p, zr=zr, hs=hs, sided=sided, r=r)
if ( mode=="integer" ) a2 <- ceiling(a2)
}
c1 <- c2 - .5
a1 <- xewma.q(l, c1, mu, p, zr=zr, hs=hs, sided=sided, r=r)
a.error <- 1; c.error <- 1
while ( a.error>1e-6 && c.error>1e-12 ) {
c3 <- c1 + (L0 - a1)/(a2 - a1)*(c2 - c1)
a3 <- xewma.q(l, c3, mu, p, zr=zr, hs=hs, sided=sided, r=r)
if ( mode=="integer" ) a3 <- ceiling(a3)
c1 <- c2; c2 <- c3
a1 <- a2; a2 <- a3
a.error <- abs(a2 - L0); c.error <- abs(c2 - c1)
}
c3
}
XEWMA.c.of.quantile <- Vectorize("xewma.c.of.quantile", "l")
# (ii) re-calculate the thresholds and remove the standardization step
L0 <- 500
G.h.new <- XEWMA.c.of.quantile(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))Run the code above in your browser using DataLab