# =============================================================
# Example 8.1  –  Upward NSR (Nonparametric Signed-Rank) chart
#                 Qiu (2014), Section 8.2.1, Figure 8.2
# =============================================================

# --- Parameters ---
set.seed(42)

n_batch <- 30         # total batches
m       <- 10         # batch size
mu0     <- 0          # known IC median
U_limit <- 49         # control limit (Table 8.1, ARL+_0 ~ 204)

# =============================================================
# 1. Generate data
#    IC:  batches 1-20 from t_3
#    OOC: batches 21-30 from t_3 + 1  (median shift of 1)
# =============================================================

X <- matrix(NA, nrow = n_batch, ncol = m)
for (i in 1:20)  X[i, ] <- rt(m, df = 3)
for (i in 21:30) X[i, ] <- rt(m, df = 3) + 1

# =============================================================
# 2. Compute the Wilcoxon signed-rank statistic  psi_i
#    psi_i = sum of sign(X_ij - mu0) * rank(|X_ij - mu0|)
# =============================================================

psi <- numeric(n_batch)

for (i in 1:n_batch) {
  d      <- X[i, ] - mu0
  psi[i] <- sum(sign(d) * rank(abs(d)))
}

# --- OOC flags ---
ooc <- psi >= U_limit

cat("First signal at i =",
    ifelse(any(ooc), min(which(ooc)), NA), "\n")

# --- Colours ---
ic_col  <- "dodgerblue2"
ooc_col <- "firebrick2"

# =============================================================
# 3. Plot: 2 rows
#    Row 1: raw data (all observations by batch)
#    Row 2: upward NSR chart (psi_i with control limit)
# =============================================================

pdf("fig_ex8_1_nsr.pdf", width = 10.5, height = 6.5)
par(mfrow = c(2, 1), mar = c(4, 4, 2.5, 1))

# --- Panel (a): raw batch data ---
plot(rep(1:n_batch, each = m), as.vector(t(X)),
     pch = 16, col = adjustcolor("grey50", 0.5), cex = 0.7,
     xlab = "i", ylab = "", bty = "l",
     main = expression("(a) Observed data from " * t[3] *
                         " with median shift at " * i == 21))

# --- Panel (b): upward NSR chart ---
plot(1:n_batch, psi, type = "b", pch = 16, col = ic_col,
     xlab = "i", ylab = expression(psi[i]), bty = "l",
     main = expression("(b) Upward NSR chart, " * U == 49))
abline(h = U_limit, lty = 2)
points(which(ooc), psi[ooc], pch = 16, col = ooc_col)

dev.off()
cat("Saved figure: fig_ex8_1_nsr.pdf\n")



# =============================================================
# Example 8.3  –  Upward NSR CUSUM chart
#                 Qiu (2014), Section 8.2.2, Figure 8.4
# =============================================================

# --- Parameters ---
set.seed(42)

n_batch <- 30         # total batches
m       <- 6          # batch size
mu0     <- 0          # known IC median
k_cusum <- 8          # allowance constant (integer, Table 8.4)
h_cusum <- 10         # control limit (integer)
delta   <- 1          # location shift size

# =============================================================
# 1. Generate data
#    IC:  batches 1-20 from t_3  (symmetric, heavy-tailed)
#    OOC: batches 21-30 from t_3 + 1
# =============================================================

X <- matrix(NA, nrow = n_batch, ncol = m)
for (i in 1:20)  X[i, ] <- rt(m, df = 3)
for (i in 21:30) X[i, ] <- rt(m, df = 3) + delta

# =============================================================
# 2. Compute Wilcoxon signed-rank statistic psi_n
#    and the upward CUSUM recursion C_n^+
# =============================================================

psi <- numeric(n_batch)
Cp  <- numeric(n_batch)

