Controle Estatístico de Processos: Gráficos CUSUM e MMEP (EWMA)

Mapa do Capítulo 9

Agenda

• Limitação estrutural do gráfico de Shewhart: memória curta
• Gráfico CUSUM: motivação, forma tabular, parâmetros \(K\) e \(H\)
• Planejamento do CUSUM via ARL (Hawkins, Siegmund), CUSUM padronizado, RIR/headstart
• Gráfico MMEP (EWMA): definição recursiva, limites exatos e assintóticos
• Planejamento do MMEP (Lucas & Saccucci) e robustez à não normalidade
• Gráfico da Média Móvel (não ponderada) — Seção 9.3
• Exercícios aplicados no R com o pacote qcc

Configuração do Ambiente em R

library(ggplot2)
library(qcc)
library(patchwork)

pal <- c(azul   = "#7ecff7",
         verde  = "#4f9d6e",
         ambar  = "#f7c87e",
         roxo   = "#c87ef7",
         coral  = "#f07050",
         neutro = "#aab4c8")

tema_cep <- function(base_size = 13) {
  theme_minimal(base_size = base_size) +
    theme(plot.title = element_text(face = "bold"),
          panel.grid.minor = element_blank())
}

plot_carta <- function(stat, center, lic, lsc, titulo, ylab = "Estatística") {
  df <- data.frame(t = seq_along(stat), y = stat)
  df$fora <- df$y > lsc | df$y < lic
  ggplot(df, aes(t, y)) +
    geom_line(color = pal["neutro"], linewidth = .6) +
    geom_point(aes(color = fora), size = 2.6) +
    scale_color_manual(values = c(`FALSE` = unname(pal["azul"]),
                                   `TRUE`  = unname(pal["coral"])), guide = "none") +
    geom_hline(yintercept = c(lsc, center, lic),
               color = c(pal["coral"], pal["verde"], pal["coral"]),
               linetype = c("dashed", "solid", "dashed"), linewidth = .8) +
    labs(title = titulo, x = "Amostra", y = ylab) +
    tema_cep()
}

Funções Auxiliares: CUSUM e MMEP

# Gráfico do CUSUM tabular (C+ e C- em torno de zero)
plot_cusum <- function(Cmais, Cmenos, H, titulo = "Gráfico CUSUM Tabular") {
  df <- data.frame(
    t     = rep(seq_along(Cmais), 2),
    valor = c(Cmais, -Cmenos),
    lado  = rep(c("C+", "C-"), each = length(Cmais))
  )
  df$fora <- (df$lado == "C+" & df$valor >  H) |
             (df$lado == "C-" & df$valor < -H)
  ggplot(df, aes(t, valor, color = lado)) +
    geom_line(linewidth = .6) +
    geom_point(aes(shape = fora), size = 2.4) +
    geom_hline(yintercept = c(H, -H, 0),
               color = c(pal["coral"], pal["coral"], pal["neutro"]),
               linetype = c("dashed", "dashed", "solid")) +
    scale_color_manual(values = c("C+" = unname(pal["azul"]), "C-" = unname(pal["roxo"]))) +
    scale_shape_manual(values = c(`FALSE` = 16, `TRUE` = 17), guide = "none") +
    labs(title = titulo, x = "Amostra", y = expression(C[i]), color = NULL) +
    tema_cep()
}

# Gráfico MMEP com limites de controle que variam com i (limites exatos)
plot_mmep <- function(z, centro, lic, lsc, titulo = "Gráfico MMEP") {
  df <- data.frame(t = seq_along(z), z = z, lic = lic, lsc = lsc)
  df$fora <- df$z > df$lsc | df$z < df$lic
  ggplot(df, aes(t, z)) +
    geom_ribbon(aes(ymin = lic, ymax = lsc), fill = pal["neutro"], alpha = .15) +
    geom_line(aes(y = lic), color = pal["coral"], linetype = "dashed", linewidth = .6) +
    geom_line(aes(y = lsc), color = pal["coral"], linetype = "dashed", linewidth = .6) +
    geom_hline(yintercept = centro, color = pal["verde"], linewidth = .8) +
    geom_line(color = pal["neutro"], linewidth = .6) +
    geom_point(aes(color = fora), size = 2.6) +
    scale_color_manual(values = c(`FALSE` = unname(pal["azul"]), `TRUE` = unname(pal["coral"])),
                        guide = "none") +
    labs(title = titulo, x = "Amostra", y = expression(z[i])) +
    tema_cep()
}

A Limitação Estrutural do Gráfico de Shewhart

• Usa apenas a informação do último ponto plotado
• Ignora toda a informação contida na sequência histórica de pontos
• Consequência: insensível a mudanças pequenas no processo (tipicamente \(< 1{,}5\sigma\))
• Testes de sequência e limites de alerta são paliativos — reduzem a simplicidade de interpretação e o CMS sob controle
• Duas alternativas eficazes para a Fase II de monitoramento: CUSUM e MMEP

Dados Motivadores (Tabela 9.1 do Livro)

• 20 primeiras observações: \(X \sim N(\mu=10,\ \sigma=1)\) — processo sob controle
• 10 últimas observações: \(X \sim N(\mu=11,\ \sigma=1)\) — mudança de \(1\sigma\) na média

x_tab91 <- c(9.45, 7.99, 9.29, 11.66, 12.16, 10.18, 8.04, 11.46, 9.20, 10.34,
             9.03, 11.47, 10.51, 9.40, 10.08, 9.37, 10.62, 10.31, 8.52, 10.84,
             10.90, 9.33, 12.29, 11.50, 10.60, 11.08, 10.38, 11.62, 11.31, 10.52)
mu0 <- 10; sigma <- 1
length(x_tab91)

[1] 30

O Gráfico de Shewhart Falha em Detectar a Mudança

plot_carta(x_tab91, center = mu0, lic = mu0 - 3*sigma, lsc = mu0 + 3*sigma,
           titulo = "Carta de Shewhart (Individual) — Tabela 9.1", ylab = "x")

• Nenhum ponto ultrapassa os limites \(3\sigma\) — mesmo havendo, de fato, uma mudança real na média
• Motivo: o gráfico de Shewhart para médias só é eficaz para mudanças de \(1{,}5\sigma\) a \(2\sigma\) ou mais

