ewma_qcc <- function(amostra_inicial,
lambda,
amostra_subsequente = NULL,
nsigmas = 1.96) {
arguments <- list(
amostra_inicial,
lambda = lambda,
nsigmas = nsigmas,
plot = F
)
if (!is.null(amostra_subsequente)) {
arguments$newdata <- amostra_subsequente
}
ew <- do.call(ewma, arguments)
registros <- if (is.null(amostra_subsequente)) {
seq(1, length(amostra_inicial))
} else {
seq(length(amostra_inicial) + 1, length(amostra_inicial) + length(amostra_subsequente))
}
ewma <- as.numeric(ew$y[registros])
LI <- ew$limits[, 1][registros]
LS <- ew$limits[, 2][registros]
fora_de_controle <- ewma < LI | ewma > LS
total_fora_de_controle <- sum(fora_de_controle, na.rm = TRUE)
fracao_fora_de_controle <- total_fora_de_controle / length(ewma)
list(
ewma = ewma,
LI = LI,
LS = LS,
fora_de_controle = fora_de_controle,
total_fora_de_controle = total_fora_de_controle,
fracao_fora_de_controle = fracao_fora_de_controle
)
}2 CEQ usando resíduos
2.1 EWMA
O EWMA (Exponential Weighted Moving Average – Média Móvel Exponencialmente Ponderada) é um método de controle de processos que utiliza uma média móvel ponderada exponencialmente para detectar mudanças no processo (MONTGOMERY, 2009)5.
Ele parte de uma ideia simples: ao invés de usarmos apenas a última amostra, usamos uma média ponderada de todas as amostras anteriores. A ponderação é exponencial, o que significa que amostras mais antigas têm menos peso que amostras mais recentes.
Matematicamente, a média EWMA é dada por:
\[ z_i = \lambda y_i + (1 - \lambda) z_{i-1} \]
Onde, \(z_i\) é a estatística de controle no instante \(i\), \(y_i\) é a observação no instante \(i\), \(\lambda\) é o fator de suavização e \(z_{i-1}\) é a estatística de controle no instante anterior. O valor inicial de \(z_0\) é definido como a média do processo, tal que \(z_0 = \mu_0\).
Outro fato importante, é que, uma vez que o EWMA é uma média ponderada de todas as amostras anteriores, ele é pouco sensível à suposição de normalidade dos dados.
A estatística, \(z_i\), é, então, comparada com os limites de controle, \(\text{LCS}\) e \(\text{LCI}\), definidos como:
\[ \text{LC} = \mu_0 \pm L \sigma \sqrt{\frac{\lambda}{2 - \lambda} \left[1 - \left(1 - \lambda\right)^{2i}\right]} \]
Aqui, utilizaremos o pacote qcc4 para implementar o EWMA.
2.1.1 O fator λ
O \(\lambda\) é uma constante definida no intervalo \((0, 1]\), quanto mais próximo de 1, mais peso é dado à amostra mais recente, tanto que, quando \(\lambda = 1\), teremos a carta de controle de Shewhart, pois a média EWMA será igual à média das amostras.
Para Montgomery (2009)5, em geral, valores de \(\lambda\) entre 0.05 e 0.25 são recomendados. No entanto, esta escolha depende do tipo de processo e do grau de sensibilidade desejado.
2.1.2 Simulação
Vamos buscar o melhor valor de \(\lambda\) para o EWMA, para um nível de significância constante, \(\alpha = 0.05\).
ewma_monte_carlo <- cache_dados("simulacao-ewma", function() {
gerador_monte_carlo(parametros = list(lambda = seq(0.1, 0.9, by = 0.1))) %>%
mutate(
qcc = list(ewma_qcc(f1_residuos, lambda, f2_residuos)),
ewma = list(qcc$ewma),
LI = list(qcc$LI),
LS = list(qcc$LS),
fora_de_controle = list(qcc$fora_de_controle),
total_fora_de_controle = qcc$total_fora_de_controle,
fracao_fora_de_controle = qcc$fracao_fora_de_controle
) %>%
select(-qcc)
})
ewma_monte_carlo$lambda <- as.factor(ewma_monte_carlo$lambda)2.1.2.1 Resumo
Resumo da fração de pontos fora de controle (FPFC) para diferentes valores de \(\lambda\) e \(\Phi_2\).
2.1.2.2 Cartas para λ = 0.8
Vamos analisar o EWMA para \(\lambda = 0.8\). Aqui, vamos considerar apenas a primeira execução.
2.2 EWMA-AR
Nota do autor: O EWMA-AR é uma extensão do EWMA que utiliza um modelo autorregressivo para prever a próxima observação. Como será apresentado a seguir, o EWMA-AR possuim um poder muito menor que o EWMA para detectar mudanças no processo. Por consequência disso, o EWMA-AR é tratado aqui apenas como uma curiosidade.
Segundo Montgomery (2009)5, o EWMA-AR é uma extensão do EWMA que utiliza um modelo autorregressivo para prever a próxima observação.
Assim, temos que \(\lambda \in (0, 1]\), sendo que, a previsão para a observação \(x_{t+1}\) é dada por \(\hat{x}_{t+1}(t)=z_{t} = \lambda x_t + (1 - \lambda) z_{t-1}\).
2.2.1 Otimizando λ
De acordo com Montgomery (2009)5, podemos encontrar um valor ótimo para \(\lambda\) através da minimização da soma dos quadrados dos resíduos.
E, ainda, temos que os erros de previsão são dados por \(e_t = x_t - \hat{x}_t(t-1)\) conforme a Eq. 10.16 (Montgomery, 2009)5.
Ou, de outra forma
\[ \hat{\lambda} = \min \left( \underset{\lambda\in(0,1]}{\arg\min}\,\text{Err}(\lambda) \right) \]
onde \(\hat{\lambda}\) é o \(\lambda\) ótimo por simulação, e \(\text{Err}(\lambda) = \sum_{t=1}^{n} e_t^2\).
Segundo, também Montgomery (2009)5, podemos estimar o valor de \(\sigma^2\) para o modelo EWMA-AR como \(\sigma^2 = \frac{\sum (\text{err}_i^2|_\lambda)}{n}\). Onde \(\text{err}|_\lambda\) são os resíduos do modelo EWMA-AR para o melhor valor de \(\lambda\) encontrado.
# Função para computar resíduos do EWMA-AR
ewma_ar_residuos <- function(dados, ewma) {
n <- length(dados)
residuos <- numeric(n)
residuos[1] <- dados[1]
for (i in 2:n) {
residuos[i] <- dados[i] - ewma[i - 1]
}
return(residuos)
}
# Função para computar o EWMA-AR
ewma_ar <- function(dados,
lambda,
x0 = NULL,
desvio = NULL) {
desv <- ifelse(is.null(desvio), sqrt(sd(dados)), desvio)
n <- length(dados)
final <- ifelse(is.null(x0), n, n + 1)
ewma_serie <- numeric(final)
ewma_serie[1] <- ifelse(is.null(x0), dados[1], x0)
for (i in 2:final) {
ewma_serie[i] <- lambda * dados[i] + (1 - lambda) * ewma_serie[i - 1]
}
ewma <- ewma_serie[ifelse(is.null(x0), 1, 2):final]
LI <- ewma - 1.96 * desv
LS <- ewma + 1.96 * desv
fora_de_controle <- dados < LI | dados > LS
total_fora_de_controle <- sum(fora_de_controle, na.rm = TRUE)
fracao_fora_de_controle <- total_fora_de_controle / length(ewma)
return(
list(
ewma = ewma,
residuos = ewma_ar_residuos(dados, ewma_serie),
LI = LI,
LS = LS,
fora_de_controle = fora_de_controle,
total_fora_de_controle = total_fora_de_controle,
fracao_fora_de_controle = fracao_fora_de_controle
)
)
}lambda_otimo_optimize <- function(amostra_inicial) {
lambda <- NA
minimo <- Inf
for (x in c(seq(0.01, 0.1, by = 0.01), seq(0.1, 0.9, by = 0.1))) {
ew <- ewma_ar(amostra_inicial, lambda = x)
erros <- ew$residuos
soma_quadrados <- sum(erros^2)
if (soma_quadrados < minimo) {
minimo <- soma_quadrados
lambda <- x
}
}
return(list(lambda = lambda, minimo = minimo))
}
ewmar_lambda_otimo <- cache_dados("ewmaar-lambda-otimo", \() {
data.frame(id = 1:100) %>%
rowwise() %>%
mutate(
f1_amostra = list(barma.sim(n = n1, phi = phi_parametro)),
f1_phi = list(barma.phi_estimado(f1_amostra)),
f1_residuos = list(barma.residuos(f1_amostra, f1_phi)),
ewma_ar = list(lambda_otimo_optimize(f1_residuos)),
lambda = ewma_ar$lambda,
minimo = ewma_ar$minimo
) %>%
select(-ewma_ar)
})Por simulação encontramos um \(\lambda\) ótimo de \(\sim 0.07\). Com erro médio de \(\sim 22\), o que nos dá um desvio de \(22 / 100 \simeq 0.2\).
2.2.2 Simulação
ewmaar_monte_carlo <- cache_dados("simulacao-ewma-ar", function() {
gerador_monte_carlo(parametros = list(lambda = c(0.1, 0.07, 0.03))) %>%
mutate(
ewma_ar = list(ewma_ar(
f2_controle, lambda, f1_amostras[n1], sd(f1_amostras)
)),
ewma = list(ewma_ar$ewma),
LI = list(ewma_ar$LI),
LS = list(ewma_ar$LS),
fora_de_controle = list(ewma_ar$fora_de_controle),
total_fora_de_controle = ewma_ar$total_fora_de_controle,
fracao_fora_de_controle = ewma_ar$fracao_fora_de_controle
) %>%
select(-ewma_ar)
})
ewmaar_monte_carlo$lambda <- as.factor(ewmaar_monte_carlo$lambda)
ewmaar_monte_carlo_resumo <- ewmaar_monte_carlo %>%
group_by(lambda, f2_phi) %>%
summarise(
mean = mean(fracao_fora_de_controle),
min = min(fracao_fora_de_controle),
max = max(fracao_fora_de_controle),
.groups = "drop"
)
datatable(
ewmaar_monte_carlo_resumo %>%
mutate(across(where(is.numeric), \(x) round(x, digits = 4))),
caption = "Pontos fora de controle",
colnames = c(
"Valores de λ", "Valores de Φ₂", "Média", "Mínimo", "Máximo"
)
)Vamos analisar o EWMA-AR para \(\lambda = 0.07\). O resultado encontrado anteriormente.
Vamos, além disso, comparar com o EWMA realizado a partir dos resíduos do modelo βAR(1).
2.3 CUSUM
O CUSUM (Cumulative Sum – Soma Acumulada) é um método de controle de processos que utiliza a soma acumulada dos resíduos para detectar mudanças no processo (MONTGOMERY, 2009)5.
Neste método são definidas duas constantes, \(H\) e \(K\), que representam o tamanho da mudança que se deseja detectar e a sensibilidade do método, respectivamente. As estatísticas de controle, sendo duas, são definidas como:
\[ \begin{matrix} \text{C}^+_i & = & \max\left[0, x_i - (\mu_0 + K) + C^+_{i-1}\right] \\ \text{C}^-_i & = & \max\left[0, (\mu_0 - K) - x_i + C^-_{i-1}\right] \\ \end{matrix} \]
onde \(x_i\) é a observação no instante \(i\), \(\mu_0\) é a média do processo, \(C^+_0 = C^-_0 = 0\) e \(i = 1, 2, \ldots, n\).
A constante \(K\), chamada de valor de referência, geralmente, segundo MONTGOMERY (2009)5, é definida como a metade entre o \(\mu_0\) alvo e o valor fora de controle de média \(\mu_1\) que queremos detectar.
Já a constante \(H\), chamada de limite de decisão, é definida como o valor que a estatística de controle deve atingir para podermos detectar a mudança no processo. Para MONTGOMERY (2009)5, um valor razoável para \(H\) é \(5\sigma\). Assim, se \(\text{C}^+_i \geq H\) ou \(\text{C}^-_i \leq -H\), então o processo está fora de controle.
Aqui, utilizaremos o pacote qcc4 para implementar o CUSUM.
cusum_qcc <- function(amostra_inicial,
desvio_detectavel = 1,
intervalo_de_decisao = 5,
amostra_subsequente = NULL) {
arguments <- list(
data = amostra_inicial,
se.shift = desvio_detectavel,
decision.interval = intervalo_de_decisao,
plot = F
)
if (!is.null(amostra_subsequente)) {
arguments$newdata <- amostra_subsequente
}
cu <- do.call(cusum, arguments)
# registros <- if (is.null(amostra_subsequente)) {
# seq(1, nH0)
# } else {
# seq(nH0 + 1, nH0 + nH1)
# }
registros <- if (is.null(amostra_subsequente)) {
seq(1, length(amostra_inicial))
} else {
seq(length(amostra_inicial) + 1, length(amostra_inicial) + length(amostra_subsequente))
}
pos <- cu[["pos"]][registros]
neg <- cu[["neg"]][registros]
fora_de_controle_pos <- pos > intervalo_de_decisao
fora_de_controle_neg <- neg < -intervalo_de_decisao
fora_de_controle <- ifelse(
fora_de_controle_pos,
fora_de_controle_pos,
fora_de_controle_neg
)
total_fora_de_controle <- sum(fora_de_controle)
fracao_fora_de_controle <- total_fora_de_controle / length(registros)
list(
pos = pos,
neg = neg,
fora_de_controle_pos = fora_de_controle_pos,
fora_de_controle_neg = fora_de_controle_neg,
fora_de_controle = fora_de_controle,
total_fora_de_controle = total_fora_de_controle,
fracao_fora_de_controle = fracao_fora_de_controle
)
}2.3.1 Simulação
Vamos buscar o melhor valor de \(K\) e \(H\) para o CUSUM, para um nível de significância constante.
comb_cusum <- cache_dados("simulacao-cusum-k-h", function() {
gerador_monte_carlo(parametros = expand.grid(
desvio_detectavel = seq(0.6, 1.8, by = 0.2),
intervalo_de_decisao = seq(3, 6, by = 1)
)) %>%
mutate(
qcc = list(
cusum_qcc(
f1_residuos,
desvio_detectavel = desvio_detectavel,
intervalo_de_decisao = intervalo_de_decisao,
amostra_subsequente = f2_residuos
)
),
pos = list(qcc$pos),
neg = list(qcc$neg),
fora_de_controle_pos = list(qcc$fora_de_controle_pos),
fora_de_controle_neg = list(qcc$fora_de_controle_neg),
fora_de_controle = list(qcc$fora_de_controle),
total_fora_de_controle = qcc$total_fora_de_controle,
fracao_fora_de_controle = qcc$fracao_fora_de_controle
) %>%
select(-qcc)
})
comb_cusum$parametro <- as.factor(paste0(comb_cusum$desvio_detectavel, ";", comb_cusum$intervalo_de_decisao))2.3.1.1 Resumo
Resumo da fração de pontos fora de controle (FPFC) para diferentes valores de \(K\) e \(H\).
2.3.1.2 Cartas para K = 1 e H = 4
2.4 EWMA vs. CUSUM
Vamos comparar os melhores valores de \(\lambda\) e \(K/H\) para o EWMA e CUSUM, respectivamente.
Observa-se que o CUSUM possui mais poder de detecção que o EWMA, mas, em compensação, possui uma variabilidade maior.
2.4.1 ARL
O ARL (Average Run Length – Comprimento Médio de Execução) é uma medida de desempenho de um método de controle de processo. É definido como o número médio de observações necessárias para detectar uma mudança no processo. E pode ser calculado como
\[ \text{ARL} = \frac{1}{\text{FPR}} \]
onde FPR é a taxa de falsos positivos, ou seja, a fração de vezes que o método detecta uma mudança no processo quando não há. Ou seja, é a taxa de pontos fora de controle quando o processo está sob controle.
Vamos ver como o ARL se comporta para os melhores valores de \(\lambda\) e \(K/H\) encontrados.
arl <- function(fpr) {
return(1 / fpr)
}
arl <- dados_resumo_ %>%
filter(f2_phi == 0.2) %>%
group_by(parametro, algoritmo) %>%
summarise(
mean = mean(mean),
min = min(min),
max = max(max),
.groups = "drop"
) %>%
mutate(
arl = ceiling(1 / mean)
)
datatable(
arl %>%
mutate(across(where(is.numeric), \(x) round(x, digits = 4))),
caption = "ARL",
colnames = c(
"Parâmetro",
"Algoritmo",
"Média",
"Mínimo",
"Máximo",
"ARL"
)
)2.5 Referências
1. ROCHA, A. V., CRIBARI-NETO, F. Beta autoregressive moving average models. TEST 18, 529–545 (2009). DOI: 10.1007/s11749-008-0112-z
2. ROCHA, A. V., CRIBARI-NETO, F. Erratum to: Beta autoregressive moving average models. TEST 26, 451–459 (2017). DOI: 10.1007/s11749-017-0528-4
3. PRASS, T. S., et. al. BTSR: Bounded Time Series Regression. R package version 0.1.5. 2023-09-22. DOI: 10.32614/CRAN.package.BTSR
4. SCRUCCA, L., et. al. qcc: Quality Control Charts. R package version 2.7. 2017-07-09. DOI: 10.32614/CRAN.package.qcc
5. MONTGOMERY, D. C. Introduction to Statistical Quality Control. 2013. John Wiley & Sons.