for (n in 1:n_batch) {
  d      <- X[n, ] - mu0
  psi[n] <- sum(sign(d) * rank(abs(d)))
  
  prev   <- ifelse(n == 1, 0, Cp[n - 1])
  Cp[n]  <- max(0, prev + psi[n] - k_cusum)
}

# --- OOC flags ---
ooc <- Cp > h_cusum

cat("First signal at n =",
    ifelse(any(ooc), min(which(ooc)), NA), "\n")

# --- Colours ---
ic_col  <- "dodgerblue2"
ooc_col <- "firebrick2"

# =============================================================
# 3. Plot: 2 rows
#    Row 1: raw data (all observations by batch)
#    Row 2: upward NSR CUSUM (C_n^+ with control limit)
# =============================================================

pdf("fig_ex8_3_nsr_cusum.pdf", width = 10.5, height = 6.5)
par(mfrow = c(2, 1), mar = c(4, 4, 2.5, 1))

# --- Panel (a): raw batch data ---
plot(rep(1:n_batch, each = m), as.vector(t(X)),
     pch = 16, col = adjustcolor("grey50", 0.5), cex = 0.7,
     xlab = "n", ylab = "", bty = "l",
     main = expression("(a) Phase II batch data, " * m == 6))

# --- Panel (b): upward NSR CUSUM ---
plot(1:n_batch, Cp, type = "b", pch = 16, col = ic_col,
     xlab = "n", ylab = expression(C[n]^"+"), bty = "l",
     main = expression("(b) Upward NSR CUSUM, " *
                         italic(k) * " = 8, " * italic(h) * " = 10"))
abline(h = h_cusum, lty = 2)
points(which(ooc), Cp[ooc], pch = 16, col = ooc_col)

dev.off()
cat("Saved figure: fig_ex8_3_nsr_cusum.pdf\n")


# =============================================================
# EWMA-Lepage chart  –  joint location and scale monitoring
# Mukherjee and Chakraborti (2012), QREI, 28, 335-352
# =============================================================

# --- Parameters ---
set.seed(42)

M          <- 100      # reference sample size
n_batch    <- 50       # Phase II batches
m          <- 10       # batch size
shift_time <- 31       # shift starts here
mu_oc      <- 1        # OC mean (IC = 0)
sigma_oc   <- 2        # OC standard deviation (IC = 1)
lambda_ew  <- 0.05     # EWMA smoothing constant

# =============================================================
# 1. Generate reference sample (IC)
# =============================================================

Y_ref <- rnorm(M, mean = 0, sd = 1)

# =============================================================
# 2. Generate Phase II batch data
#    IC:  batches 1-30 from N(0,1)
#    OOC: batches 31-50 from N(1,4)  (joint shift)
# =============================================================

X <- matrix(NA, nrow = n_batch, ncol = m)
for (i in 1:(shift_time - 1)) X[i, ] <- rnorm(m, 0, 1)
for (i in shift_time:n_batch)  X[i, ] <- rnorm(m, mu_oc, sigma_oc)

# =============================================================
# 3. Ansari-Bradley scores and IC moments
#    Score for rank r in combined sample of size N: min(r, N+1-r)
# =============================================================

N_pool <- M + m

# Wilcoxon rank-sum moments
mu_W    <- m * (N_pool + 1) / 2
sigma_W <- sqrt(m * M * (N_pool + 1) / 12)

# Ansari-Bradley moments
ab_scores_all <- pmin(1:N_pool, (N_pool + 1) - 1:N_pool)
mu_A    <- m * mean(ab_scores_all)
sigma_A <- sqrt(m * M / (N_pool * (N_pool - 1)) *
                  sum((ab_scores_all - mean(ab_scores_all))^2))

# =============================================================
# 4. Compute Lepage statistic and EWMA recursion
# =============================================================

Ln <- numeric(n_batch)
En <- numeric(n_batch)