Duas Alternativas ao Gráfico de Shewhart

• CUSUM (Cumulative Sum) — soma cumulativa dos desvios em relação ao alvo
• MMEP / EWMA — média móvel exponencialmente ponderada de todas as observações passadas
• Ambos combinam informação de várias amostras, e não apenas da última
• São chamados, em conjunto, de gráficos de controle ponderados pelo tempo
• Particularmente eficazes com amostras de tamanho \(n=1\) — ideais para indústrias químicas, de processo contínuo e manufatura com medição automática

Gráfico de Controle da Soma Cumulativa (CUSUM)

\[C_i = \sum_{j=1}^{i}(\bar x_j - \mu_0) = (\bar x_i - \mu_0) + C_{i-1}\]

• Proposto por Page (1954); estudado por Ewan (1963), Lucas (1976, 1982), Hawkins (1981, 1993)
• Se o processo permanece sob controle em \(\mu_0\): \(C_i\) é um passeio aleatório com média zero
• Se a média se desloca para \(\mu_1 > \mu_0\): tendência positiva se desenvolve em \(C_i\)
• Se a média se desloca para \(\mu_1 < \mu_0\): tendência negativa se desenvolve em \(C_i\)

CUSUM Simples Aplicado à Tabela 9.1

Ci <- cumsum(x_tab91 - mu0)
df_c <- data.frame(t = seq_along(Ci), Ci = Ci)
ggplot(df_c, aes(t, Ci)) +
  geom_line(color = pal["azul"], linewidth = .8) +
  geom_point(color = pal["azul"], size = 2) +
  geom_hline(yintercept = 0, color = pal["neutro"], linetype = "dashed") +
  labs(title = "Soma Cumulativa dos Desvios em Relação ao Alvo",
       x = "Amostra", y = expression(C[i])) +
  tema_cep()

• Nas 20 primeiras observações, \(C_i\) flutua perto de zero; nas 10 últimas, desenvolve-se uma tendência de subida clara

Duas Representações do CUSUM

• A partir daqui, o foco é inteiramente na forma tabular

O CUSUM Tabular: Ideia

• Acumula os desvios de \(\mu_0\) acima do alvo em \(C^+\) (cusum unilateral superior)
• Acumula os desvios de \(\mu_0\) abaixo do alvo em \(C^-\) (cusum unilateral inferior)
• Ambas as estatísticas são “resetadas” em zero quando ficam negativas — só acumulam evidência de mudança na direção de interesse

\[C_i^+ = \max\bigl[0,\ x_i - (\mu_0+K) + C_{i-1}^+\bigr] \qquad C_i^- = \max\bigl[0,\ (\mu_0-K) - x_i + C_{i-1}^-\bigr]\]

com \(C_0^+ = C_0^- = 0\).

O Parâmetro \(K\): Valor de Referência

\[\mu_1 = \mu_0 + \delta\sigma \qquad \delta = \frac{|\mu_1-\mu_0|}{\sigma} \qquad K = \frac{\delta}{2}\sigma = k\sigma\]

• \(K\) representa o ponto médio entre o valor-alvo \(\mu_0\) e o valor fora de controle \(\mu_1\) que se deseja detectar rapidamente
• Escolha típica na prática: \(k = 1/2\), ou seja, projetado para detectar uma mudança de \(1\sigma\)

O Parâmetro \(H\): Intervalo de Decisão

\[H = h\sigma\]

• Se \(C_i^+\) ou \(C_i^-\) excederem \(H\), o processo é declarado fora de controle
• Valor comumente usado na prática: \(h = 4\) ou \(h = 5\)
• \(N^+\) e \(N^-\): número de períodos consecutivos em que \(C^+\) e \(C^-\) permaneceram não nulos — úteis para estimar em qual amostra a mudança ocorreu

Exemplo: CUSUM Tabular (Montgomery)

K <- 0.5; H <- 5

• Os cálculos manuais das duas primeiras amostras (Montgomery, p. 397):

\[C_1^+ = \max[0;\ 9{,}45-(10{,}5)+0] = 0 \qquad C_1^- = \max[0;\ (9{,}5)-9{,}45+0] = 0{,}05\] \[C_2^+ = \max[0;\ 7{,}99-(10{,}5)+0] = 0 \qquad C_2^- = \max[0;\ (9{,}5)-7{,}99+0{,}05] = 1{,}56\]

Implementando o CUSUM Tabular no R

Cmais <- Cmenos <- numeric(length(x_tab91))
Cmais[1]  <- max(0,  x_tab91[1] - (mu0 + K))
Cmenos[1] <- max(0, (mu0 - K)  - x_tab91[1])
for (i in 2:length(x_tab91)) {
  Cmais[i]  <- max(0,  x_tab91[i] - (mu0 + K) + Cmais[i-1])
  Cmenos[i] <- max(0, (mu0 - K)  - x_tab91[i] + Cmenos[i-1])
}
round(cbind(amostra = 1:30, x = x_tab91, `C+` = Cmais, `C-` = Cmenos), 3)[1:8, ]

     amostra     x   C+   C-
[1,]       1  9.45 0.00 0.05
[2,]       2  7.99 0.00 1.56
[3,]       3  9.29 0.00 1.77
[4,]       4 11.66 1.16 0.00
[5,]       5 12.16 2.82 0.00
[6,]       6 10.18 2.50 0.00
[7,]       7  8.04 0.04 1.46
[8,]       8 11.46 1.00 0.00

Gráfico CUSUM Tabular Resultante

plot_cusum(Cmais, Cmenos, H, titulo = "CUSUM Tabular — Dados da Tabela 9.1 (K=0,5 ; H=5)")

• \(C^+\) ultrapassa \(H=5\) por volta da amostra 26–27 — sinal de processo fora de controle, coerente com a mudança real introduzida na amostra 21

Confirmando com o Pacote qcc

q_cusum <- cusum(x_tab91, center = mu0, std.dev = sigma,
                  se.shift = 2*K, decision.interval = H, plot = FALSE)
summary(q_cusum)

Call:
cusum(data = x_tab91, center = mu0, std.dev = sigma, decision.interval = H,     se.shift = 2 * K, plot = FALSE)

cusum chart for x_tab91 

