## 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
if (FALSE) 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 )
if (FALSE) 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"))
Run the code above in your browser using DataLab