for (n in 1:n_batch) {
  # Pool reference + batch, rank everything
  r_all   <- rank(c(Y_ref, X[n, ]))
  r_batch <- r_all[(M + 1):N_pool]
  
  # Wilcoxon rank-sum (location)
  W_n <- sum(r_batch)
  
  # Ansari-Bradley statistic (scale)
  A_n <- sum(pmin(r_batch, (N_pool + 1) - r_batch))
  
  # Lepage = Z_W^2 + Z_A^2
  Ln[n] <- ((W_n - mu_W) / sigma_W)^2 + ((A_n - mu_A) / sigma_A)^2
  
  # EWMA update
  prev  <- ifelse(n == 1, 0, En[n - 1])
  En[n] <- lambda_ew * Ln[n] + (1 - lambda_ew) * prev
}

# =============================================================
# 5. Control limit h_L  (bisection + Monte Carlo, ARL_0 ~ 200)
# =============================================================

target_arl <- 200
n_sim      <- 5000
max_len    <- 2000
Y_ref_cal  <- rnorm(M, 0, 1)

ab_sc <- pmin(1:N_pool, (N_pool + 1) - 1:N_pool)

sim_rl <- function(h) {
  E_prev <- 0
  for (t in 1:max_len) {
    r_all <- rank(c(Y_ref_cal, rnorm(m, 0, 1)))
    r_b   <- r_all[(M + 1):N_pool]
    L_t   <- ((sum(r_b) - mu_W) / sigma_W)^2 +
      ((sum(pmin(r_b, (N_pool + 1) - r_b)) - mu_A) / sigma_A)^2
    E_cur <- lambda_ew * L_t + (1 - lambda_ew) * E_prev
    if (E_cur > h) return(t)
    E_prev <- E_cur
  }
  return(max_len)
}

h_lo <- 1.0;  h_hi <- 8.0

cat("Calibrating control limit (this may take a few minutes)...\n")
for (iter in 1:20) {
  h_mid   <- (h_lo + h_hi) / 2
  arl_est <- mean(replicate(n_sim, sim_rl(h_mid)))
  cat(sprintf("  iter %2d: h = %.4f, ARL = %.1f\n", iter, h_mid, arl_est))
  if (arl_est < target_arl) h_lo <- h_mid else h_hi <- h_mid
  if (abs(h_hi - h_lo) < 0.01) break
}
h_L <- (h_lo + h_hi) / 2
cat(sprintf("Calibrated h_L = %.4f\n", h_L))

# --- OOC flags ---
ooc <- En > h_L

cat("First signal at n =",
    ifelse(any(ooc), min(which(ooc)), NA), "\n")

# --- Colours ---
ic_col  <- "dodgerblue2"
ooc_col <- "firebrick2"

# =============================================================
# 6. Plot: 2 rows
#    Row 1: raw data (all observations by batch)
#    Row 2: EWMA-Lepage chart (E_n with control limit)
# =============================================================

pdf("fig_ewma_lepage.pdf", width = 10.5, height = 6.5)
par(mfrow = c(2, 1), mar = c(4, 4, 2.5, 1))

# --- Panel (a): raw batch data ---
plot(rep(1:n_batch, each = m), as.vector(t(X)),
     pch = 16, col = adjustcolor("grey50", 0.5), cex = 0.7,
     xlab = "n", ylab = "", bty = "l",
     main = expression("(a) Phase II data: " * mu *
                         " shifts to 1, " * sigma *
                         " shifts to 2 at " * n == 31))
abline(v = shift_time - 0.5, lty = 3, col = "grey40")

# --- Panel (b): EWMA-Lepage chart ---
plot(1:n_batch, En, type = "b", pch = 16, col = ic_col,
     xlab = "n", ylab = expression(E[n]), bty = "l",
     main = expression("(b) EWMA-Lepage chart, " *
                         lambda * " = 0.05, ARL"[0] * " = 200"))
abline(h = h_L, lty = 2)
points(which(ooc), En[ooc], pch = 16, col = ooc_col)

dev.off()
cat("Saved figure: fig_ewma_lepage.pdf\n")