Summary of group statistics:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  7.990   9.377  10.445  10.315  11.252  12.290 

Group sample size:  1
Number of groups:  30
Center of group statistics:  10
Standard deviation:  1 

Decision interval (std.err.): 5 
Shift detection  (std. err.): 1 

• se.shift no qcc corresponde a \(\delta\) (mudança em unidades de \(\sigma\), e não a \(K=\delta/2\))
• decision.interval corresponde a \(h\)

Recomendações para o Planejamento do CUSUM

• Parâmetros \(K=k\sigma\) e \(H=h\sigma\) devem ser escolhidos para fornecer um bom CMS (ARL)
• Meta usual: \(ARL_0 \approx 370\) (mesmo desempenho sob controle do gráfico \(3\sigma\) de Shewhart)
• Na prática, bons resultados são obtidos com \(h=4\) ou \(h=5\) e \(k=1/2\)

Tabela de Hawkins (1993): \((k,h)\) para \(ARL_0 = 370\)

tab_hawkins <- data.frame(
  k = c(0.25, 0.5, 0.75, 1.0, 1.25, 1.5),
  h = c(8.01, 4.77, 3.34, 2.52, 1.99, 1.61)
)
tab_hawkins

     k    h
1 0.25 8.01
2 0.50 4.77
3 0.75 3.34
4 1.00 2.52
5 1.25 1.99
6 1.50 1.61

• Quanto maior \(k\) (mudança maior a detectar), menor o \(h\) necessário para manter o mesmo \(ARL_0\)

Aproximação de Siegmund (1985) para o ARL Unilateral

\[ARL \approx \frac{\exp(-2\Delta b) - 1 + 2\Delta b}{2\Delta^2}, \qquad b = h + 1{,}166\]

• Para \(C^+\): \(\Delta = \delta^\ast - k\); para \(C^-\): \(\Delta = -\delta^\ast - k\)
• Se \(\Delta = 0\), usa-se \(ARL = b^2\)
• \(\delta^\ast\) é a mudança na média (em unidades de \(\sigma\)) para a qual se deseja calcular o ARL

arl_siegmund <- function(delta_star, k, h, lado = c("+", "-")) {
  lado <- match.arg(lado)
  Delta <- if (lado == "+") delta_star - k else -delta_star - k
  b <- h + 1.166
  if (Delta == 0) return(b^2)
  (exp(-2 * Delta * b) - 1 + 2 * Delta * b) / (2 * Delta^2)
}

Exemplo: ARL Bilateral sob Controle (\(k=1/2,\ h=5,\ \delta^\ast=0\))

arl_mais <- arl_siegmund(delta_star = 0, k = 0.5, h = 5, lado = "+")
arl_menos <- arl_siegmund(delta_star = 0, k = 0.5, h = 5, lado = "-")
arl0 <- 1 / (1/arl_mais + 1/arl_menos)
cat(sprintf("ARL+ = ARL- = %.1f\nARL0 (bilateral) = %.1f\n", arl_mais, arl0))

ARL+ = ARL- = 938.2
ARL0 (bilateral) = 469.1

• Resultado bate com a Tabela 9.3 do livro: \(ARL_0 \approx 469\) (a aproximação de Siegmund fica próxima do valor exato, \(370{,}4\), mas superestima um pouco)

ARL: CUSUM vs. Gráfico de Shewhart

deltas <- seq(0, 4, by = 0.25)
arl_cusum_h5 <- sapply(deltas, function(d)
  1 / (1/arl_siegmund(d, 0.5, 5, "+") + 1/arl_siegmund(d, 0.5, 5, "-")))
arl_shewhart <- sapply(deltas, function(d)
  1 / (1 - (pnorm(3 - d) - pnorm(-3 - d))))

df_arl <- data.frame(
  delta = rep(deltas, 2),
  arl   = c(arl_cusum_h5, arl_shewhart),
  metodo = rep(c("CUSUM (h=5, k=0,5)", "Shewhart (3-sigma)"), each = length(deltas))
)
ggplot(df_arl, aes(delta, arl, color = metodo)) +
  geom_line(linewidth = 1) + geom_point(size = 1.6) +
  scale_y_log10() +
  scale_color_manual(values = c(unname(pal["azul"]), unname(pal["coral"]))) +
  labs(title = "Comprimento Medio da Sequencia: CUSUM vs. Shewhart",
       x = expression(delta~"(mudanca em unidades de sigma)"),
       y = "ARL (escala log)", color = NULL) +
  tema_cep()

Leitura do Gráfico de ARL

• Para \(\sigma=1\) (mudanças pequenas), o CUSUM detecta muito mais rápido que o Shewhart
• À medida que \(\delta\) (tamanho da mudança) aumenta, os dois métodos convergem
• Conclusão prática: CUSUM vale a pena justamente quando o interesse está em mudanças pequenas a moderadas

O CUSUM Padronizado

\[y_i = \frac{x_i - \mu_0}{\sigma}\]

\[C_i^+ = \max\bigl[0,\ y_i - k + C_{i-1}^+\bigr] \qquad C_i^- = \max\bigl[0,\ -k - y_i + C_{i-1}^-\bigr]\]

• Vantagem principal: os mesmos valores de \(k\) e \(h\) podem ser usados em diferentes gráficos CUSUM, pois a escolha desses parâmetros deixa de depender da escala da variável

Subgrupos Racionais no CUSUM

• A extensão para subgrupos de tamanho \(n>1\) é imediata: basta trocar \(x_i\) por \(\bar x_i\), e \(\sigma\) por \(\sigma_{\bar x}=\sigma/\sqrt{n}\)
• Importante (Montgomery): ao contrário do gráfico de Shewhart, usar \(n>1\) não melhora o desempenho do CUSUM
• Exemplo: entre amostrar \(n=1\) a cada 30 min ou \(n=5\) a cada 2,5h, o CUSUM funciona melhor, em geral, com \(n=1\)
• Subgrupos maiores só se justificam por economia de escala ou outra razão prática relevante

Melhorando o CUSUM para Grandes Mudanças

• Procedimento combinado CUSUM–Shewhart (Lucas, 1982)
• Construir um gráfico de Shewhart sobre \(C^+\) e \(C^-\), com limites de \(3{,}5\sigma\) (o acréscimo de \(0{,}5\sigma\) compensa o interesse específico em mudanças grandes)
• O CUSUM fica “responsável” pelas mudanças pequenas; o Shewhart, pelas grandes
• Um sinal fora de controle em qualquer um dos dois gráficos aciona a investigação

