# ============================================================= # Example: Q-chart for individual Normal observations # Quesenberry (1991), J. Qual. Technol., 23, 213-224 # ============================================================= # --- Parameters --- set.seed(42) n_obs <- 20 # total observations outlier_pos <- 11 # position of the planted outlier outlier_val <- 4.5 # outlier magnitude (in sigma units) L <- 3.023 # control limit (ARL_0 = 370) # ============================================================= # 1. Generate data # All observations from N(0,1); one outlier at n = 11 # ============================================================= X <- rnorm(n_obs, mean = 0, sd = 1) X[outlier_pos] <- outlier_val # ============================================================= # 2. Compute Q-statistics (defined for n >= 3) # T_n = (X_n - Xbar_{n-1}) / s_{n-1} # Q_n = Phi^{-1}[ F_{t_{n-2}}( sqrt((n-1)/n) * T_n ) ] # ============================================================= Q <- rep(NA_real_, n_obs) for (n in 3:n_obs) { xbar_prev <- mean(X[1:(n - 1)]) s_prev <- sd(X[1:(n - 1)]) T_n <- (X[n] - xbar_prev) / s_prev t_scaled <- sqrt((n - 1) / n) * T_n p_n <- pt(t_scaled, df = n - 2) p_n <- pmin(pmax(p_n, .Machine$double.eps), 1 - .Machine$double.eps) Q[n] <- qnorm(p_n) } # --- Signal detection --- signal <- !is.na(Q) & abs(Q) > L cat("Signal(s) at observation(s):", ifelse(any(signal), paste(which(signal), collapse = ", "), "none"), "\n") # --- Colours --- ic_col <- "dodgerblue2" ooc_col <- "firebrick2" # ============================================================= # 3. Plot: 2 rows # Row 1: raw observations (outlier in red) # Row 2: Q-chart with control limits # ============================================================= pdf("fig_qchart_normal.pdf", width = 10.5, height = 6.5) par(mfrow = c(2, 1), mar = c(4, 4, 2.5, 1)) # --- Panel (a): raw observations --- pt_col <- ifelse(1:n_obs == outlier_pos, ooc_col, ic_col) plot(1:n_obs, X, type = "b", pch = 16, col = pt_col, cex = 0.9, xlab = "n", ylab = "X", bty = "l", main = expression("(a) Individual observations: " * italic(N)(0*","*1) * " with a 4.5" * sigma * " outlier at " * italic(n) == 11)) abline(h = 0, lty = 1, col = "grey70") # --- Panel (b): Q-chart --- pt_col_q <- ifelse(signal & !is.na(Q), ooc_col, ic_col) plot(1:n_obs, Q, type = "b", pch = 16, col = pt_col_q, cex = 0.9, xlab = "n", ylab = "Q", bty = "l", main = bquote("(b) Q-chart, " * italic(L) == .(L) * " (ARL"[0] * " = 370)")) abline(h = L, lty = 2) abline(h = -L, lty = 2) abline(h = 0, lty = 1, col = "grey70") dev.off() cat("Saved: fig_qchart_normal.pdf\n") # ============================================================= # Example: Self-starting CUSUM for individual Normal observations # Both mu_0 and sigma unknown, estimated via Q-statistics # Hawkins and Olwell (1998) # ============================================================= # --- Parameters --- n_obs <- 30 # total observations tau <- 11 # change-point (shift starts here) delta <- 1 # mean shift size (in sigma units) k <- 0.25 # CUSUM reference value h <- 4.5 # control limit # ============================================================= # 1. Helper functions # ============================================================= # --- Q-statistics: Q_n ~ N(0,1) under IC --- compute_Q <- function(X) { n <- length(X) Q <- rep(NA_real_, n) for (i in 3:n) { xbar <- mean(X[1:(i - 1)]) s <- sd(X[1:(i - 1)]) T_i <- (X[i] - xbar) / s p <- pt(sqrt((i - 1) / i) * T_i, df = i - 2) p <- pmin(pmax(p, .Machine$double.eps), 1 - .Machine$double.eps) Q[i] <- qnorm(p) } return(Q) } # --- One-sided CUSUM on Q-statistics --- compute_cusum <- function(Q, k) { n <- length(Q) C <- rep(NA_real_, n) C_prev <- 0 for (i in seq_len(n)) { if (!is.na(Q[i])) { C[i] <- max(0, C_prev + Q[i] - k) C_prev <- C[i] } } return(C) } # ============================================================= # 2. Find a seed with no false alarm before tau, # and at least one signal at or after tau # ============================================================= chosen_seed <- NA_integer_ for (s in 1:50000) { set.seed(s) X_try <- c(rnorm(tau - 1, 0, 1), rnorm(n_obs - tau + 1, delta, 1)) C_try <- compute_cusum(compute_Q(X_try), k) sig <- !is.na(C_try) & C_try > h if (!any(sig[1:(tau - 1)], na.rm = TRUE) && any(sig[tau:n_obs], na.rm = TRUE)) { chosen_seed <- s break } } if (is.na(chosen_seed)) stop("No suitable seed found. Try lowering h.") cat("Chosen seed:", chosen_seed, "\n") # ============================================================= # 3. Generate final dataset # ============================================================= set.seed(chosen_seed) X <- c(rnorm(tau - 1, 0, 1), rnorm(n_obs - tau + 1, delta, 1)) # ============================================================= # 4. Compute Q-statistics and self-starting CUSUM # ============================================================= Q <- compute_Q(X) C <- compute_cusum(Q, k) # --- Signal detection --- signal <- !is.na(C) & C > h cat("Signal(s) at observation(s):", ifelse(any(signal), paste(which(signal), collapse = ", "), "none"), "\n") # --- Colours --- ic_col <- "dodgerblue2" ooc_col <- "firebrick2" # ============================================================= # 5. Plot: 2 rows # Row 1: raw observations (IC blue, OOC red) # Row 2: self-starting CUSUM with control limit # ============================================================= pdf("fig_ss_cusum.pdf", width = 10.5, height = 6.5) par(mfrow = c(2, 1), mar = c(4, 4, 2.5, 1)) # --- Panel (a): raw observations --- pt_col <- ifelse(1:n_obs < tau, ic_col, ooc_col) plot(1:n_obs, X, type = "b", pch = 16, col = pt_col, cex = 0.9, xlab = "n", ylab = "X", bty = "l", main = paste0("(a) Individual observations: N(0,1) up to n = ", tau - 1, ", then N(", delta, ",1)")) abline(h = 0, lty = 1, col = "grey70") abline(v = tau - 0.5, lty = 3, col = "grey40") # --- Panel (b): self-starting CUSUM --- pt_col_c <- ifelse(signal & !is.na(C), ooc_col, ic_col) plot(1:n_obs, C, type = "b", pch = 16, col = pt_col_c, cex = 0.9, xlab = "n", ylab = "C", bty = "l", main = paste0("(b) Self-starting CUSUM, k = ", k, ", h = ", h)) abline(h = h, lty = 2) abline(h = 0, lty = 1, col = "grey70") abline(v = tau - 0.5, lty = 3, col = "grey40") dev.off() cat("Saved: fig_ss_cusum.pdf\n")