Resposta Inicial Rápida (RIR) — Headstart

• Procedimento de Lucas & Crosier (1982) para melhorar a sensibilidade do CUSUM logo no início do monitoramento
• Ideia: em vez de \(C_0^+=C_0^-=0\), inicia-se com um valor não nulo, tipicamente \(H/2\) — headstart de 50%

\[C_0^+ = C_0^- = \frac{H}{2}\]

• Se o processo já começa sob controle, o CUSUM cai rapidamente a zero (headstart tem pouco efeito)
• Se o processo já começa fora do alvo, o headstart permite detectar isso muito mais rápido

Implementando o Headstart no R

FS <- H / 2  # headstart de 50%
Cmais_rir <- Cmenos_rir <- numeric(length(x_tab91))
Cmais_rir[1]  <- max(0,  x_tab91[1] - (mu0 + K) + FS)
Cmenos_rir[1] <- max(0, (mu0 - K)  - x_tab91[1] + FS)
for (i in 2:length(x_tab91)) {
  Cmais_rir[i]  <- max(0,  x_tab91[i] - (mu0 + K) + Cmais_rir[i-1])
  Cmenos_rir[i] <- max(0, (mu0 - K)  - x_tab91[i] + Cmenos_rir[i-1])
}

plot_cusum(Cmais_rir, Cmenos_rir, H, titulo = "CUSUM com Headstart de 50% (FS = H/2)")

CUSUM Unilateral

• Há situações em que apenas um dos lados (\(C^+\) ou \(C^-\)) interessa
• Exemplo clássico: viscosidade de um produto químico — quedas abaixo do alvo não são problema, mas qualquer aumento deve ser detectado rapidamente
• Nesse caso, monitora-se apenas \(C^+\), e o ARL correspondente é obtido diretamente pela aproximação de Siegmund unilateral

CUSUM com Sensibilidades Diferentes

• É possível planejar \(k^+ \neq k^-\) (ou \(h^+ \neq h^-\)) quando um desvio (acima ou abaixo do alvo) é mais crítico que o outro
• Útil, por exemplo, quando uma mudança para cima tem consequências de segurança ou custo muito mais graves que uma mudança para baixo

CUSUM para a Variabilidade do Processo

• É possível estender o CUSUM tabular para monitorar \(\sigma\), e não apenas a média (Montgomery, 1981)
• A ideia básica é semelhante: transformar as observações (por exemplo, via \(|x_i - \mu_0|\) ou estatísticas padronizadas de variância) e aplicar a mesma lógica de acumulação unilateral
• Também é possível construir CUSUMs para variáveis Poisson e Binomiais (não conformidades, fração não conforme)

Exercício Aplicado: Viscosidade da Tinta

x_visc <- c(33.75, 33.05, 34.00, 33.81, 33.46, 34.02, 33.68, 33.27,
            33.49, 33.20, 33.62, 33.00, 33.54, 33.12, 33.84)
mu0_v <- 33.5; K_v <- 0.5; H_v <- 4.77
mean(x_visc)

[1] 33.52333

CUSUM Tabular — Viscosidade da Tinta

Cm_v <- Cn_v <- numeric(length(x_visc))
Cm_v[1] <- max(0,  x_visc[1] - (mu0_v + K_v))
Cn_v[1] <- max(0, (mu0_v - K_v) - x_visc[1])
for (i in 2:length(x_visc)) {
  Cm_v[i] <- max(0,  x_visc[i] - (mu0_v + K_v) + Cm_v[i-1])
  Cn_v[i] <- max(0, (mu0_v - K_v) - x_visc[i] + Cn_v[i-1])
}

plot_cusum(Cm_v, Cn_v, H_v, titulo = "CUSUM Tabular — Viscosidade da Tinta")

Viscosidade da Tinta — Zoom em \(C^+\) e \(C^-\)

• Na escala de \(H=4{,}77\), \(C^+\) e \(C^-\) parecem “colados” em zero — mas eles não são identicamente zero em toda a série
• Restringindo o eixo \(y\) a \([-0{,}1,\ 0{,}1]\), aparece um pequeno pico em \(C^+\) na amostra 6 (\(x_6=34{,}02 > \mu_0+K=34{,}0\))

plot_cusum(Cm_v, Cn_v, H_v, titulo = "CUSUM Tabular — Viscosidade da Tinta (zoom)") +
  coord_cartesian(ylim = c(-.1, .1))

• O pico é resetado já na amostra seguinte (\(x_7=33{,}68 < 34{,}0\)) — evidência de ruído normal em torno do alvo, não de uma mudança real no processo

Mesmo Exercício, Agora com Headstart de 50%

FS_v <- H_v / 2
Cm_vh <- Cn_vh <- numeric(length(x_visc))
Cm_vh[1] <- max(0,  x_visc[1] - (mu0_v + K_v) + FS_v)
Cn_vh[1] <- max(0, (mu0_v - K_v) - x_visc[1] + FS_v)
for (i in 2:length(x_visc)) {
  Cm_vh[i] <- max(0,  x_visc[i] - (mu0_v + K_v) + Cm_vh[i-1])
  Cn_vh[i] <- max(0, (mu0_v - K_v) - x_visc[i] + Cn_vh[i-1])
}

plot_cusum(Cm_vh, Cn_vh, H_v, titulo = "CUSUM com Headstart — Viscosidade da Tinta")

• Nenhum dos dois gráficos sinaliza — coerente com um processo estável em torno do alvo

Transição: Do CUSUM ao MMEP

• O MMEP (Média Móvel Exponencialmente Ponderada) é outra alternativa eficaz aos gráficos de Shewhart para pequenas mudanças
• Desempenho equivalente ao CUSUM tabular em muitos cenários
• De certa forma, mais fácil de estabelecer e operar
• Introduzido por S. Roberts (1959)
• Assim como o CUSUM, é tipicamente usado para observações individuais (\(n=1\)), mas se estende a subgrupos racionais

Definição do Gráfico MMEP

\[z_i = \lambda x_i + (1-\lambda) z_{i-1}, \qquad 0 < \lambda \le 1\]

• Valor inicial: \(z_0 = \mu_0\) (ou a média amostral, se o alvo não é conhecido)
• \(\lambda\) controla a “memória” do gráfico: quanto menor \(\lambda\), mais peso para o passado; quanto maior, mais parecido com o gráfico de Shewhart (\(\lambda=1\))

MMEP Como Média Ponderada de Todo o Passado

Substituindo \(z_{i-1}, z_{i-2}, \dots\) recursivamente:

\[z_i = \lambda\sum_{j=0}^{i-1}(1-\lambda)^j x_{i-j} + (1-\lambda)^i z_0\]

• Os pesos \(\lambda(1-\lambda)^j\) decrescem geometricamente com a idade da observação
• Como o MMEP pondera todas as observações passadas e correntes, o gráfico é insensível à hipótese de normalidade — ideal para observações individuais

Variância do MMEP

Se as observações \(x_i\) são independentes com variância \(\sigma^2\):

\[\text{Var}(z_i) = \sigma^2 \left(\frac{\lambda}{2-\lambda}\right)\Bigl[1-(1-\lambda)^{2i}\Bigr]\]

• A variância de \(z_i\) cresce com \(i\) até estabilizar em um valor de regime permanente
• Essa expressão é a base para os limites de controle exatos do MMEP

Limites de Controle Exatos do MMEP

\[LSC_i = \mu_0 + L\,\sigma\sqrt{\left(\frac{\lambda}{2-\lambda}\right)\Bigl[1-(1-\lambda)^{2i}\Bigr]}\]

\[LM = \mu_0\]

\[LIC_i = \mu_0 - L\,\sigma\sqrt{\left(\frac{\lambda}{2-\lambda}\right)\Bigl[1-(1-\lambda)^{2i}\Bigr]}\]

• Os limites variam com \(i\) — mais estreitos no início, alargando até o regime permanente
• \(L\) e \(\lambda\) são os parâmetros de planejamento do gráfico

Limites Assintóticos (Regime Permanente)

Quando \(i \to \infty\), \([1-(1-\lambda)^{2i}] \to 1\), e os limites se aproximam de:

\[LSC = \mu_0 + L\sigma\sqrt{\frac{\lambda}{2-\lambda}} \qquad LIC = \mu_0 - L\sigma\sqrt{\frac{\lambda}{2-\lambda}}\]

• Uso comum na prática por simplicidade — mas Montgomery recomenda enfaticamente usar os limites exatos, especialmente nas primeiras amostras

Exemplo MMEP: Cálculo Manual (1ª Amostra)

\[z_1 = 0{,}1(9{,}45) + (1-0{,}1)(10) = 9{,}945\]

\[LIC_1 = 10 - 2{,}7\sqrt{\tfrac{0,1}{1,9}\bigl[1-(0,9)^2\bigr]} = 9{,}73 \qquad LSC_1 = 10 + 2{,}7\sqrt{\tfrac{0,1}{1,9}\bigl[1-(0,9)^2\bigr]} = 10{,}27\]

Implementando o MMEP no R

lambda <- 0.1; L <- 2.7
z <- numeric(length(x_tab91))
z[1] <- lambda * x_tab91[1] + (1 - lambda) * mu0
for (i in 2:length(x_tab91)) {
  z[i] <- lambda * x_tab91[i] + (1 - lambda) * z[i-1]
}

i_seq <- seq_along(x_tab91)
fator <- sqrt((lambda / (2 - lambda)) * (1 - (1 - lambda)^(2 * i_seq)))
lic_mmep <- mu0 - L * sigma * fator
lsc_mmep <- mu0 + L * sigma * fator
round(cbind(amostra = 1:5, z = z[1:5], LIC = lic_mmep[1:5], LSC = lsc_mmep[1:5]), 3)

     amostra      z   LIC    LSC
[1,]       1  9.945 9.730 10.270
[2,]       2  9.749 9.637 10.363
[3,]       3  9.704 9.576 10.424
[4,]       4  9.899 9.533 10.467
[5,]       5 10.125 9.500 10.500

Gráfico MMEP Resultante — Tabela 9.1

plot_mmep(z, mu0, lic_mmep, lsc_mmep,
          titulo = "Gráfico MMEP (lambda = 0,1 ; L = 2,7) — Tabela 9.1")

• Assim como o CUSUM, o MMEP sinaliza a mudança de média perto do final da série — reagindo mais rápido que o Shewhart, mas com uma resposta mais suave (efeito de “alisamento”)

Confirmando com o Pacote qcc

q_ewma <- ewma(x_tab91, center = mu0, std.dev = sigma,
               lambda = 0.1, nsigmas = 2.7, plot = FALSE)
summary(q_ewma)

Call:
ewma(data = x_tab91, center = mu0, std.dev = sigma, lambda = 0.1,     nsigmas = 2.7, plot = FALSE)

ewma chart for x_tab91 

Summary of group statistics:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 7.9900  9.3775 10.4450 10.3150 11.2525 12.2900 

Group sample size:  1
Number of groups:  30
Center of group statistics:  10
Standard deviation:  1 

Smoothing parameter: 0.1 
Control limits:
           LCL      UCL
[1,]  9.730000 10.27000
[2,]  9.636752 10.36325
...                    
[30,] 9.381134 10.61887

Planejamento do Gráfico MMEP

• Os parâmetros de planejamento são \(L\) (multiplicador dos limites) e \(\lambda\) (memória)
• Podem ser escolhidos para obter um desempenho de ARL comparável ao do CUSUM
• Estudo de referência: Lucas & Saccucci (1990)

ARL do MMEP — Lucas & Saccucci (1990)

tab_ls <- data.frame(
  delta = c(0, 0.25, 0.5, 0.75, 1.00, 1.50, 2.00, 2.50, 3.00, 4.00),
  `lambda_04`  = c(500.0, 224.0, 71.2, 28.4, 14.3, 5.9, 3.5, 2.5, 2.0, 1.4),
  `lambda_025` = c(500.0, 170.0, 48.2, 20.1, 11.1, 5.5, 3.6, 2.7, 2.3, 1.7),
  `lambda_02`  = c(500.0, 150.0, 41.8, 18.2, 10.5, 5.5, 3.7, 2.9, 2.4, 1.9),
  `lambda_01`  = c(500.0, 106.0, 31.3, 15.9, 10.3, 6.1, 4.4, 3.4, 2.9, 2.2),
  `lambda_005` = c(500.0, 84.1, 28.2, 16.4, 11.4, 7.1, 5.2, 4.2, 3.5, 2.7),
  check.names = FALSE
)
tab_ls

   delta lambda_04 lambda_025 lambda_02 lambda_01 lambda_005
1   0.00     500.0      500.0     500.0     500.0      500.0
2   0.25     224.0      170.0     150.0     106.0       84.1
3   0.50      71.2       48.2      41.8      31.3       28.2
4   0.75      28.4       20.1      18.2      15.9       16.4
5   1.00      14.3       11.1      10.5      10.3       11.4
6   1.50       5.9        5.5       5.5       6.1        7.1
7   2.00       3.5        3.6       3.7       4.4        5.2
8   2.50       2.5        2.7       2.9       3.4        4.2
9   3.00       2.0        2.3       2.4       2.9        3.5
10  4.00       1.4        1.7       1.9       2.2        2.7

Visualizando o ARL do MMEP

df_ls <- data.frame(
  delta  = rep(tab_ls$delta, 5),
  arl    = c(tab_ls$lambda_04, tab_ls$lambda_025, tab_ls$lambda_02,
             tab_ls$lambda_01, tab_ls$lambda_005),
  lambda = rep(c("0,40", "0,25", "0,20", "0,10", "0,05"), each = nrow(tab_ls))
)
ggplot(df_ls, aes(delta, arl, color = lambda, group = lambda)) +
  geom_line(linewidth = .9) + geom_point(size = 1.6) +
  scale_y_log10() +
  scale_color_manual(values = unname(pal[c("azul","verde","ambar","roxo","coral")])) +
  labs(title = "ARL do MMEP para Varios Valores de lambda",
       x = expression(delta~"(mudanca em unidades de sigma)"),
       y = "ARL (escala log)", color = expression(lambda)) +
  tema_cep()

Conclusões do Estudo de Lucas & Saccucci

• Valores \(0{,}05 \le \lambda \le 0{,}25\) funcionam bem na prática; \(\lambda \in \{0{,}05;\ 0{,}1;\ 0{,}2\}\) são escolhas populares
• Usar \(\lambda\) menor para detectar mudanças pequenas
• \(L=3\) funciona razoavelmente bem com \(\lambda > 0{,}25\)
• Quando \(\lambda \le 0{,}10\), deve-se trabalhar com \(2{,}6 \le L \le 2{,}7\)

MMEP–Shewhart Combinado para Grandes Mudanças

• Assim como no CUSUM, o MMEP tem desempenho fraco contra mudanças grandes
• Solução: procedimento combinado MMEP–Shewhart (Lucas, 1982)
• Construir um gráfico de Shewhart sobre \(z_i\), com limites de \(3{,}25\sigma\) ou \(3{,}5\sigma\)
• O MMEP cuida das mudanças pequenas; o Shewhart, das grandes; sinal em qualquer um dos dois aciona investigação

Robustez do MMEP à Não Normalidade

• O gráfico de Shewhart para observações individuais é muito sensível à não normalidade — gera excesso de alarmes falsos
• Borror, Montgomery & Runger (1999) compararam o \(ARL_0\) de Shewhart vs. MMEP usando:
  • Distribuição Gama (assimétrica)
  • Distribuição t-Student (simétrica, caudas pesadas)

\(ARL_0\) sob Distribuições Assimétricas (Gama)

tab_gama <- data.frame(
  distribuicao = c("Normal", "Gama(4,1)", "Gama(3,1)", "Gama(2,1)", "Gama(1,1)", "Gama(0,5;1)"),
  `lambda_005` = c(370, 372, 372, 372, 369, 357),
  `lambda_01`  = c(371, 341, 332, 315, 274, 229),
  `lambda_02`  = c(371, 259, 238, 208, 163, 131),
  `Shewhart_lambda_1` = c(370, 97, 85, 71, 55, 45),
  check.names = FALSE
)
tab_gama

  distribuicao lambda_005 lambda_01 lambda_02 Shewhart_lambda_1
1       Normal        370       371       371               370
2    Gama(4,1)        372       341       259                97
3    Gama(3,1)        372       332       238                85
4    Gama(2,1)        372       315       208                71
5    Gama(1,1)        369       274       163                55
6  Gama(0,5;1)        357       229       131                45

\(ARL_0\) sob Distribuições Simétricas (t-Student)

tab_t <- data.frame(
  distribuicao = c("Normal", "t(50)", "t(40)", "t(30)", "t(20)", "t(15)", "t(10)", "t(8)", "t(6)", "t(4)"),
  `lambda_005` = c(370, 369, 369, 368, 367, 365, 361, 358, 351, 343),
  `lambda_01`  = c(371, 365, 363, 361, 355, 349, 335, 324, 305, 274),
  `lambda_02`  = c(371, 353, 348, 341, 325, 310, 280, 259, 229, 188),
  `Shewhart_lambda_1` = c(370, 283, 266, 242, 204, 176, 137, 117, 96, 76),
  check.names = FALSE
)
tab_t

   distribuicao lambda_005 lambda_01 lambda_02 Shewhart_lambda_1
1        Normal        370       371       371               370
2         t(50)        369       365       353               283
3         t(40)        369       363       348               266
4         t(30)        368       361       341               242
5         t(20)        367       355       325               204
6         t(15)        365       349       310               176
7         t(10)        361       335       280               137
8          t(8)        358       324       259               117
9          t(6)        351       305       229                96
10         t(4)        343       274       188                76

Visualizando a Robustez: MMEP vs. Shewhart

df_rob <- data.frame(
  gl  = c(50, 40, 30, 20, 15, 10, 8, 6, 4),
  mmep_005 = tab_t$lambda_005[-1],
  shewhart = tab_t$Shewhart_lambda_1[-1]
)
df_rob_long <- data.frame(
  gl = rep(df_rob$gl, 2),
  arl0 = c(df_rob$mmep_005, df_rob$shewhart),
  metodo = rep(c("MMEP (lambda=0,05)", "Shewhart (individual)"), each = nrow(df_rob))
)
ggplot(df_rob_long, aes(gl, arl0, color = metodo)) +
  geom_line(linewidth = 1) + geom_point(size = 2) +
  geom_hline(yintercept = 370, linetype = "dotted", color = pal["neutro"]) +
  scale_color_manual(values = c(unname(pal["azul"]), unname(pal["coral"]))) +
  labs(title = "ARL0 sob t-Student: Queda mais suave no MMEP",
       x = "Graus de liberdade da t-Student", y = expression(ARL[0]),
       color = NULL) +
  tema_cep()

Conclusões sobre Robustez

• Distribuições não normais reduzem drasticamente o \(ARL_0\) do Shewhart para observações individuais → excesso de alarmes falsos
• Um MMEP bem planejado (com \(\lambda\) pequeno) mantém desempenho satisfatório tanto sob normalidade quanto sob não normalidade
• Recomendação prática: priorizar o MMEP bem planejado como gráfico de controle padrão para medidas individuais

Subgrupos Racionais no MMEP

• Extensão imediata a \(n>1\): substituir \(x_i\) por \(\bar x_i\) e \(\sigma\) por \(\sigma_{\bar x} = \sigma/\sqrt{n}\)
• As mesmas considerações de eficiência do CUSUM se aplicam — subgrupos maiores nem sempre compensam a perda de frequência de amostragem

Extensões do MMEP

• MMEP multivariado (Lowry et al., 1992) — monitoramento simultâneo de vários parâmetros
• MMEP para monitorar variabilidade do processo, e não apenas a média
• MMEP com limites variáveis no tempo e resposta inicial rápida (Steiner, 1999)
• Aplicações em processos autocorrelacionados — tema do Capítulo 10

Seção 9.3: O Gráfico da Média Móvel (Não Ponderada)

\[\bar M_i = \frac{x_i + x_{i-1} + \cdots + x_{i-w+1}}{w}\]

• Média móvel simples, com janela de largura \(w\) — todas as observações na janela recebem o mesmo peso
• Diferente do MMEP, cujos pesos decrescem geometricamente
• Menos usado na prática que CUSUM e MMEP, mas conceitualmente mais simples de explicar

Limites de Controle da Média Móvel

\[LSC = \mu_0 + 3\frac{\sigma}{\sqrt{w}} \qquad LM = \mu_0 \qquad LIC = \mu_0 - 3\frac{\sigma}{\sqrt{w}}\]

• Nas primeiras \(w-1\) amostras, a variância de \(\bar M_i\) é maior (janela ainda incompleta) — os limites, na prática exata, se ajustam a cada \(i < w\)
• Assim como no MMEP, \(\bar M_i\) tende a suavizar ruído — mas sem a robustez à não normalidade que o MMEP oferece

Média Móvel Aplicada à Tabela 9.1

w <- 5
M <- stats::filter(x_tab91, rep(1/w, w), sides = 1)
lic_mm <- mu0 - 3 * sigma / sqrt(w)
lsc_mm <- mu0 + 3 * sigma / sqrt(w)

df_mm <- data.frame(t = seq_along(x_tab91), M = as.numeric(M))
df_mm <- df_mm[!is.na(df_mm$M), ]
plot_carta(df_mm$M, center = mu0, lic = lic_mm, lsc = lsc_mm,
           titulo = "Gráfico da Média Móvel (w=5) — Tabela 9.1", ylab = expression(bar(M)[i]))

Comparando os Três Gráficos: Shewhart, CUSUM e MMEP

p1 <- plot_carta(x_tab91, mu0, mu0 - 3*sigma, mu0 + 3*sigma, "Shewhart")
p2 <- plot_cusum(Cmais, Cmenos, H, "CUSUM Tabular")
p3 <- plot_mmep(z, mu0, lic_mmep, lsc_mmep, "MMEP")
p1 / p2 / p3

Leitura Comparativa

• O Shewhart não sinaliza a mudança de \(1\sigma\)
• O CUSUM sinaliza claramente a partir da amostra 26–27
• O MMEP também sinaliza, com trajetória mais suave — resposta ligeiramente mais lenta que o CUSUM neste exemplo específico, mas de magnitude comparável
• Em geral, CUSUM e MMEP bem planejados têm desempenho muito próximo

Exercício Aplicado 1: Enchimento de Latas de Óleo

x_oleo <- c(8.00, 8.01, 8.02, 8.01, 8.00, 8.01, 8.06, 8.07, 8.01, 8.04, 8.02, 8.01,
            8.05, 8.04, 8.03, 8.05, 8.06, 8.04, 8.05, 8.06, 8.04, 8.02, 8.03, 8.05)
sigma_oleo <- 0.05
mu0_oleo <- 8.02  # alvo do processo

(a) CUSUM Tabular — Latas de Óleo

K_o <- 0.5 * sigma_oleo; H_o <- 4.77 * sigma_oleo
Cm_o <- Cn_o <- numeric(length(x_oleo))
Cm_o[1] <- max(0,  x_oleo[1] - (mu0_oleo + K_o))
Cn_o[1] <- max(0, (mu0_oleo - K_o) - x_oleo[1])
for (i in 2:length(x_oleo)) {
  Cm_o[i] <- max(0,  x_oleo[i] - (mu0_oleo + K_o) + Cm_o[i-1])
  Cn_o[i] <- max(0, (mu0_oleo - K_o) - x_oleo[i] + Cn_o[i-1])
}

plot_cusum(Cm_o, Cn_o, H_o, "CUSUM Tabular — Latas de Óleo (mu0=8,02)")

(b) CUSUM RIR — Ajuste ao Alvo \(\mu=8{,}0\)

mu0_ajuste <- 8.0
FS_o <- H_o / 2
Cm_or <- Cn_or <- numeric(length(x_oleo))
Cm_or[1] <- max(0,  x_oleo[1] - (mu0_ajuste + K_o) + FS_o)
Cn_or[1] <- max(0, (mu0_ajuste - K_o) - x_oleo[1] + FS_o)
for (i in 2:length(x_oleo)) {
  Cm_or[i] <- max(0,  x_oleo[i] - (mu0_ajuste + K_o) + Cm_or[i-1])
  Cn_or[i] <- max(0, (mu0_ajuste - K_o) - x_oleo[i] + Cn_or[i-1])
}

plot_cusum(Cm_or, Cn_or, H_o, "CUSUM RIR (headstart 50%) — Alvo mu = 8,0")

Interpretando o CUSUM RIR

• Com o headstart, \(C^+\) já parte de \(H_o/2\)
• Se as primeiras observações confirmam que o processo retornou ao alvo, \(C^+\) cai rapidamente para zero
• Se o ajuste não foi bem-sucedido (processo ainda deslocado), o sinal aparece muito mais cedo do que sem o headstart

(c) Gráfico MMEP — Latas de Óleo

lambda_o <- 0.2; L_o <- 3
z_o <- numeric(length(x_oleo))
z_o[1] <- lambda_o * x_oleo[1] + (1 - lambda_o) * mu0_oleo
for (i in 2:length(x_oleo)) {
  z_o[i] <- lambda_o * x_oleo[i] + (1 - lambda_o) * z_o[i-1]
}
i_seq_o <- seq_along(x_oleo)
fator_o <- sqrt((lambda_o / (2 - lambda_o)) * (1 - (1 - lambda_o)^(2 * i_seq_o)))
lic_o <- mu0_oleo - L_o * sigma_oleo * fator_o
lsc_o <- mu0_oleo + L_o * sigma_oleo * fator_o

plot_mmep(z_o, mu0_oleo, lic_o, lsc_o, "MMEP (lambda=0,2 ; L=3) — Latas de Óleo")

Comparando os Três Gráficos — Latas de Óleo

q1 <- plot_cusum(Cm_o, Cn_o, H_o, "(a) CUSUM Tabular")
q2 <- plot_cusum(Cm_or, Cn_or, H_o, "(b) CUSUM RIR (alvo ajustado)")
q3 <- plot_mmep(z_o, mu0_oleo, lic_o, lsc_o, "(c) MMEP")
q1 / q2 / q3

• Nenhum dos três gráficos aponta um processo fora de controle — consistente com um processo estável e bem ajustado ao alvo

Exercício Aplicado 2: Simulação Comparativa

set.seed(13)
x_sim <- c(rnorm(20, 10, 1), rnorm(20, 9, 1))
mu0_sim <- 10; sigma_sim <- 1

Simulação: Aplicando o CUSUM

K_s <- 0.5 * sigma_sim; H_s <- 4.774 * sigma_sim
Cm_s <- Cn_s <- numeric(length(x_sim))
Cm_s[1] <- max(0,  x_sim[1] - (mu0_sim + K_s))
Cn_s[1] <- max(0, (mu0_sim - K_s) - x_sim[1])
for (i in 2:length(x_sim)) {
  Cm_s[i] <- max(0,  x_sim[i] - (mu0_sim + K_s) + Cm_s[i-1])
  Cn_s[i] <- max(0, (mu0_sim - K_s) - x_sim[i] + Cn_s[i-1])
}
sinal_cusum <- which(Cm_s > H_s | Cn_s > H_s)[1]
sinal_cusum

[1] 28

Simulação: Aplicando o MMEP

lambda_s <- 0.20; L_s <- 2.859
z_s <- numeric(length(x_sim))
z_s[1] <- lambda_s * x_sim[1] + (1 - lambda_s) * mu0_sim
for (i in 2:length(x_sim)) {
  z_s[i] <- lambda_s * x_sim[i] + (1 - lambda_s) * z_s[i-1]
}
i_seq_s <- seq_along(x_sim)
fator_s <- sqrt((lambda_s / (2 - lambda_s)) * (1 - (1 - lambda_s)^(2 * i_seq_s)))
lic_s <- mu0_sim - L_s * sigma_sim * fator_s
lsc_s <- mu0_sim + L_s * sigma_sim * fator_s
sinal_mmep <- which(z_s > lsc_s | z_s < lic_s)[1]
sinal_mmep

[1] 28

Comparando os Dois Sinais

cat(sprintf("CUSUM sinalizou na amostra: %s\nMMEP sinalizou na amostra: %s\n",
            ifelse(is.na(sinal_cusum), "nao sinalizou", sinal_cusum),
            ifelse(is.na(sinal_mmep), "nao sinalizou", sinal_mmep)))

CUSUM sinalizou na amostra: 28
MMEP sinalizou na amostra: 28

r1 <- plot_cusum(Cm_s, Cn_s, H_s, "CUSUM (h=4,774 ; k=0,5)")
r2 <- plot_mmep(z_s, mu0_sim, lic_s, lsc_s, "MMEP (lambda=0,20 ; L=2,859)")
r1 / r2

Interpretando a Simulação

• Ambos os métodos foram planejados para o mesmo \(ARL_0 \approx 370\) — a comparação é justa
• Com parâmetros equivalentes, CUSUM e MMEP tendem a sinalizar em amostras muito próximas
• Pequenas diferenças de uma simulação para outra são esperadas — o ideal é comparar o desempenho médio via Monte Carlo (muitas replicações), não uma única simulação

Resumo de Fórmulas — Capítulo 9

Comparando CUSUM e MMEP: Quando Usar Cada Um

• CUSUM: mais natural quando o interesse é detectar uma mudança de magnitude específica e conhecida (\(\delta\))
• MMEP: mais flexível para diferentes magnitudes de mudança; mais robusto à não normalidade
• Ambos superam o Shewhart isolado para mudanças pequenas; ambos perdem para o Shewhart em mudanças muito grandes (daí os esquemas combinados)
• Na prática, a escolha frequentemente se resume à familiaridade da equipe e à disponibilidade de suporte (ex.: pacote qcc, Minitab)

Glossário Rápido de Siglas

Conclusão

• CUSUM e MMEP incorporam toda a sequência histórica de observações — muito mais sensíveis a pequenas mudanças que o gráfico de Shewhart
• O planejamento de ambos gira em torno do ARL: escolher parâmetros que equilibrem falsos alarmes (\(ARL_0\) alto) e detecção rápida (\(ARL_1\) baixo)
• O MMEP se destaca pela robustez à não normalidade, sendo recomendado como padrão para observações individuais
• Esquemas combinados com Shewhart cobrem tanto mudanças pequenas quanto grandes
• O pacote qcc no R implementa cusum() e ewma() prontos para uso em produção