5  Otros Métodos Estadísticos de Control de Procesos

Los capítulos 3 y 4 se centraron en el desarrollo de diagrama de control tipo Shewhart para la vigilancia de una característica de la calidad, ya que se basaron en los principios de los diagramas de control desarrollados por el Dr. Samuel Walter Shewhart. Una de las principales desventajas de cualquier carta de control de Shewhart es que sólo utiliza la información del proceso contenida en el último punto graficado e ignora cualquier información ofrecida por la secuencia completa de puntos. Esta característica hace que la carta de control de Shewhart sea relativamente insensible a los corrimientos pequeños del proceso (por ejemplo, del orden de aproximadamente \(1,\!5\sigma\) o menos). Desde luego, es posible aplicar otros criterios en las cartas de Shewhart, como las pruebas de corridas y el uso de límites de advertencia, cuya intención es incorporar información del conjunto completo de puntos en el procedimiento de decisión. Sin embargo, el uso de reglas de sensibilidad adicionales reduce la simplicidad y facilidad de interpretación de las cartas de control de Shewhart. Además, el uso de estas reglas de sensibilidad puede ocasionar una reducción drástica de la longitud media de corrida de la carta de control cuando el proceso está bajo control, lo cual puede ser indeseable. Debido a lo anteriormente expuesto en este núcleo temático se presentan dos alternativas, muy efectivas, a las cartas de control de Shewhart cuando son de interés los corrimientos pequeños: el diagrama de control de suma acumulada (CUSUM) y el diagrama de control del promedio móvil ponderado exponencialmente (EWMA).

Por otro lado, los diagramas de control tipo Shewhart expuestos en estos capítulos se ocupan del monitoreo y control de procesos en la perspectiva básica de una sola variable, es decir, se ha supuesto que hay una sola variable de salida o característica de la calidad del proceso de interés. En la práctica, sin embargo, muchos, de no ser que la mayoría, de los escenarios de monitoreo y control incluyen varias variables relacionadas. Aun cuando una posible solución es la aplicación de cartas de control de una sola variable a cada variable individual, este enfoque es ineficiente y puede llevar a conclusiones equivocadas. Se necesitan métodos multivariados que consideren a las variables en conjunto, este es el caso del diagrama multivariado \(T^2\) de Hotelling.

5.1 Diagramas de Control CUSUM para Monitorear la Media del proceso

Este diagrama de control fue propuesto por primera vez por Page (1954), y el nombre de CUSUM se debe a que en este diagrama de control el estadístico que se grafica es la suma acumulada de las desviaciones con respecto a la media global. El estadístico CUSUM se define de la siguiente manera, suponer que se colectan muestras de tamaño \(n \geq 1\), y que \(\bar{X}_j\) es el promedio de la \(j\)-ésima muestra. Entonces, si \(\mu_0\), es el objetivo para la media del proceso, el diagrama de control de suma acumulada se construye graficando la cantidad

CUSUM

\[ \begin{equation} C_i = \sum_{j = 1}^i \left(\bar{X}_j -\mu_0 \right) \end{equation} \tag{5.1}\]

contra la muestra \(i\). A \(C_i\) se le llama la suma acumulada hasta la \(i\)-ésima muestra, incluyéndola. Debido a que combina información de varias muestras, las cartas de suma acumulada son más efectivas con tamaño de muestra \(n = 1\). Esto hace que el diagrama de control de suma acumulada sea un buen candidato para usarse en la industria química y de procesamiento, donde es frecuente que los subgrupos sean de tamaño uno, y en fabricación de piezas discretas con medición automática de cada pieza y control en línea mediante el uso directo de una microcomputadora en el centro de trabajo.

Dado que \(C_i\) acumula las desviaciones de cada media muestral con respecto a la media es fácil interpretar que si el proceso se mantiene bajo control en el valor objetivo \(\mu_0\), las medias muestral \(\bar{X}_j\) fluctuarán de manera aleatoria alrededor del valor objetivo \(\mu_0\), por lo tanto la suma acumulada se alternará de manera aleatoria alrededor del cero. Sin embargo si la media experimenta un corrimiento ascendente a un valor \(\mu_1 > \mu_0\), por ejemplo, entonces se producirá una corrida ascendente o positiva en la suma acumulada \(C_i\). Recíprocamente, si la media experimenta un corrimiento descendente a un valor \(\mu_1 < \mu_0\), entonces \(C_i\) describirá una corrida descendente o negativa. Por lo tanto, si se desarrolla una tendencia en los puntos graficados, sea ascendente o descendente, ésta deberá considerarse como evidencia de que la media del proceso se ha corrido y deberá realizarse la búsqueda de alguna causa atribuible.

Fíjese que la ecuación 5.1 se puede reescribir como

CUSUM recursiva

\[ \begin{equation} C_i = \left(\bar{X}_i -\mu_0 \right) + \sum_{j = 1}^{i - 1}\left(\bar{X}_j -\mu_0 \right) = \left(\bar{X}_i -\mu_0 \right) + C_{i -1} \end{equation}. \tag{5.2}\]

En la ecuación anterior, el valor inicial de la CUSUM, \(C_0\), se toma como cero. Esta forma permite calcular \(C_i\) de manera recursiva.

Para mostrar una aplicación de esta teoría, considere los datos de la Tabla 5.1. Las primeras 20 observaciones se generaron de una variable aleatoria con distribución normal con media \(\mu = 100\) y desviación estándar \(\sigma = 5\), mientras que las últimas 10 se generaron de una variable aleatoria con distribución normal con media \(\mu = 105\) y desviación estándar \(\sigma = 5\). Por consiguiente, podría considerarse que estas 10 últimas observaciones se sacaron del proceso cuando estaba fuera de control (es decir, después que el proceso ha experimentado un corrimiento en la media de \(1\sigma\)). Para aplicar la CUSUM de la ecuación 5.2 a estas observaciones, se tomará \(x_i = \bar{x}_i\) (puesto que el tamaño de la muestra es \(n = 1\)), y se toma como valor objetivo \(\mu_0 = 100\). Por lo tanto la CUSUM es

\[ \begin{align*} C_i & = \left({x}_i -100 \right) + \sum_{j = 1}^{i - 1}\left({x}_j - 100 \right)\\ & =\left({x}_i -100 \right) + C_{i-1} \end{align*} \:. \]

La columna (b) de la Tabla 5.1 contiene las diferencias \(x_i-100\), y en la columna (c) se calculan las sumas acumuladas.

Código 
```{r}
#| label: tbl-datos-cusum
#| tbl-cap: "Datos para el ejemplo del diagrama CUSUM"
#| table-caption: left
#| class: output

set.seed(347)
data_cusum <- tibble::tibble(
  `muestra` = 1:30,
  `x_i` = c(
    sample(
      rnorm(n = 10000, mean = 100, sd = 5),
      size = 20, replace = TRUE
    ),
    sample(
      rnorm(n = 10000, mean = 105, sd = 5),
      size = 10, replace = TRUE
    )
  )
)
data_cusum <- data_cusum |>
  dplyr::mutate(
    `col_b` = x_i - 100,
    `col_c` = cumsum(col_b)
  ) |>
  data.table::data.table()

knitr::kable(
  data_cusum,
  # table.attr = 'data-quarto-disable-processing="true"',
  # booktabs = TRUE,
  format = "markdown",
  digits = 4,
  format.args = list(decimal.mark = ",", big.mark = "."),
  col.names = c(
    "$Muestra \\: (i)$", "$(a)\\:x_i$", "$(b)\\:x_i - 100$",
    "$(c)\\: \\left(x_i - 100 \\right) + C_{i -1}$"
  ),
  align = "cccc",
  escape = FALSE
)
```
Tabla 5.1: Datos para el ejemplo del diagrama CUSUM
\(Muestra \: (i)\) \((a)\:x_i\) \((b)\:x_i - 100\) \((c)\: \left(x_i - 100 \right) + C_{i -1}\)
1 96,9663 -3,0337 -3,0337
2 97,2494 -2,7506 -5,7843
3 103,6353 3,6353 -2,1490
4 96,7796 -3,2204 -5,3694
5 110,3248 10,3248 4,9554
6 97,7647 -2,2353 2,7202
7 97,4743 -2,5257 0,1945
8 100,8277 0,8277 1,0221
9 91,8081 -8,1919 -7,1698
10 107,2495 7,2495 0,0797
11 94,6188 -5,3812 -5,3015
12 105,6507 5,6507 0,3492
13 94,9607 -5,0393 -4,6901
14 107,0322 7,0322 2,3422
15 102,3922 2,3922 4,7343
16 100,0102 0,0102 4,7445
17 96,8418 -3,1582 1,5863
18 97,6166 -2,3834 -0,7971
19 97,9844 -2,0156 -2,8127
20 99,4052 -0,5948 -3,4075
21 110,7501 10,7501 7,3426
22 102,4706 2,4706 9,8132
23 108,5581 8,5581 18,3713
24 113,1497 13,1497 31,5209
25 112,3129 12,3129 43,8338
26 107,0480 7,0480 50,8818
27 110,5134 10,5134 61,3952
28 95,7660 -4,2340 57,1612
29 107,5611 7,5611 64,7223
30 98,1360 -1,8640 62,8584


La Figura 5.1 (a) muestra la gráfica de la CUSUM dada en la columna (c) de la Tabla 5.1. Obsérvese que para las primeras 20 observaciones, para las que \(\mu = 100\), la CUSUM tiende a presentar una desalineación lenta, en este caso manteniendo valores próximos a cero. Sin embargo, en las 10 últimas observaciones, para la que la media se ha corrido a \(\mu = 105\), se desarrolla una marcada tendencia ascendente, lo que representa una evidencia de que la media del proceso se ha corrido y deberá realizarse la búsqueda de alguna causa asignable.

La Figura 5.1 (b) presenta el diagrama tipo Shewhart estándar para las observaciones individuales mostradas en la columna (a) de la Tabla 5.1. Como se puede notar, ninguno de estos puntos se sale de los límites de control, por lo que no se tiene evidencia sólida de que el proceso está fuera de control. Obsérvese que hay un indicio de un corrimiento en el nivel medio del proceso en los 10 últimos puntos, ya que con excepción de dos puntos todos los demás se localizan arriba de la línea central. Sin embargo, si se confía en la señal tradicional de un proceso fuera de control (es decir, uno o más puntos que rebasen uno de los límites de control 3 sigma), entonces la carta de control de Shewhart ha fallado para detectar el corrimiento.

Código 
```{r}
#| label: fig-cusum-x
#| fig-cap: "Gráfica de la CUSUM y diagrama de control $x$"
#| fig-subcap:
#|   - "Gráfica de la CUSUM"
#|   - "Diagrama $x$"
#| layout-ncol: 2

# Gráfica de la cusum
ggplot(data = data_cusum, aes(x = muestra, y = col_c)) +
  geom_line(
    mapping = aes(x = muestra, y = 0), color = "red", lty = "dashed"
  ) +
  geom_point(
    mapping = aes(x = muestra, y = col_c), color = "black"
  ) +
  geom_line(
    mapping = aes(x = muestra, y = col_c), color = "black"
  ) +
  geom_vline(
    xintercept = (20 + 21) / 2, linetype = 2, color = "blue"
  ) +
  annotate(
    geom = "segment", x = 8, y = -15, xend = 1, yend = -15,
    arrow = arrow(length = unit(2, "mm"))
  ) +
  annotate(
    geom = "segment", x = 12, y = -15, xend = 20.5, yend = -15,
    arrow = arrow(length = unit(2, "mm"))
  ) +
  annotate(
    geom = "segment", x = 23, y = -15, xend = 20.5, yend = -15,
    arrow = arrow(length = unit(2, "mm"))
  ) +
  annotate(
    geom = "segment", x = 27, y = -15, xend = 30, yend = -15,
    arrow = arrow(length = unit(2, "mm"))
  ) +
  annotate(
    "text",
    x = 10, y = -15,
    label = latex2exp::TeX(
      "$\\mu = 100$", bold = TRUE, italic = TRUE
    )
  ) +
  annotate(
    "text",
    x = 25, y = -15,
    label = latex2exp::TeX(
      "$\\mu = 105$", bold = TRUE, italic = TRUE
    )
  ) +
  labs(
    x = latex2exp::TeX(
      "Número de muestra", bold = TRUE, italic = TRUE
    ),
    y = latex2exp::TeX(
      "$C_i$", bold = TRUE, italic = TRUE
    )
  ) +
  scale_x_continuous(
    breaks = seq(from = 0, to = 30, by = 2),
  ) +
  theme_bw()


# diagrama x
ggplot(
  data = data_cusum,
  aes(x = muestra, y = x_i)
) +
  geom_step(
    mapping = aes(x = muestra, y = 100 - 3 * 5),
    direction = "hv", color = "red", lty = "dashed"
  ) +
  geom_step(
    mapping = aes(x = muestra, y = 100),
    direction = "hv", color = "blue"
  ) +
  geom_step(
    mapping = aes(x = muestra, y = 100 + 3 * 5),
    direction = "hv", color = "red", lty = "dashed"
  ) +
  geom_point(
    mapping = aes(x = muestra, y = x_i), size = 2
  ) +
  geom_line(
    mapping = aes(x = muestra, y = x_i), color = "black"
  ) +
  pammtools::geom_stepribbon(
    aes(ymin = 100 - 3 * 5, ymax = 100),
    fill = "green", alpha = 0.2, direction = "hv"
  ) +
  pammtools::geom_stepribbon(
    aes(ymin = 100, ymax = 100 + 3 * 5),
    fill = "green", alpha = 0.2, direction = "hv"
  ) +
  geom_vline(xintercept = (20 + 21) / 2, linetype = 2, color = "blue") +
  annotate(
    "text",
    x = 30 + 0.9, y = 100 - 3 * 5,
    label = latex2exp::TeX("$LCI$"),
    colour = "red", parse = TRUE
  ) +
  annotate(
    "text",
    x = 30 + 0.75, y = 100,
    label = latex2exp::TeX("$LC$"),
    colour = "blue"
  ) +
  annotate(
    "text",
    x = 30 + 0.9, y = 100 + 3 * 5,
    label = latex2exp::TeX("$LCS$"),
    colour = "red"
  ) +
  annotate(
    geom = "segment", x = 8, y = 90, xend = 1, yend = 90,
    arrow = arrow(length = unit(2, "mm"))
  ) +
  annotate(
    geom = "segment", x = 12, y = 90, xend = 20.5, yend = 90,
    arrow = arrow(length = unit(2, "mm"))
  ) +
  annotate(
    geom = "segment", x = 23, y = 90, xend = 20.5, yend = 90,
    arrow = arrow(length = unit(2, "mm"))
  ) +
  annotate(
    geom = "segment", x = 27, y = 90, xend = 30, yend = 90,
    arrow = arrow(length = unit(2, "mm"))
  ) +
  annotate(
    "text",
    x = 10, y = 90,
    label = latex2exp::TeX(
      "$\\mu = 100$",
      bold = TRUE, italic = TRUE
    )
  ) +
  annotate(
    "text",
    x = 25, y = 90,
    label = latex2exp::TeX(
      "$\\mu = 105$",
      bold = TRUE, italic = TRUE
    )
  ) +
  labs(
    x = latex2exp::TeX(
      "Número de muestra",
      bold = TRUE, italic = TRUE
    ),
    y = latex2exp::TeX(
      "$x_i$",
      bold = TRUE, italic = TRUE
    )
  ) +
  scale_x_continuous(
    breaks = seq(from = 0, to = 30, by = 2),
    limits = c(0, 31)
  ) +
  theme_bw()
```
(a) Gráfica de la CUSUM
(b) Diagrama \(x\)
Figura 5.1: Gráfica de la CUSUM y diagrama de control \(x\)

Desde luego, la gráfica de la suma acumulada de la Figura 5.1 (a) no es una carta de control, ya que carece de límites de control estadísticos. Hay dos formas de representar gráficos de control CUSUM, la CUSUM tabular (o algorítmica) y en la forma de máscara V de la CUSUM. De estas dos representaciones, es preferible la CUSUM tabular.

5.1.1 La CUSUM Tabular o Algorítmica para Monitorear la Media del Proceso

Es posible construir CUSUM tanto para observaciones individuales como para los promedios de los subgrupos racionales. El caso de observaciones individuales ocurre con mucha frecuencia en la práctica, por lo que se tratará primero esta situación. Después se verá cómo modificar este resultado para subgrupos racionales.

Sea \(x_i\) la \(i\)-ésima observación del proceso. Supóngase, que cuando el proceso está bajo control, \(x_i\) tiene una distribución normal con media \(\mu_0\) y desviación estándar \(\sigma\). Se supone que \(\sigma\) es conocida o que se cuenta con una estimación de la misma. Más adelante se discutirá el monitoreo de \(\sigma\) con una CUSUM.

En ocasiones \(\mu_0\) se considera como el valor “objetivo” para la característica de la calidad \(X\).

La CUSUM tabular funciona acumulando las desviaciones de \(\mu_0\) que están arriba del objetivo con un estadístico \(C^+\) y acumulando las desviaciones de \(\mu_0\) que están abajo del objetivo con otro estadístico \(C^-\). A los estadísticos \(C^+\) y \(C^-\) se les denomina CUSUM unilateral superior e inferior, respectivamente. Y se calculan a través de las siguientes ecuaciones

CUSUM unilateral

\[ \begin{align*} C_i^+ & = \mathrm{máx}\left\{0, \, X_i - \left (\mu_0 + K \right ) + C_{i - 1}^+ \right\}\\ C_i^- & = \mathrm{máx}\left\{0, \, \left (\mu_0 - K \right ) - X_i + C_{i - 1}^ - \right\}\\ \end{align*}\,, \tag{5.3}\] donde los valores iniciales son \(C_0^+= C_0^- = 0\)

En la ecuación 5.3, a \(K\) suele llamársele el valor de referencia (o tolerancia, o valor laxo), y con frecuencia se establece como la mitad de la distancia entre el objetivo \(\mu_0\) y el valor fuera de control de la media \(\mu_1\) que el analista está interesado en detectar con rapidez.

Por tanto, si el cambio se expresa en unidades de desviación estándar como \(\mu_1 = \mu_0 + \delta\sigma\) (o \(\delta = \left|\mu_1 - \mu_0 \right|/\sigma\)), entonces \(K\) es la mitad del corrimiento o

Valor de referencia

\[ \begin{align*} K = \frac{\delta}{2} \ \sigma= \frac{\left|\mu_1 - \mu_0 \right|}{2} \end{align*}. \tag{5.4}\]

Se debe tener en cuenta que \(C_i^+\) y \(C_i^-\) acumulan desviaciones del valor objetivo \(\mu_0\) que son mayores que \(K\), reiniciando ambas cantidades en cero cuando se hacen negativas. Si \(C_i^+\) o \(C_i^-\) excede el intervalo de decisión \(H\), se considera que el proceso está fuera de control.

Se ha mencionado brevemente cómo elegir \(K\), pero ¿cómo seleccionar \(H\)? De hecho, la selección adecuada de estos dos parámetros es muy importante, ya que tiene un impacto sustancial sobre el desempeño de la CUSUM. Se abundará al respecto más adelante, pero un valor razonable de \(H\) es cinco veces la desviación estándar del proceso, es decir, \(H = 5\sigma\).

Es útil construir la representación gráfica de la CUSUM tabular. Esta representación gráfica suele llamarse diagramas de control CUSUM del status. Se construyen graficando \(C_i^+\) y \(-C_i^-\) contra el número de muestra, además, se grafíca el intervalo de decisión, acotado por las líneas horizontales \(H = 5\sigma\) y \(-H = -5\sigma\). Algunas variantes de este diagrama incluyen: barras verticales, en cuyo extremo superior se ubica el valor de \(C_i^+\) y en el extremo inferior se ubica el valor \(C_i^-\) en el periodo \(i\). También, suele agregarse la gráfica de las desviaciones \(x_i - \mu_0\) para cada periodo, aunque otras variantes incluyen los puntos muestrales. Con frecuencia esto ayuda al usuario de la carta de control a visualizar el desempeño real del proceso que ha llevado a un valor particular de la carta CUSUM.

Cuando se aplica el diagrama de control CUSUM, la acción que se toma luego de una señal fuera de control es similar a la de cualquier diagrama de control. Se debe buscar la causa asignable, tomar medidas correctivas si es necesario y luego reiniciar el CUSUM en cero. La CUSUM es particularmente útil para determinar cuándo ha ocurrido la causa asignable, tan solo se cuenta hacia atrás desde la señal fuera de control hasta el periodo de tiempo en la que la CUSUM superó el cero para encontrar el primer periodo del proceso posterior al corrimiento. Los Contadores \(N^+\) y \(N^-\) se usan para este fin.

En situaciones en que se requiere el ajuste de una variable manipulable para devolver el proceso en el objetivo \(\mu_0\), puede ser útil tener una estimación de la nueva media del proceso después del corrimiento. Ésta puede calcularse a partir de

Estimación de la media del proceso después del corrimiento

\[ \hat{\mu } = \begin{cases} \mu_{0} + K + \frac{C_{i}^{+}}{N^{+}}, & \text{ si }\: C_{i}^{+} > 5 \\ \mu_{0} - K - \frac{C_{i}^{-}}{N^{-}}, & \text{ si }\: C_{i}^{-} > 5 \end{cases} \tag{5.5}\]

Cabe hacer notar que las pruebas para corridas, y otras reglas de sensibilidad como las reglas de zona, no pueden aplicarse con seguridad a la CUSUM, ya que los valores sucesivos de \(C_i^+\) y \(C_i^-\) no son independientes. De hecho, la CUSUM puede considerarse como un promedio ponderado, en la que las ponderaciones son estocásticas o aleatorias.

Ejemplo 5.1 (Aplicación de la CUSUM tabular)
Diseñe la CUSUM tabular para los datos de la Tabla 5.1.

Para diseñar la CUSUM tabular es necesario recordar que el valor objetivo es \(\mu_0 = 100\), el tamaño del subgrupo es \(n = 1\), la desviación estándar del proceso es \(\sigma = 5\), y se supondrá que la magnitud del corrimiento que se desea detectar es \(\delta \sigma = 1(5)=5\). Por lo tanto, el valor fuera de control de la media del proceso es \(\mu_1 = \mu_0 + \delta \sigma = 100 + 5 = 105\). Se usará una CUSUM tabular con \(K=5/2\) (ya que el tamaño del corrimiento es \(1\sigma\) y \(\sigma = 5\)) y \(H = 25\) (ya que el valor recomendado del intervalo decisión es \(H = 5\sigma = 5(5) = 25\)).

En la Tabla 5.2 se presenta el esquema de CUSUM tabular. Las ecuaciones para \(C_i^+\) y \(C_i^-\) son

\[ \begin{align*} C_i^+ & = \mathrm{máx}\left\{0, \, x_i - \left (100 + \frac{5}{2} \right ) + C_{i - 1}^+ \right\} = \mathrm{máx}\left\{0, \, x_i - 102,5 + C_{i - 1}^+ \right\}\\ C_i^- & = \mathrm{máx}\left\{0, \, \left (\mu_0 - \frac{5}{2} \right ) - x_i + C_{i - 1}^ - \right\} = \mathrm{máx}\left\{0, \, 97,5 - x_i + C_{i - 1}^ - \right\} \end{align*}. \]

Para ilustrar los cálculos de la CUSUM unilateral, considérese los periodos 1 y 2.

  • Periodo 1:

\[ \begin{align*} C_1^+ & = \mathrm{máx}\left\{0, \, 96,9663 - 102,5 + 0 \right\} = \mathrm{máx}\left\{0, \, -5,5337 \right\} = 0\\ C_1^- & = \mathrm{máx}\left\{0, \, 97,5 - 96,9663 + 0 \right\}= \mathrm{máx}\left\{0, \, 0,5337 \right\} = 0,5337 \end{align*}. \]

  • Periodo 2:

\[ \begin{align*} C_2^+ & = \mathrm{máx}\left\{0, \, 97,2494 - 102,5 + 0 \right\} = \mathrm{máx}\left\{0, \, -5,2506 \right\} = 0\\ C_2^- & = \mathrm{máx}\left\{0, \, 97,5 - 97,2494 + 0,5337 \right\} = \mathrm{máx}\left\{0, \, 0,7843 \right\} = 0,7843 \end{align*}. \]

En los paneles (a) y (b) de la Tabla 5.2 se resumen los cálculos restantes. Las cantidades \(N^+\) y \(N^-\) de la Tabla 5.2 indican el número de periodos consecutivos en los que los \(C_i^+\) y \(C_i^-\) han sido diferentes de cero.

Código 
```{r}
#| label: tbl-cusum-tabular
#| tbl-cap:  "La CUSUM tabular para los datos de la @tbl-datos-cusum"
#| table-caption: left
#| class: output

cusum_stats <- function(x, mu_0 = 0, K, C_0 = 0) {
  c_mas <- numeric(length = length(x))
  n_mas <- numeric(length = length(x))
  c_menos <- numeric(length = length(x))
  n_menos <- numeric(length = length(x))
  c_mas[1] <- max(0, x[1] - (mu_0 + K) + C_0)
  n_mas[1] <- ifelse(c_mas[1] == 0, 0, 1)
  c_menos[1] <- max(0, (mu_0 - K) - x[1] + C_0)
  n_menos[1] <- ifelse(c_menos[1] == 0, 0, 1)
  for (i in 2:length(x)) {
    c_mas[i] <- max(0, x[i] - (mu_0 + K) + c_mas[i - 1])
    n_mas[i] <- ifelse(c_mas[i] != 0, 1 + n_mas[i - 1], 0)
    c_menos[i] <- max(0, (mu_0 - K) - x[i] + c_menos[i - 1])
    n_menos[i] <- ifelse(c_menos[i] != 0, 1 + n_menos[i - 1], 0)
  }
  return(list(c_mas = c_mas, n_mas = n_mas, c_menos = c_menos, n_menos = n_menos))
}

data_cusum_tabular <- data_cusum |>
  dplyr::select(
    muestra, x_i
  ) |>
  dplyr::mutate(
    col_b = x_i - (100 + abs(105 - 100) / 2),
    c_mas = cusum_stats(x_i, mu_0 = 100, K = abs(105 - 100) / 2)[[1]],
    n_mas = cusum_stats(x_i, mu_0 = 100, K = abs(105 - 100) / 2)[[2]],
    col_c = (100 - abs(105 - 100) / 2) - x_i,
    c_menos = cusum_stats(x_i, mu_0 = 100, K = abs(105 - 100) / 2)[[3]],
    n_menos = cusum_stats(x_i, mu_0 = 100, K = abs(105 - 100) / 2)[[4]]
  ) |>
  data.table::data.table()

knitr::kable(
  data_cusum_tabular,
  format = "markdown",
  digits = 4,
  format.args = list(decimal.mark = ",", big.mark = "."),
  col.names = c(
    "Periodo ($i$)", "$x_i$", "$x_i - 102,5$",
    "$C_i^+$", "$N^+$", "$97,5 - x_i$", "$C_i^-$",
    "$N^-$"
  ),
  align = "cccccccc",
  escape = FALSE
)
```
Tabla 5.2: La CUSUM tabular para los datos de la Tabla 5.1
Periodo (\(i\)) \(x_i\) \(x_i - 102,5\) \(C_i^+\) \(N^+\) \(97,5 - x_i\) \(C_i^-\) \(N^-\)
1 96,9663 -5,5337 0,0000 0 0,5337 0,5337 1
2 97,2494 -5,2506 0,0000 0 0,2506 0,7843 2
3 103,6353 1,1353 1,1353 1 -6,1353 0,0000 0
4 96,7796 -5,7204 0,0000 0 0,7204 0,7204 1
5 110,3248 7,8248 7,8248 1 -12,8248 0,0000 0
6 97,7647 -4,7353 3,0896 2 -0,2647 0,0000 0
7 97,4743 -5,0257 0,0000 0 0,0257 0,0257 1
8 100,8277 -1,6723 0,0000 0 -3,3277 0,0000 0
9 91,8081 -10,6919 0,0000 0 5,6919 5,6919 1
10 107,2495 4,7495 4,7495 1 -9,7495 0,0000 0
11 94,6188 -7,8812 0,0000 0 2,8812 2,8812 1
12 105,6507 3,1507 3,1507 1 -8,1507 0,0000 0
13 94,9607 -7,5393 0,0000 0 2,5393 2,5393 1
14 107,0322 4,5322 4,5322 1 -9,5322 0,0000 0
15 102,3922 -0,1078 4,4244 2 -4,8922 0,0000 0
16 100,0102 -2,4898 1,9345 3 -2,5102 0,0000 0
17 96,8418 -5,6582 0,0000 0 0,6582 0,6582 1
18 97,6166 -4,8834 0,0000 0 -0,1166 0,5416 2
19 97,9844 -4,5156 0,0000 0 -0,4844 0,0572 3
20 99,4052 -3,0948 0,0000 0 -1,9052 0,0000 0
21 110,7501 8,2501 8,2501 1 -13,2501 0,0000 0
22 102,4706 -0,0294 8,2207 2 -4,9706 0,0000 0
23 108,5581 6,0581 14,2788 3 -11,0581 0,0000 0
24 113,1497 10,6497 24,9284 4 -15,6497 0,0000 0
25 112,3129 9,8129 34,7413 5 -14,8129 0,0000 0
26 107,0480 4,5480 39,2893 6 -9,5480 0,0000 0
27 110,5134 8,0134 47,3027 7 -13,0134 0,0000 0
28 95,7660 -6,7340 40,5687 8 1,7340 1,7340 1
29 107,5611 5,0611 45,6298 9 -10,0611 0,0000 0
30 98,1360 -4,3640 41,2659 10 -0,6360 0,0000 0


Los cálculos de las CUSUM de la Tabla 5.2 indican que la CUSUM del lado superior en el periodo 25 es \(C_{25}^+ = 34,7413\). Puesto que se trata del primer periodo en el que \(C_i^+ > H = 25\), se concluirá que el proceso está fuera de control en ese punto. La CUSUM tabular también indica el momento probable en que ocurrió el corrimiento. El contador \(N^+\) registra el número de periodos consecutivos desde que la CUSUM del lado superior \(C_i^+\) subió por encima del valor cero. Puesto que \(N^ + = 5\) en el periodo 25, se concluirá que la última vez que el proceso estuvo bajo control fue en el periodo \(25 - 5 = 20\), por lo que el corrimiento ocurrió posiblemente en el periodo 20 y 21.

La gráfica de la CUSUM del status se muestra en la Figura 5.2, y en esta se hace evidente que el estado fuera de control se detecta en la muestra 25.

Código 
```{r}
#| label: fig-cusum-tabular
#| fig-cap: "Diagrama de control CUSUM del status para los datos de la @tbl-cusum-tabular"

color <- ifelse(
  data_cusum_tabular$c_mas > 25 | -data_cusum_tabular$c_menos < -25, "#FF0000", "#0000FF"
)

ggplot(
  data = data_cusum_tabular,
  aes(x = muestra, y = -c_menos, yend = c_mas)
) +
  geom_line(y = 25, color = "red", lty = "dashed") +
  geom_line(y = -25, color = "red", lty = "dashed") +
  geom_line(
    y = abs(105 - 100) / 2, color = "blue", lty = "dashed"
  ) +
  geom_line(
    y = -abs(105 - 100) / 2, color = "blue", lty = "dashed"
  ) +
  geom_errorbar(
    aes(ymin = -c_menos, ymax = c_mas),
    color = as.factor(color), linewidth = 1
  ) +
  geom_point(
    aes(x = muestra, y = x_i - 100),
    shape = "square", color = factor(color)
  ) +
  pammtools::geom_stepribbon(
    aes(ymin = -25, ymax = 25),
    fill = "green",
    alpha = 0.2, direction = "hv"
  ) +
  geom_vline(
    xintercept = (20 + 21) / 2, linetype = 2, color = "blue"
  ) +
  annotate(
    "text",
    x = 30 + 0.9, y = 25,
    label = latex2exp::TeX("$H$"),
    colour = "red",
    parse = TRUE
  ) +
  annotate(
    "text",
    x = 30 + 0.65, y = -25,
    label = latex2exp::TeX("$-H$"),
    colour = "red"
  ) +
  annotate(
    "text",
    x = 30 + 0.9, y = 2.5,
    label = latex2exp::TeX("$K$"),
    colour = "red",
    parse = TRUE
  ) +
  annotate(
    "text",
    x = 30 + 0.65, y = -2.5,
    label = latex2exp::TeX("$-K$"),
    colour = "red"
  ) +
  labs(
    x = latex2exp::TeX(
      "Periodo ($i$)",
      bold = TRUE, italic = TRUE
    ),
    y = latex2exp::TeX(
      "$CUSUM \\, Unilateral$",
      bold = TRUE, italic = TRUE
    )
  ) +
  scale_x_continuous(
    breaks = seq(from = 1, to = 30, by = 3),
    limits = c(0, 31)
  ) +
  scale_y_continuous(
    breaks = seq(from = -25, to = 50, by = 10),
  ) +
  theme_bw()
```
Figura 5.2: Diagrama de control CUSUM del status para los datos de la Tabla 5.2

Para ilustrar la aplicación de la ecuación 5.5 en este ejemplo, considere la CUSUM en el periodo 25, es decir, \(C_{25}^+ = 34,7413\), con \(N^+ = 5\). Por la ecuación 5.5, el nuevo promedio del proceso se estimará como

\[ \hat{\mu}=\mu_{0} - K - \frac{C_{}^{+}}{N^{+}} = 100 + 2,5 + \frac{34,7413}{5} = 109,4483\,. \]

La Figura 5.3 muestra el diagrama de control CUSUM del status con el paquete highcharter.

Código 
```{r}
#| label: fig-cusum-tabular-highcharter
#| fig-cap: "Diagrama de control CUSUM del statuscon el paquete `highcharter`"

highchart() |>
  hc_add_series(
    data = data_cusum_tabular,
    hcaes(x = muestra, low = 0, high = 25),
    marker = list(enabled = FALSE, visible = FALSE),
    type = "arearange", fillOpacity = 0.2, 
    name = "", color = "green",
    showInLegend = FALSE, 
    tooltip = list(pointFormat = "{NULL}")
  ) |>
  hc_add_series(
    data = data_cusum_tabular,
    hcaes(x = muestra, low = -25, high = 0),
    marker = list(enabled = FALSE, visible = FALSE),
    type = "arearange", fillOpacity = 0.2, 
    name = "", color = "blue",
    showInLegend = FALSE, 
    tooltip = list(pointFormat = "{NULL}")
  ) |>
  hc_add_series(
    data = data_cusum_tabular,
    type = "line", hcaes(x = muestra, y = 25), 
    useHTML = TRUE, color = "red", 
    dashStyle = "Dash", name = "<b><i>H</i></b>",
    step = "hv", marker = list(enabled = FALSE),
    tooltip = list(valueDecimals = 0)
  ) |>
  hc_add_series(
    data = data_cusum_tabular,
    type = "line", hcaes(x = muestra, y = -25), useHTML = TRUE,
    color = "red", dashStyle = "Dash", 
    name = "<b><i>-H</i></b>",
    step = "hv", marker = list(enabled = FALSE),
    tooltip = list(valueDecimals = 2)
  ) |>
  hc_add_series(
    data = data_cusum_tabular,
    type = "line", hcaes(x = muestra, y = 2.5), 
    useHTML = TRUE,
    color = "black", dashStyle = "Dash", 
    name = "<b><i>K</i></b>",
    step = "hv", marker = list(enabled = FALSE),
    tooltip = list(valueDecimals = 1)
  ) |>
  hc_add_series(
    data = data_cusum_tabular,
    type = "line", hcaes(x = muestra, y = -2.5), 
    useHTML = TRUE,
    color = "black", dashStyle = "Dash", 
    name = "<b><i>-K</i></b>",
    step = "hv", marker = list(enabled = FALSE),
    tooltip = list(valueDecimals = 1)
  ) |>
  hc_add_series(
    data = data_cusum_tabular, type = "errorbar", 
    linewhich = 4,
    hcaes(
      x = muestra, low = -c_menos, 
      high = c_mas, color = factor(color)
    ),
    showInLegend = TRUE, lineWidth = 3, useHTML = TRUE,
    name = "<b><i><strong>-C<sub>i</sub><sup>-</sup> - C<sub>i</sub><sup>+</sup></strong></i></b>",
    marker = list(enabled = TRUE),
    tooltip = list(valueDecimals = 2, useHTML = TRUE)
  ) |>
  hc_add_series(
    data = data_cusum_tabular,
    hcaes(
      x = muestra, y = x_i - 100, color = factor(color)
    ),
    showInLegend = TRUE, type = "scatter", useHTML = TRUE,
    name = "<b><i><p>x<sub>i</sub> - &mu;<sub>0</p></sub></i></b>",
    marker = list(
      enabled = TRUE, symbol = "circle",
      color = factor(color), radius = 3
    ),
    tooltip = list(
      valueDecimals = 2, useHTML = TRUE, crosshairs = TRUE,
      shared = TRUE, headerFormat = "<b><i>{point.x}</i></b><br>",
      pointFormat = "<b><i><p>x<sub>i</sub> - &mu;<sub>0</p></sub></i></b>:  {point.y}"
    )
  ) |>
  hc_xAxis(
    title = list(
      useHTML = TRUE,
      text = "<b><i>Periodo</i></b><br>"
    ),
    plotLines = list(
      list(
        color = "blue", value = 20.5, 
        width = 2, dashStyle = "Dash"
      )
    )
  ) |>
  hc_yAxis(
    title = list(
      useHTML = TRUE,
      text = "<b><i>CUSUM</i></b><br>"
    )
  ) |>
  hc_tooltip(
    crosshairs = TRUE, shared = TRUE, sort = TRUE
  ) |>
  hc_exporting(enabled = TRUE) |>
  hc_plotOptions(series = list(animation = FALSE))
```
Figura 5.3: Diagrama de control CUSUM del statuscon el paquete highcharter

La Figura 5.4 muestra otra versión del diagrama de control CUSUM del status con el paquete highcharter, en este casos las CUSUM unilateral se han representado con barras y se han agregados las observaciones \(x_i\).

Código 
```{r}
#| label: fig-cusum-tabular-highcharter2
#| fig-cap: "Diagrama de control CUSUM del status con el paquete `highcharter`"

highchart() |>
  hc_add_series(
    data = data_cusum_tabular,
    type = "column",
    hcaes(
      x = muestra, y = c_mas, color = color
    ),
    showInLegend = TRUE, lineWidth = 3, useHTML = TRUE,
    name = "<b><i>C<sub>i</sub><sup>+</sup></i></b>",
    marker = list(enabled = TRUE),
    tooltip = list(
      valueDecimals = 2, useHTML = TRUE, crosshairs = TRUE,
      shared = TRUE, headerFormat = "<b><i>{point.x}</i></b><br>",
      pointFormat = "<b><i>C<sub>i</sub><sup>+</sup></i></b>: {point.y}"
    ),
    min = min(data_cusum_tabular$c_mas),
    max = max(data_cusum_tabular$c_mas),
    plotOptions = list(stacking = "normal")
  ) |>
  hc_add_series(
    data = data_cusum_tabular,
    type = "column",
    hcaes(
      x = muestra, y = -c_menos
    ),
    color = "blue",
    showInLegend = TRUE, lineWidth = 3, useHTML = TRUE,
    name = "<b><i>C<sub>i</sub><sup>-</sup></i></b>",
    marker = list(enabled = TRUE),
    tooltip = list(
      valueDecimals = 2, useHTML = TRUE, crosshairs = TRUE,
      shared = TRUE, headerFormat = "<b><i>{point.x}</i></b><br>",
      pointFormat = "<b><i>C<sub>i</sub><sup>-</sup></i></b>: {point.y}"
    ),
    plotOptions = list(stacking = "normal")
  ) |>
  # hc_plotOptions(series = list(stacking = "normal")) |>
  hc_add_series(
    data = data_cusum_tabular,
    hcaes(x = muestra, low = 0, high = 25),
    marker = list(enabled = FALSE, visible = FALSE),
    type = "arearange", fillOpacity = 0.2, 
    name = "", color = "green",
    showInLegend = FALSE, tooltip = list(pointFormat = "{NULL}")
  ) |>
  hc_add_series(
    data = data_cusum_tabular,
    hcaes(x = muestra, low = -25, high = 0),
    marker = list(enabled = FALSE, visible = FALSE),
    type = "arearange", fillOpacity = 0.2, 
    name = "", color = "blue",
    showInLegend = FALSE, tooltip = list(pointFormat = "{NULL}")
  ) |>
  hc_add_series(
    data = data_cusum_tabular,
    type = "line", hcaes(x = muestra, y = 25), useHTML = TRUE,
    color = "red", dashStyle = "Dash", 
    name = "<b><i>H</i></b>",
    step = "hv", marker = list(enabled = FALSE),
    tooltip = list(valueDecimals = 0)
  ) |>
  hc_add_series(
    data = data_cusum_tabular,
    type = "line", hcaes(x = muestra, y = -25), useHTML = TRUE,
    color = "red", dashStyle = "Dash", 
    name = "<b><i>-H</i></b>",
    step = "hv", marker = list(enabled = FALSE),
    tooltip = list(valueDecimals = 2)
  ) |>
  hc_add_series(
    data = data_cusum_tabular,
    type = "line", hcaes(x = muestra, y = 2.5), useHTML = TRUE,
    color = "black", dashStyle = "Dash", 
    name = "<b><i>K</i></b>",
    step = "hv", marker = list(enabled = FALSE),
    tooltip = list(valueDecimals = 1)
  ) |>
  hc_add_series(
    data = data_cusum_tabular,
    type = "line", hcaes(x = muestra, y = -2.5), 
    useHTML = TRUE, color = "black", 
    dashStyle = "Dash", name = "<b><i>-K</i></b>",
    step = "hv", marker = list(enabled = FALSE),
    tooltip = list(valueDecimals = 1)
  ) |>
  hc_yAxis_multiples(
    list(
      title = list(
        text = "<p><strong>⟵</strong><em><strong>C<sub>i</sub><sup>-</sup> &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;C<sub>i</sub><sup>+</sup></strong></em><strong>⟶</strong>&nbsp; &nbsp; &nbsp;&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;</p>" # ,
        # align = "center", y = 0, margin = 0
      ),
      style = list(fontWeight = "bold", fontSize = "auto"),
      labels = list(
        format = "{value}",
        style = list(fontSize = "auto", color = "black")
      ),
      opposite = FALSE,
      # min = min(-data_cusum_tabular$c_menos),
      # max = max(data_cusum_tabular$c_mas),
      lineWidth = 2
    ),
    list(
      title = list(text = "<p><em><strong>x<sub>i</sub></strong></em></p>"),
      style = list(fontWeight = "bold", fontSize = "auto"),
      labels = list(
        format = "{value}"
      ),
      min = min(data_cusum_tabular$x_i),
      max = max(data_cusum_tabular$x_i),
      opposite = FALSE,
      lineWidth = 2
    )
  ) |>
  hc_add_series(
    data = data_cusum_tabular,
    hcaes(x = muestra, y = x_i),
    showInLegend = TRUE, type = "scatter", useHTML = TRUE,
    name = "<b><i>x<sub>i</sub></i></b>",
    marker = list(
      enabled = TRUE, radius = 3, lineColor = "black"
    ),
    tooltip = list(
      valueDecimals = 2, useHTML = TRUE, crosshairs = TRUE,
      shared = TRUE, headerFormat = "<b><i>{point.x}</i></b><br>",
      pointFormat = "<b><i>x<sub>i</sub></i></b>: {point.y}"
    ),
    yAxis = 1, xAxis = 0,
    plotOptions = list(stacking = "normal")
  ) |>
  hc_tooltip(split = TRUE) |>
  hc_xAxis(
    title = list(
      useHTML = TRUE,
      text = "<b><i>Periodo</i></b><br>"
    ),
    plotLines = list(
      list(
        color = "blue", value = 20.5, 
        width = 2, dashStyle = "Dash"
      )
    )
  ) |>
  hc_plotOptions(
    area = list(stacking = "normal"),
    column = list(stacking = "normal"),
    line = list(
      stacking = "normal",
      color = "blue",
      marker = list(
        stacking = "normal",
        fillColor = color,
        lineWidth = 2,
        lineColor = "blue"
      )
    )
  ) # |>
  # hc_add_dependency("plugins/multicolor_series.js")
```
Figura 5.4: Diagrama de control CUSUM del status con el paquete highcharter

La Figura 5.5 muestra el diagrama de control CUSUM del status con el paquete plotly.

Código 
```{r}
#| label: fig-cusum-tabular-plotly
#| fig-cap: "Diagrama de control CUSUM del status con el paquete `plotly`"

vline <- function(x = 0, color = "blue") {
  list(
    type = "line",
    y0 = 0,
    y1 = 1,
    yref = "paper",
    x0 = x,
    x1 = x,
    line = list(color = color, dash = "dash")
  )
}

data_cusum_tabular$sd <- abs(
  (data_cusum_tabular$c_mas + data_cusum_tabular$c_menos) / 2
)

data_cusum_tabular$y <- (data_cusum_tabular$c_mas - data_cusum_tabular$c_menos) / 2
color <- ifelse(
  data_cusum_tabular$c_mas > 25 | -data_cusum_tabular$c_menos < -25, "#FF0000", "#0000FF"
)

# library(plotly)
data_cusum_tabular |> 
  plot_ly(x = ~muestra) |> 
  add_trace(
    x = ~muestra,
    y = ~y,
    type = "scatter",
    mode = "none", # markers
    # marker = list(color = color),
    color = color,
    name = "<b><i>C<sub>i</sub><sup>+</sup> - C<sub>i</sub><sup>-</sup></i></b>",
    # name = plotly::TeX("C_{i}^{+} - C_{i}^{-}"),
    customdata = ~c_mas,
    text = ~-c_menos,
    hovertemplate = paste(
      "Periodo: %{x}<br>",
      "<b><i>C<sub>i</sub><sup>+</sup></i></b>",
      "%{customdata:.2f}<br>",
      "<b><i>C<sub>i</sub><sup>-</sup></i></b>",
      "%{text:.2f}",
      "<extra></extra>"
    ),
    error_y = list(
      array = ~sd,
      color = color
    )
  ) |> 
  add_trace(
    y = ~ x_i - 100,
    color = color,
    name = "<b><i>x<sub>i</sub> -  &mu;<sub>0</sub></i></b>",
    type = "scatter",
    mode = "markers",
    marker = list(color = color),
    # marker = list(color = factor(color)),
    hovertemplate = paste(
      "Periodo: %{x}<br>",
      "<b><i>x<sub>i</sub> -  &mu;<sub>0</sub></i></b> = %{y:.2f}<extra></extra>"
    )
  ) |> 
  add_lines(
    y = 0,
    name = "<b><i>LC</i></b>",
    line = list(
      shape = "linear",
      color = "blue",
      dash = "dash"
    ),
    hovertemplate = paste("<b><i>H = %{y:.0f}</i></b><br><extra></extra>")
  ) |> 
  add_lines(
    y = 25,
    name = "<b><i>H</i></b>",
    fill = "tozeroy",
    fillcolor = "green",
    showlegend = TRUE,
    hovertemplate = paste("<b><i>H = %{y:.0f}</i></b><br><extra></extra>"),
    line = list(
      shape = "linear",
      color = "red",
      dash = "dash"
    )
  ) |> 
  add_lines(
    y = -25,
    name = "<b><i>-H</i></b>",
    fill = "tozeroy",
    fillcolor = "blue",
    showlegend = TRUE,
    hovertemplate = paste("<b><i>-H = %{y:.0f}</i></b><br><extra></extra>"),
    line = list(
      shape = "linear",
      color = "red",
      dash = "dash"
    )
  ) |> 
  add_lines(
    y = 2.5,
    name = "<b><i>K</i></b>",
    line = list(
      shape = "linear",
      color = "black",
      dash = "dot"
    ),
    hovertemplate = paste("<b><i>K = %{y:.1f}</i></b><br><extra></extra>")
  ) |> 
  add_lines(
    y = -2.5,
    name = "<b><i>-K</i></b>",
    line = list(
      shape = "linear",
      color = "black",
      dash = "dot"
    ),
    hovertemplate = paste("<b><i>-K = %{y:.1f}</i></b><br><extra></extra>")
  ) |> 
  layout(
    legend = list(
      x = 0.5,
      y = -0.23,
      xref = "paper",
      yref = "paper",
      orientation = "h",
      xanchor = "center",
      showarrow = FALSE
    ),
    shapes = list(vline(20.5)),
    xaxis = list(
      title = "<b><i>Periodo (i)</i></b><br>",
      zeroline = FALSE
    ),
    yaxis = list(
      title = "<b><i>CUSUM</i></b><br>",
      zeroline = FALSE
    )
  ) |> 
  plotly::config(locale = "es", mathjax = "cdn") |> 
  fillOpacity(alpha = 0.3)
```
Figura 5.5: Diagrama de control CUSUM del status con el paquete plotly

5.1.2 Recomendaciones para el Diseño de la CUSUM Tabular

La CUSUM tabular se diseña eligiendo valores para el valor de referencia \(K\) y para el intervalo de decisión \(H\). Es común recomendar que estos parámetros se seleccionen para producir un buen desempeño de la longitud media de la corrida. Se han realizado diversos estudios analíticos del desempeño de la \(\textit{LMC}\) para la CUSUM. Con base en estos estudios, es posible ofrecer algunas recomendaciones generales para seleccionar \(H\) y \(K\).

Definir \(H = h\sigma\) y \(K = k\sigma\), donde \(\sigma\) es la desviación estándar de la variable muestral usada para formar la CUSUM. Utilizando \(h = 4\) o \(h = 5\) y \(k=1/2\) por lo general se obtendrá una CUSUM que tiene propiedades convenientes de la \(\textit{LMC}\) contra un corrimiento de aproximadamente \(1\sigma\) en la media del proceso.

Para ilustrar que tan bien funcionan las recomendaciones \(h = 4\) o \(h = 5\) con \(k = 1/2\), considere la longitud media de corrida bilaterales que se muestra en la Tabla 5.3, la cual fue tomada de D. C. Montgomery (2005, 415). Obsérvese que un corrimiento de \(1\sigma\) se detectará en 8,38 muestras (con \(k = 1/2\) y \(h = 4\)) o en 10,4 muestras (con \(k = 1/2\) y \(h = 5\)). En comparación, una carta de control de Shewhart para mediciones individuales requeriría 43,96 muestras, en promedio, para detectar este corrimiento.

Obsérvese así mismo por la Tabla 5.3 que \(h = 4\) resulta en una \({\textit{LMC}}_0 = 168\), mientras que \(h = 5\) resulta en una \({\textit{LMC}}_0 = 465\). Si se elige \(h = 4,77\), se obtendrá un CUSUM con , el cual coincide con el valor \({\textit{LMC}}_0\) para la diagrama de control de Shewhart con los límites \(3\sigma\) usuales.

Código 
```{r}
#| label: tbl-lmc-cusum
#| tbl-cap: "Desempeño de la $LMC$ de la CUSUM tabular con $k = 1/2$ y $h = 4$ o $h = 5$"
#| table-caption: left

lmc_cusum <- data.table::data.table(
  check.names = FALSE,
  corrimiento = c(0, 0.25, 0.50, 0.75, 1, 1.5, 2, 2.5, 3, 4),
  h_4 = c(168, 74.2, 26.6, 13.3, 8.38, 4.75, 3.34, 2.62, 2.19, 1.71),
  h_5 = c(465, 139, 38, 17, 10.4, 5.75, 4.01, 3.11, 2.57, 2.01)
)

knitr::kable(
  lmc_cusum,
  format = "markdown",
  table.attr = 'data-quarto-disable-processing="true"',
  col.names = c(
    "$\\text{Corrimiento de la media}\\, \\left(\\text{múltiplo de}\\,\\sigma \\right)$",
    "$h = 4$", "$h = 5$"
  ),
  align = "ccc",
  digits = 2,
  format.args = list(decimal.mark = ",", big.mark = "."),
  escape = FALSE
)
```
Tabla 5.3: Desempeño de la \(LMC\) de la CUSUM tabular con \(k = 1/2\) y \(h = 4\) o \(h = 5\)
\(\text{Corrimiento de la media}\, \left(\text{múltiplo de}\,\sigma \right)\) \(h = 4\) \(h = 5\)
0,00 168,00 465,00
0,25 74,20 139,00
0,50 26,60 38,00
0,75 13,30 17,00
1,00 8,38 10,40
1,50 4,75 5,75
2,00 3,34 4,01
2,50 2,62 3,11
3,00 2,19 2,57
4,00 1,71 2,01

En general, la elección de \(k\) querrá hacerse respecto del tamaño del corrimiento que quiera detectarse; es decir, \(k=1/2 \; \delta\), donde \(\delta\) es la magnitud del corrimiento en unidades de desviaciones estándar que quiere detectarse. Este enfoque es casi equivalente a minimizar el valor de \({\textit{LMC}}_1\) para detectar un corrimiento de magnitud \(\delta\) para una \(\textit{LMC}_0\) fija. Como ya se mencionó, un valor muy usado en la práctica es \(k=1/2\). Entonces, una vez que se ha seleccionado \(k\), \(h\) deberá elegirse para obtener el desempeño deseado de la \(\textit{LMC}_0\). Douglas M. Hawkins (1993a) presenta una tabla de valores de \(k\) y los valores correspondientes de \(h\) que producirán \({\textit{LMC}}_0=370\). En la Tabla 5.4 se presentan estos valores.

Código 
```{r}
#| label: tbl-lmcbc-cusum
#| tbl-cap: "Valores de $k$ y los valores correspondientes de $h$ que producen ${\\textit{LMC}}_0 = 370$ para la CUSUM tabular bilateral $[$de @hawkins1993cumulative$]$"
#| table-caption: left

lmc_cusum <- data.table::data.table(
  check.names = FALSE,
  `k` = seq(from = 0.25, to = 1.5, by = 0.25),
  `h` = c(8.01, 4.77, 3.34, 2.52, 1.99, 1.61)
)

knitr::kable(
  lmc_cusum,
  format = "markdown",
  col.names = c("$k$", "$h$"),
  digits = 2,
  align = "cc",
  format.args = list(decimal.mark = ",", big.mark = "."),
  escape = FALSE
)
```
Tabla 5.4: Valores de \(k\) y los valores correspondientes de \(h\) que producen \({\textit{LMC}}_0 = 370\) para la CUSUM tabular bilateral \([\)de Douglas M. Hawkins (1993a)\(]\)
\(k\) \(h\)
0,25 8,01
0,50 4,77
0,75 3,34
1,00 2,52
1,25 1,99
1,50 1,61

Siegmund (1985) propone unas fórmulas para aproximar la \({\textit{LMC}}\) para una CUSUM unilateral (es decir, \(C_i^+\) o \(C_i^-\)) para valores de \(h\) y \(k\) dados, la aproximación de Siegmund es

\(LMC\) para una CUSUM Unilateraal

\[ \begin{equation} LMC = \frac{e^{-2\Delta b} + 2\Delta b -1}{2\Delta ^{2}} \end{equation} \tag{5.6}\]

para \(\Delta \neq 0\), donde \(\Delta = \delta^* - k\) para la CUSUM unilateral superior \(C_i^+\), \(\Delta = -\delta^* -k\) para la \(\textit{CUSUM}\) unilateral inferior \(C_i^-\), \(b = h + 1,166\), y \(\delta^* = \left(\mu_1 - \mu_0 \right) / \sigma\). Si \(\Delta = 0\), puede usarse \(LMC=b^2\).

La cantidad \(\delta^*\) representa el corrimiento en la media, en unidades de \(\sigma\), para el que debe calcularse la \(LMC\). Por lo tanto, si \(\delta^* = 0\), \({\textit{LMC}}_0\) se calculará con la ecuación 5.6, mientras que si \(\delta^* \neq 0\), se calculará el valor de \({\textit{LMC}}_1\) correspondiente a un cambio de magnitud \(\delta^*\). Para obtener la \(LMC\) de la CUSUM bilateral a partir de las \(LMC\) de los dos estadísticos unilaterales (por ejemplo \({\textit{LMC}}^+\) y \({\textit{LMC}}^-\)) se usa

\(LMC\) de la CUSUM bilateral

\[ \begin{equation} \frac{1}{{\textit{LMC}}} = \frac{1}{{\textit{LMC}}^+} + \frac{1}{{\textit{LMC}}^-} \end{equation}. \tag{5.7}\]

Para ilustrar, considérese la CUSUM bilateral con \(k = 1/2\) y \(h = 5\). Para encontrar \({LMC}_0\), primero se calcularían los valores \({LMC}_0\) para los esquemas bilaterales (por ejemplo, \({\textit{LMC}}_0^+\) y \({\textit{LMC}}_0^-\)). Se hace \(\delta^* = 0\); entonces \(\Delta = 0 - k =0 - 1/2 = -1/2\), \(b = h + 1,166 = 5 + 1,166 = 6,166\), y por la ecuación 5.7

\[ LMC = \frac{e^{-2\left (\frac{1}{2} \right ) \left ( 6,166 \right )} + 2\left ( \frac{1}{2} \right ) \left ( 6,166 \right ) -1}{2\left ( \frac{1}{2} \right ) ^{2}} = 938,22. \]

Por simetría, se tiene que \({LMC}_0^+ = {LMC}_0^-\), de donde se obtiene por la ecuación 5.6, la \(\textit{LMC}\) bajo control de la CUSUM bilateral es

\[ \frac{1}{LMC} = \frac{1}{938,22} + \frac{1}{938,22} = \frac{100}{46.911} \Rightarrow LMC = 469,1. \]

Este valor es muy próximo al verdadero valor \({\textit{LMC}}_0\) de 465 que se presenta en la tabla Tabla 5.3. Si ocurre un corrimiento en la media de \(2\sigma\), entonces \(\delta^* = 2\), \(\delta=1,5\) para la CUSUM unilateral superior, \(\Delta = -2,5\) para la CUSUM unilateral inferior y, por las ecuaciones 5.6 y 5.7, la \({\textit{LMC}}_1\) aproximada de la CUSUM bilateral puede calcularse como ⋍ 3,89$. El valor exacto que se presenta en la Tabla 5.3 es 4,01.

5.1.3 La CUSUM Estandarizada

Muchos usuarios de la CUSUM prefieren estandarizar la variable \(X_i\) antes de realizar el cálculo de la CUSUM unilateral. Sea

\(X_i\) estandarizada

\[ \begin{equation} Z_i = \frac{X_i - \mu_0}{\sigma} \end{equation} \tag{5.8}\]

el valor estandarizado de \(X_i\). Entonces la CUSUM estandarizada se define como sigue

CUSUM unilateral estandarizada

\[ \begin{align*} C_i^+ & = \mathrm{máx}\left\{0, \, Z_i - k + C_{i - 1}^+ \right\}\\ C_i^- & = \mathrm{máx}\left\{0, \, - k - Z_i + C_{i - 1}^ - \right\}\\ \end{align*}\:, \tag{5.9}\] donde los valores iniciales son \(C_0^+= C_0^- = 0\)

La estandarización de la CUSUM ofrece dos ventajas:

  • La selección de los parámetros \(H\) y \(K\) en el diseño del diagrama de control CUSUM estandarizado no dependen de la escala (es decir, no depende de \(\sigma\)). Esto se debe a que como la media y la desviación estándar del estadístico dado en la ecuación 5.8 es 0 y 1, respectivamente. Por lo que, \(H = h\sigma_{Z_i} = h(1) = h\) y \(K = k \sigma = k(1) = k\). Generalmente, se toma \(h = 4\) o \(h = 5\) y \(k = 1/2\)

  • La CUSUM estandarizada lleva de manera natural a una CUSUM para controlar la variabilidad, como se verá en la Sección 5.1.7.

Ejemplo 5.2 (Aplicación de la CUSUM tabular estandarizada)
Diseñe la CUSUM tabular estandarizada para los datos de la Tabla 5.1.

Para que halla una correspondencia entre la CUSUM tabular diseñada en el Ejemplo 5.1 y la CUSUM tabular estandarizada que se quiere ejecutar en este ejemplo se debe establecer \(H = h = 5\) y \(K = k = 1/2\).

En la Tabla 5.2 se presenta el esquema de CUSUM tabular. Las ecuaciones para \(C_i^+\) y \(C_i^-\) son

\[ \begin{align*} C_i^+ & = \mathrm{máx}\left\{0, \, z_i - \frac{1}{2} + C_{i - 1}^+ \right\}\\ C_i^- & = \mathrm{máx}\left\{0, \, -\frac{1}{2} - z_i + C_{i - 1}^ - \right\} \end{align*}. \]

Para ilustrar los cálculos de la CUSUM tabular estandarizada, considérese los periodos 1 y 2.

  • Periodo 1:

\[ \begin{align*} C_1^+ & = \mathrm{máx}\left\{0, \, -0,6067 - 0,5 + 0 \right\} = \mathrm{máx}\left\{0, \, -1,1067 \right\} = 0\\ C_1^- & = \mathrm{máx}\left\{0, \, -0,5 - (-0,6067) + 0 \right\} = \mathrm{máx}\left\{0, \, 0,1067 \right\} = 0,1067 \end{align*}\:. \]

  • Periodo 2:

\[ \begin{align*} C_2^+ & = \mathrm{máx}\left\{0, \, -0,5501 - 0,5 + 0 \right\} = \mathrm{máx}\left\{0, \, -1,0501 \right\} = 0\\ C_2^- & = \mathrm{máx}\left\{0, \, -0,5 - (-0,5501) + 0,1067 \right\} = \mathrm{máx}\left\{0, \, 0,0501 \right\} = 0,0501 \end{align*}\:. \]

En los paneles (a) y (b) de la Tabla 5.5 se resumen los cálculos restantes. Las cantidades \(N^+\) y \(N^-\) de la Tabla 5.5 indican el número de periodos consecutivos en los que los \(C_i^+\) y \(C_i^-\) han sido diferentes de cero.

Código 
```{r}
#| label: tbl-cusum-tabular-estandarizada
#| tbl-cap:  "La CUSUM estandarizada tabular para los datos de la @tbl-datos-cusum"
#| table-caption: left
#| class: output

data_cusum_tabular_estandarizada <- data_cusum |>
  dplyr::transmute(
    muestra = muestra,
    z_i = (x_i - 100) / 5
  ) |>
  dplyr::mutate(
    col_b = z_i - 1 / 2,
    c_mas = cusum_stats(x = z_i, K = 1 / 2)[[1]],
    n_mas = cusum_stats(x = z_i, K = 1 / 2)[[2]],
    col_c = -1 / 2 - z_i,
    c_menos = cusum_stats(x = z_i, K = 1 / 2)[[3]],
    n_menos = cusum_stats(x = z_i, K = 1 / 2)[[4]]
  ) |>
  data.table::data.table()

knitr::kable(
  data_cusum_tabular_estandarizada,
  table.attr = 'data-quarto-disable-processing="true"',
  format = "markdown",
  digits = 4,
  format.args = list(decimal.mark = ",", big.mark = "."),
  col.names = c(
    "Periodo ($i$)", "$z_i$", "$z_i - \\frac{1}{2}$",
    "$C_i^+$", "$N^+$", "$- \\frac{1}{2} - z_i$",
    "$C_i^-$", "$N^-$"
  ),
  align = "cccccccc",
  escape = FALSE
)
```
Tabla 5.5: La CUSUM estandarizada tabular para los datos de la Tabla 5.1
Periodo (\(i\)) \(z_i\) \(z_i - \frac{1}{2}\) \(C_i^+\) \(N^+\) \(- \frac{1}{2} - z_i\) \(C_i^-\) \(N^-\)
1 -0,6067 -1,1067 0,0000 0 0,1067 0,1067 1
2 -0,5501 -1,0501 0,0000 0 0,0501 0,1569 2
3 0,7271 0,2271 0,2271 1 -1,2271 0,0000 0
4 -0,6441 -1,1441 0,0000 0 0,1441 0,1441 1
5 2,0650 1,5650 1,5650 1 -2,5650 0,0000 0
6 -0,4471 -0,9471 0,6179 2 -0,0529 0,0000 0
7 -0,5051 -1,0051 0,0000 0 0,0051 0,0051 1
8 0,1655 -0,3345 0,0000 0 -0,6655 0,0000 0
9 -1,6384 -2,1384 0,0000 0 1,1384 1,1384 1
10 1,4499 0,9499 0,9499 1 -1,9499 0,0000 0
11 -1,0762 -1,5762 0,0000 0 0,5762 0,5762 1
12 1,1301 0,6301 0,6301 1 -1,6301 0,0000 0
13 -1,0079 -1,5079 0,0000 0 0,5079 0,5079 1
14 1,4064 0,9064 0,9064 1 -1,9064 0,0000 0
15 0,4784 -0,0216 0,8849 2 -0,9784 0,0000 0
16 0,0020 -0,4980 0,3869 3 -0,5020 0,0000 0
17 -0,6316 -1,1316 0,0000 0 0,1316 0,1316 1
18 -0,4767 -0,9767 0,0000 0 -0,0233 0,1083 2
19 -0,4031 -0,9031 0,0000 0 -0,0969 0,0114 3
20 -0,1190 -0,6190 0,0000 0 -0,3810 0,0000 0
21 2,1500 1,6500 1,6500 1 -2,6500 0,0000 0
22 0,4941 -0,0059 1,6441 2 -0,9941 0,0000 0
23 1,7116 1,2116 2,8558 3 -2,2116 0,0000 0
24 2,6299 2,1299 4,9857 4 -3,1299 0,0000 0
25 2,4626 1,9626 6,9483 5 -2,9626 0,0000 0
26 1,4096 0,9096 7,8579 6 -1,9096 0,0000 0
27 2,1027 1,6027 9,4605 7 -2,6027 0,0000 0
28 -0,8468 -1,3468 8,1137 8 0,3468 0,3468 1
29 1,5122 1,0122 9,1260 9 -2,0122 0,0000 0
30 -0,3728 -0,8728 8,2532 10 -0,1272 0,0000 0


Los cálculos de las CUSUM estandarizadas de la Tabla 5.5 indican que la CUSUM estandarizada del lado superior en el periodo 25 es \(C_{25}^+ = 6,9483\). Puesto que se trata del primer periodo en el que \(C_i^+ > H = 5\), se concluirá que el proceso está fuera de control en ese punto. La CUSUM estandarizada tabular también indica el momento probable en que ocurrió el corrimiento. El contador \(N^+\) registra el número de periodos consecutivos desde que la CUSUM del lado superior \(C_i^+\) subió por encima del valor cero. Puesto que \(N^ + = 5\) en el periodo 25, se concluirá que la última vez que el proceso estuvo bajo control fue en el periodo \(25 - 5 = 20\), por lo que el corrimiento ocurrió posiblemente en el periodo 20 y 21.

La gráfica de la CUSUM estandarizada del status se muestra en la Figura 5.6, y en esta se hace evidente que el estado fuera de control se detecta en la muestra 25.

Código 
```{r}
#| label: fig-cusum-tabular-estandarizada
#| fig-cap: "Diagrama de control la CUSUM estandarizada del status  para los datos de la @tbl-cusum-tabular"

color <- ifelse(
  data_cusum_tabular_estandarizada$c_mas > 5 | -data_cusum_tabular_estandarizada$c_menos < -5,
  "#FF0000", "#0000FF"
)

ggplot(
  data = data_cusum_tabular_estandarizada,
  aes(x = muestra, y = -c_menos, yend = c_mas)
) +
  geom_line(y = 5, color = "red", lty = "dashed") +
  geom_line(y = -5, color = "red", lty = "dashed") +
  geom_line(y = 0.5, color = "blue", lty = "dashed") +
  geom_line(y = -0.5, color = "blue", lty = "dashed") +
  geom_errorbar(
    aes(ymin = -c_menos, ymax = c_mas),
    color = as.factor(color), linewidth = 1
  ) +
  geom_point(
    aes(x = muestra, y = z_i),
    shape = "square", color = factor(color)
  ) +
  pammtools::geom_stepribbon(
    aes(ymin = -5, ymax = 5),
    fill = "green",
    alpha = 0.2, direction = "hv"
  ) +
  geom_vline(
    xintercept = (20 + 21) / 2, linetype = 2, color = "blue"
  ) +
  annotate(
    "text",
    x = 30 + 0.9, y = 5,
    label = latex2exp::TeX("$H$"),
    colour = "red",
    parse = TRUE
  ) +
  annotate(
    "text",
    x = 30 + 0.65, y = -5,
    label = latex2exp::TeX("$-H$"),
    colour = "red"
  ) +
  annotate(
    "text",
    x = 30 + 0.9, y = 0.5,
    label = latex2exp::TeX("$K$"),
    colour = "red",
    parse = TRUE
  ) +
  annotate(
    "text",
    x = 30 + 0.65, y = -0.5,
    label = latex2exp::TeX("$-K$"),
    colour = "red"
  ) +
  labs(
    x = latex2exp::TeX(
      "Periodo ($i$)",
      bold = TRUE, italic = TRUE
    ),
    y = latex2exp::TeX(
      "$CUSUM \\, Unilateral$",
      bold = TRUE, italic = TRUE
    )
  ) +
  scale_x_continuous(
    breaks = seq(from = 1, to = 30, by = 3),
    limits = c(0, 31)
  ) +
  scale_y_continuous(
    breaks = seq(from = -5, to = 10, by = 2),
  ) +
  theme_bw()
```
Figura 5.6: Diagrama de control la CUSUM estandarizada del status para los datos de la Tabla 5.2

Para ilustrar la aplicación de la ecuación 5.5 en este ejemplo, considere la CUSUM en el periodo 25, es decir, \(C_{25}^+ = 34,7413\), con \(N^+ = 5\). Por la ecuación 5.5, el nuevo promedio del proceso se estimará como

\[ \hat{\mu}=\mu_{0} - K - \frac{C_{}^{+}}{N^{+}} = 100 + 2,5 + \frac{34,7413}{5} = 109,4483\,. \]

La siguiente gráfica muestra el diagrama de control CUSUM estandarizada del status con el paquete qcc.

Código 
```{r}
#| label: fig-cusum-tabular-estandarizada-qcc
#| fig-cap: "Diagrama de control CUSUM estandarizada del status  para los datos de la @tbl-cusum-tabular con el paquete `qcc`"

cusum <- cusum(
  data = data_cusum$x_i[1:20], sizes = 1, center = 100, 
  std.dev = 5, head.start = 0, decision.interval = 5, 
  se.shift = 1, newdata = data_cusum$x_i[21:30],
  newsizes = 1, plot = FALSE
)

plot(
  cusum, add.stats = FALSE, label.bounds = c("-H", "H"), 
  title = "", xlab = "Periodo", ylab = "CUSUM", 
  chart.all = TRUE
)
```
Figura 5.7: Diagrama de control CUSUM estandarizada del status para los datos de la Tabla 5.2 con el paquete qcc

La Figura 5.8 muestra el diagrama de control CUSUM estandarizada del status con el paquete highcharter.

Código 
```{r}
#| label: fig-cusum-tabular-estandarizada-highcharter
#| fig-cap: "Diagrama de control CUSUM estandarizado del status con el paquete `highcharter`"

highchart() |>
  hc_add_series(
    data = data_cusum_tabular_estandarizada,
    hcaes(x = muestra, low = 0, high = 5),
    marker = list(enabled = FALSE, visible = FALSE),
    type = "arearange", fillOpacity = 0.2, 
    name = "", color = "green",
    showInLegend = FALSE, tooltip = list(pointFormat = "{NULL}")
  ) |>
  hc_add_series(
    data = data_cusum_tabular_estandarizada,
    hcaes(x = muestra, low = -5, high = 0),
    marker = list(enabled = FALSE, visible = FALSE),
    type = "arearange", fillOpacity = 0.2, 
    name = "", color = "blue",
    showInLegend = FALSE, tooltip = list(pointFormat = "{NULL}")
  ) |>
  hc_add_series(
    data = data_cusum_tabular_estandarizada,
    type = "line", hcaes(x = muestra, y = 5), useHTML = TRUE,
    color = "red", dashStyle = "Dash", name = "<b><i>H</i></b>",
    step = "hv", marker = list(enabled = FALSE),
    tooltip = list(valueDecimals = 0)
  ) |>
  hc_add_series(
    data = data_cusum_tabular_estandarizada,
    type = "line", hcaes(x = muestra, y = -5), useHTML = TRUE,
    color = "red", dashStyle = "Dash", name = "<b><i>-H</i></b>",
    step = "hv", marker = list(enabled = FALSE),
    tooltip = list(valueDecimals = 0)
  ) |>
  hc_add_series(
    data = data_cusum_tabular_estandarizada,
    type = "line", hcaes(x = muestra, y = 0.5), useHTML = TRUE,
    color = "black", dashStyle = "Dash", name = "<b><i>K</i></b>",
    step = "hv", marker = list(enabled = FALSE),
    tooltip = list(valueDecimals = 1)
  ) |>
  hc_add_series(
    data = data_cusum_tabular_estandarizada,
    type = "line", hcaes(x = muestra, y = -0.5), useHTML = TRUE,
    color = "black", dashStyle = "Dash", name = "<b><i>-K</i></b>",
    step = "hv", marker = list(enabled = FALSE),
    tooltip = list(valueDecimals = 1)
  ) |>
  hc_add_series(
    data = data_cusum_tabular_estandarizada, 
    type = "errorbar", linewhich = 4,
    hcaes(
      x = muestra, low = -c_menos, high = c_mas, 
      color = factor(color)
    ),
    showInLegend = TRUE, lineWidth = 3, useHTML = TRUE,
    name = "<p><em><strong>-C<sub>i</sub><sup>-</sup>&nbsp;- C<sub>i</sub><sup>+</sup></strong></em></p>",
    marker = list(enabled = TRUE),
    tooltip = list(valueDecimals = 2, useHTML = TRUE)
  ) |>
  hc_add_series(
    data = data_cusum_tabular_estandarizada,
    hcaes(x = muestra, y = z_i, color = factor(color)),
    showInLegend = TRUE, type = "scatter", useHTML = TRUE,
    name = "<b><i><p>x<sub>i</sub> - &mu;<sub>0</p></sub></i></b>",
    marker = list(
      enabled = TRUE, symbol = "circle",
      color = factor(color), radius = 3
    ),
    tooltip = list(
      valueDecimals = 2, useHTML = TRUE, crosshairs = TRUE,
      shared = TRUE, headerFormat = "<b><i>{point.x}</i></b><br>",
      pointFormat = "<b><i><p>x<sub>i</sub> - &mu;<sub>0</p></sub></i></b>:  {point.y}"
    )
  ) |>
  hc_xAxis(
    title = list(
      useHTML = TRUE,
      text = "<b><i>Periodo</i></b><br>"
    ),
    plotLines = list(
      list(
        color = "blue", value = 21.5, width = 2, dashStyle = "Dash"
      )
    )
  ) |>
  hc_yAxis(
    title = list(
      useHTML = TRUE,
      text = "<b><i>CUSUM</i></b><br>"
    )
  ) |>
  hc_tooltip(
    crosshairs = TRUE, shared = TRUE, sort = TRUE
  ) |>
  hc_exporting(enabled = TRUE) |>
  hc_plotOptions(series = list(animation = FALSE))
```
Figura 5.8: Diagrama de control CUSUM estandarizado del status con el paquete highcharter

La Figura 5.9 muestra el diagrama de control CUSUM estandarizada del status con el paquete plotly.

Código 
```{r}
#| label: fig-cusum-tabular-estandarizada-plotly
#| fig-cap: "Diagrama de control CUSUM estandarizada del status con el paquete `plotly`"

vline <- function(x = 0, color = "blue") {
  list(
    type = "line",
    y0 = 0,
    y1 = 1,
    yref = "paper",
    x0 = x,
    x1 = x,
    line = list(color = color, dash = "dash")
  )
}

data_cusum_tabular_estandarizada$sd <- abs(
  (data_cusum_tabular_estandarizada$c_mas + data_cusum_tabular_estandarizada$c_menos) / 2
)

data_cusum_tabular_estandarizada$y <- (data_cusum_tabular_estandarizada$c_mas - data_cusum_tabular_estandarizada$c_menos) / 2

color <- ifelse(
  data_cusum_tabular_estandarizada$c_mas > 5 | -data_cusum_tabular_estandarizada$c_menos < -5,
  "#FF0000", "#0000FF"
)

# library(plotly)
data_cusum_tabular_estandarizada |>
  plot_ly(x = ~muestra) |>
  add_trace(
    x = ~muestra, y = ~y, type = "scatter", 
    mode = "none", # markers
    # marker = list(color = color),
    color = color,
    name = "<b><i>C<sub>i</sub><sup>+</sup> - C<sub>i</sub><sup>-</sup></i></b>",
    # name = plotly::TeX("C_{i}^{+} - C_{i}^{-}"),
    customdata = ~c_mas,
    text = ~-c_menos,
    hovertemplate = paste(
      "Periodo: %{x}<br>",
      "<b><i>C<sub>i</sub><sup>+</sup></i></b>",
      "%{customdata:.2f}<br>",
      "<b><i>C<sub>i</sub><sup>-</sup></i></b>",
      "%{text:.2f}",
      "<extra></extra>"
    ),
    error_y = list(
      array = ~sd,
      color = color
    )
  ) |>
  add_trace(
    y = ~z_i, 
    color = color,
    name = "<b><i>z<sub>i</sub></i></b>",
    # name = plotly::TeX("z_{i}"), 
    type = "scatter",
    mode = "markers",
    marker = list(color = color),
    hovertemplate = paste(
      "<b><i>Periodo = %{x}</i></b>",
      "<br><b><i>z<sub>i</sub></i></b> = %{y:.2f}<br><extra></extra>"
    )
  ) |>
  add_lines(
    y = 0, name = "<b><i>LC</i></b>", 
    line = list(
      shape = "linear", color = "blue", dash = "dash"
    ),
    hovertemplate = paste("<b><i>H = %{y:.0f}</i></b><br><extra></extra>")
  ) |>
  add_lines(
    y = 5, name = "<b><i>H</i></b>", fill = "tozeroy",
    fillcolor = "green", showlegend = TRUE,
    hovertemplate = paste("<b><i>H = %{y:.0f}</i></b><br><extra></extra>"),
    line = list(
      shape = "linear", color = "red", dash = "dash"
    )
  ) |>
  add_lines(
    y = -5, name = "<b><i>-H</i></b>", fill = "tozeroy",
    fillcolor = "blue", showlegend = TRUE,
    hovertemplate = paste(
      "<b><i>-H = %{y:.0f}</i></b><br><extra></extra>"
    ),
    line = list(
      shape = "linear", color = "red", dash = "dash"
    )
  ) |>
  add_lines(
    y = 0.5, name = "<b><i>K</i></b>",
    line = list(
      shape = "linear", color = "black", dash = "dot"
    ),
    hovertemplate = paste("<b><i>K = %{y:.1f}</i></b><br><extra></extra>")
  ) |>
  add_lines(
    y = -0.5, name = "<b><i>-K</i></b>",
    line = list(
      shape = "linear", color = "black", dash = "dot"
    ),
    hovertemplate = paste("<b><i>-K = %{y:.1f}</i></b><br><extra></extra>")
  ) |>
  layout(
    legend = list(
      x = 0.5, y = -0.23, xref = "paper", yref = "paper",
      orientation = "h", xanchor = "center", showarrow = FALSE
    ),
    shapes = list(vline(20.5)),
    xaxis = list(
      title = "<b><i>Periodo (i)</i></b><br>",
      zeroline = FALSE
    ),
    yaxis = list(
      title = "<b><i>CUSUM</i></b><br>",
      zeroline = FALSE
    )
  ) |>
  plotly::config(locale = "es", mathjax = "cdn") |>
  fillOpacity(alpha = 0.3)
```
Figura 5.9: Diagrama de control CUSUM estandarizada del status con el paquete plotly

5.1.4 Subgrupos Racionales

Aun cuando se ha presentado el desarrollo de la CUSUM tabular para el caso de observaciones individuales (\(n = 1\)), este resultado puede extenderse con facilidad al promedio de subgrupos racionales (\(n > 1\)). Simplemente se sustituye \(X_i\) con \(\bar{X}_i\) (el promedio muestral o de los subgrupos) en las formulas anteriores, y \(σ\) se sustituye con \({\sigma}_{\bar{X}_i} = \sigma/\sqrt{n}\).

En las cartas de Shewhart, el uso de promedios de subgrupos racionales mejora sustancialmente el desempeño de la carta de control. Sin embargo, no siempre ocurre lo mismo con la CUSUM. Si, por ejemplo, se tiene la alternativa de tomar una muestra de tamaño \(n = 1\) cada media hora y una muestra de subgrupos racionales de tamaño \(n = 5\) cada 2,5 horas (obsérvese que en cada caso se tiene la misma intensidad de muestreo), la CUSUM con frecuencia funcionará mejor con la elección de \(n = 1\) cada media hora. Sólo si existe alguna razón significativa de índole económica o alguna otra razón válida para tomar muestras de tamaño mayor que uno deberán usarse subgrupos de tamaño mayor que uno con la CUSUM.

5.1.5 Mejoramiento de la Respuesta de la CUSUM para Corrimientos Grandes

Se ha señalado que el diagrama de control CUSUM es muy efectivo para determinar corrimientos pequeños. Sin embargo, el diagrama de control CUSUM no es tan efectivo como el diagrama de Shewhart para detectar corrimientos grandes. Un enfoque para mejorar la habilidad del diagrama CUSUM para detectar los corrimientos grandes en el proceso es usar un procedimiento combinado CUSUM-Shewhart para el control en línea. Agregar la carta de control de Shewhart es una modificación muy simple del procedimiento de control de suma acumulada. Los límites de control de Shewhart deberán localizarse aproximadamente a 3,5 desviaciones estándar de la línea central o del valor objetivo \(\mu_0\). Una señal de fuera de control en cualquiera de las dos gráficas de control (o en ambas) constituye una señal de acción. Lucas (1982) presenta un buen análisis de esta técnica. En la columna (a) de la Tabla 5.6 se presentan las LMC de la CUSUM básica con \(k = 1/2\) y \(h = 5\). En la columna (b) de la Tabla 5.6 se presentan las \(\textit{LMC}\) de la CUSUM con límites de control de Shewhart agregados a las mediciones individuales. Como se sugirió antes, los límites de Shewhart están a \(3,5\sigma\). Al examinar estos valores de la \(\textit{LMC}\) se observa que la adición de los límites de control de Shewhart ha mejorado la habilidad del procedimiento para detectar corrimientos grandes y ha reducido apenas ligeramente la \({\textit{LMC}}_0\). De lo anterior, se concluye que un procedimiento combinado CUSUM-Shewhart es una forma efectiva de mejorar la respuesta de la CUSUM a los corrimientos grandes.

5.1.6 La Respuesta Inicial Rápida o Característica de Ventaja Anticipada

Este procedimiento fue propuesto por Lucas y Crosier (1982) para mejorar la sensibilidad de una CUSUM cuando arranca el proceso. Una mayor sensibilidad al inicio del proceso sería deseable si la acción correctiva no restableciera la media al valor objetivo. En esencia, la respuesta inicial rápida (FIR, por sus siglas en inglés) o ventaja anticipada consiste esencialmente en establecer los valores iniciales \(C_0^+\) en algún valor distinto de cero, generalmente en H/2. Se llama ventaja anticipada del 50%.

Código 
```{r}
#| label: tbl-cusum-LMC
#| tbl-cap: "Valores de $\\textit{LMC}$ para algunas modificaciones de la **CUSUM** básica con $k = 1/2$ y $h = 5$ (si se usan subgrupos de tamaño $n > 1$, entonces $\\sigma = {\\sigma}_{\\bar{X}} = \\sigma/\\sqrt{n}$)"
#| table-caption: left

cusum_ARL <- data.frame(
  corrimiento = c(0, 0.25, 0.5, 0.75, 1, 1.5, 2, 2.5, 3, 4),
  a = c(465, 139, 38, 17, 10.4, 5.75, 4.01, 3.11, 2.57, 2.01),
  b = c(
    391, 130.9, 37.2, 16.8, 10.2, 5.58, 3.77, 2.77, 2.1, 1.34
  ),
  c = c(430, 122, 28.7, 11.2, 6.35, 3.37, 2.36, 1.86, 1.54, 1.16),
  d = c(
    360, 113.9, 28.1, 11.2, 6.32, 3.37, 2.36, 1.86, 1.54, 1.16
  )
) |>
  data.table::data.table()

knitr::kable(
  cusum_ARL,
  format = "markdown",
  table.attr = 'data-quarto-disable-processing="true"',
  escape = FALSE,
  digits = 2,
  format.args = list(decimal.mark = ",", big.mark = "."),
  col.names = kableExtra::linebreak(
    c(
      "Corrimiento en la media (múltiplo de $\\sigma$)",
      "(a) CUSUM básica",
      "(b) CUSUM-Shewhart (límites de Shewhart en 3,5$\\sigma$)",
      "(c) CUSUM con FIR",
      "(d) CUSUM-Shewhart FIRT (límites de Shewhart en 3,5$\\sigma$)"
    )
  ),
  align = "ccccc"
)
```
Tabla 5.6: Valores de \(\textit{LMC}\) para algunas modificaciones de la CUSUM básica con \(k = 1/2\) y \(h = 5\) (si se usan subgrupos de tamaño \(n > 1\), entonces \(\sigma = {\sigma}_{\bar{X}} = \sigma/\sqrt{n}\))
Corrimiento en la media (múltiplo de \(\sigma\)) (a) CUSUM básica (b) CUSUM-Shewhart (límites de Shewhart en 3,5\(\sigma\)) (c) CUSUM con FIR (d) CUSUM-Shewhart FIRT (límites de Shewhart en 3,5\(\sigma\))
0,00 465,00 391,00 430,00 360,00
0,25 139,00 130,90 122,00 113,90
0,50 38,00 37,20 28,70 28,10
0,75 17,00 16,80 11,20 11,20
1,00 10,40 10,20 6,35 6,32
1,50 5,75 5,58 3,37 3,37
2,00 4,01 3,77 2,36 2,36
2,50 3,11 2,77 1,86 1,86
3,00 2,57 2,10 1,54 1,54
4,00 2,01 1,34 1,16 1,16

Ejemplo 5.3 (Aplicación de la CUSUM FIR)
Para ilustrar el procedimiento de ventaja anticipada, considere los datos de la tabla Tabla 5.7. Estos datos tienen un valor objetivo de 10 y una desviación estándar de 1, \(K = 1/2\) y \(H = 5\). Se usará un valor de ventaja anticipada del 50% de \(C_0^+ = C_0^- = H⁄2 = 5⁄2 = 2,\!5\). Las 10 primeras observaciones están bajo control con media igual al valor objetivo de 10. Dado que \(x_1 = 9,\!39\), las CUSUM del primer periodo serán:

\[ \begin{align*} C_1^+ & = \mathrm{máx}\left\{0, \, 9,\!39 - 10,\!5 + 2,\!5 \right\} = \mathrm{máx}\left\{0, \, 1,\!39 \right\} = 1,\!39\\ C_1^- & = \mathrm{máx}\left\{0, \, 9,\!5 - 9,\!39 + 2,\!5 \right\}= \mathrm{máx}\left\{0, \, 2,61 \right\} = 2,\!61 \end{align*}. \]

Código 
```{r}
#| label: tbl-cusum_FIR_0
#| tbl-cap: "CUSUM con una ventaja anticipada, $\\mu_0 = 10, K = 1/2$"
#| table-caption: left

set.seed(347)
data_cusum_FIR_0 <- tibble::tibble(
  `muestra` = 1:10,
  `x_i` = c(
    sample(
      rnorm(n = 10000, mean = 10, sd = 1),
      size = 10, replace = TRUE
    )
  )
)
data_cusum_FIR_0 <- data_cusum_FIR_0 |>
  dplyr::mutate(
    col_b = x_i - (10 + 1 / 2),
    c_mas = cusum_stats(x = x_i, mu_0 = 10, K = 1 / 2, C_0 = 2.5)[[1]],
    n_mas = cusum_stats(x = x_i, mu_0 = 10, K = 1 / 2, C_0 = 2.5)[[2]],
    col_c = (10 - 1 / 2) - x_i,
    c_menos = cusum_stats(x = x_i, mu_0 = 10, K = 1 / 2, C_0 = 2.5)[[3]],
    n_menos = cusum_stats(x = x_i, mu_0 = 10, K = 1 / 2, C_0 = 2.5)[[4]]
  ) |>
  data.table::data.table()

knitr::kable(
  data_cusum_FIR_0,
  # booktabs = TRUE,
  format = "markdown",
  digits = 2,
  format.args = list(decimal.mark = ",", big.mark = "."),
  col.names = c(
    "Periodo ($i$)", "$x_i$", "$x_i - 10,5$",
    "$C_i^+$", "$N^+$", "$9,5 - x_i$", "$C_i^-$",
    "$N^-$"
  ),
  align = "cccccccc",
  escape = FALSE
)
```
Tabla 5.7: CUSUM con una ventaja anticipada, \(\mu_0 = 10, K = 1/2\)
Periodo (\(i\)) \(x_i\) \(x_i - 10,5\) \(C_i^+\) \(N^+\) \(9,5 - x_i\) \(C_i^-\) \(N^-\)
1 9,39 -1,11 1,39 1 0,11 2,61 1
2 9,45 -1,05 0,34 2 0,05 2,66 2
3 10,73 0,23 0,57 3 -1,23 1,43 3
4 9,36 -1,14 0,00 0 0,14 1,57 4
5 12,06 1,56 1,56 1 -2,56 0,00 0
6 9,55 -0,95 0,62 2 -0,05 0,00 0
7 9,49 -1,01 0,00 0 0,01 0,01 1
8 10,17 -0,33 0,00 0 -0,67 0,00 0
9 8,36 -2,14 0,00 0 1,14 1,14 1
10 11,45 0,95 0,95 1 -1,95 0,00 0

Como se puede observar en los paneles (a) y (b) de la Tabla 5.7 ambos CUSUM disminuyen rápidamente desde el valor inicial. De hecho, a partir del período 4 en adelante \(C_i^+\) no se ve afectado por la ventaja anticipada, y a partir del período 5, \(C_i^-\) tampoco se ve afectado por la ventaja anticipada. Esto ocurre porque el proceso está bajo control en el valor objetivo de 10 y se observan varias muestras consecutivas cerca del valor objetivo.

Ahora supongamos que el proceso estuviera fuera de control al arranque del proceso, con una media de 11 y una desviación estándar de 1. La Tabla 5.8 presenta los datos que habría producido este proceso y las CUSUM resultantes. Observa que en la quinta muestra \(C_5^+\) excede el límite \(H = 5\). Si no se hubiera utilizado ninguna ventaja anticipada, es decir, \(C_0^+ = C_0^- = 0\), la CUSUM no habría superado \(H\) hasta la muestra número 8, ver el diagrama CUSUM básico (Figura 5.11).

En este ejemplo se han implementado el diagrama CUSUM FIR y el diagrama CUSUM básico (sin FIR) para detectar un cambio en el valor objetivo (\(\mu_0 = 10\)) hacia \(\mu_1 = 11\), lo que representa un cambio de \(1\sigma\). En el caso del diagrama CUSUM FIR el estado fuera de control es detectado en la muestra 5, una vez que éste ha ocurrido; este valor no difiere mucho de la \(\textit{LMC}\) indicada en la tabla Tabla 5.6 para este caso (\(\textit{LMC} = 6,\!35\)). Por otro lado, el diagrama CUSUM básico logra detectar el estado fuera de control en la muestra 8, después de haber ocurrido; la Tabla 5.6 indica una \(\textit{LMC} = 10,\!40\).

Código 
```{r}
#| label: tbl-cusum_FIR_1
#| tbl-cap: "CUSUM con una ventaja anticipada, $\\mu_0 = 11, K = 1/2$"
#| table-caption: left

set.seed(347)
data_cusum_FIR_1 <- tibble::tibble(
  `muestra` = 1:10,
  `x_i` = c(
    sample(
      rnorm(n = 10000, mean = 11, sd = 1),
      size = 10, replace = TRUE
    )
  )
)
data_cusum_FIR_1 <- data_cusum_FIR_1 |>
  dplyr::mutate(
    col_b = x_i - (10 + 1 / 2),
    c_mas = cusum_stats(x = x_i, mu_0 = 10, K = 1 / 2, C_0 = 2.5)[[1]],
    n_mas = cusum_stats(x = x_i, mu_0 = 10, K = 1 / 2, C_0 = 2.5)[[2]],
    col_c = (10 - 1 / 2) - x_i,
    c_menos = cusum_stats(x = x_i, mu_0 = 10, K = 1 / 2, C_0 = 2.5)[[3]],
    n_menos = cusum_stats(x = x_i, mu_0 = 10, K = 1 / 2, C_0 = 2.5)[[4]]
  ) |>
  data.table::data.table()

knitr::kable(
  data_cusum_FIR_1,
  format = "markdown",
  digits = 2,
  format.args = list(decimal.mark = ",", big.mark = "."),
  col.names = c(
    "Periodo ($i$)", "$x_i$", "$x_i - 10,5$",
    "$C_i^+$", "$N^+$", "$9,5 - x_i$", "$C_i^-$",
    "$N^-$"
  ),
  align = "cccccccc",
  escape = FALSE
)
```
Tabla 5.8: CUSUM con una ventaja anticipada, \(\mu_0 = 11, K = 1/2\)
Periodo (\(i\)) \(x_i\) \(x_i - 10,5\) \(C_i^+\) \(N^+\) \(9,5 - x_i\) \(C_i^-\) \(N^-\)
1 10,39 -0,11 2,39 1 -0,89 1,61 1
2 10,45 -0,05 2,34 2 -0,95 0,66 2
3 11,73 1,23 3,57 3 -2,23 0,00 0
4 10,36 -0,14 3,43 4 -0,86 0,00 0
5 13,06 2,56 5,99 5 -3,56 0,00 0
6 10,55 0,05 6,04 6 -1,05 0,00 0
7 10,49 -0,01 6,04 7 -0,99 0,00 0
8 11,17 0,67 6,70 8 -1,67 0,00 0
9 9,36 -1,14 5,57 9 0,14 0,14 1
10 12,45 1,95 7,52 10 -2,95 0,00 0

Este ejemplo demuestra los beneficios de la ventaja anticipada. Si el proceso comienza bajo control en el valor objetivo, los CUSUM caerán rápidamente a cero y la ventaja anticipada tendrá poco efecto en el rendimiento del procedimiento CUSUM. La Figura 5.10 ilustra esta propiedad de la ventaja utilizando los datos de las tablas 5.7 y 5.8. El gráfico CUSUM se produjo utilizando el paquete highcharter. Sin embargo, si el proceso comienza en algún nivel diferente al valor objetivo, la ventaja permitirá que el CUSUM lo detecte más rápidamente, lo que resulta en valores ARL fuera de control más cortos.

La columna (c) de la Tabla 5.6 presenta el desempeño de la \(\textit{LMC}\) del CUSUM básico con la característica de ventaja o FIR. Las \(\textit{LMC}\) se calcularon utilizando una ventaja anticipada del 50%. Observa que los valores \(\textit{LMC}\) para la CUSUM FIR son válidos para el caso en que el proceso está fuera de control en el momento en que se reinician las CUSUM. Cuando el proceso está bajo control, el valor de la ventaja disminuye rápidamente a cero. Por tanto, si el proceso está bajo control cuando la CUSUM se reinicializa pero más tarde se corre a la condición de fuera de control, la \(\textit{LMC}\) más apropiada para este caso deberá leerse de la columna (a), es decir, la CUSUM sin la característica FIR.

Código 
```{r}
#| label: fig-cusum-FIR-highcharter
#| fig-cap: "Diagrama CUSUM del status con la ventaja anticipada (FIR) con el paquete `highcharter`"

data_cusum_FIR <- dplyr::bind_rows(data_cusum_FIR_0, data_cusum_FIR_1) |> 
  dplyr::mutate(muestra = 1:20)

color <- ifelse(
  data_cusum_FIR$c_mas > 5 | -data_cusum_FIR$c_menos < -5,
  "#FF0000", "#0000FF"
)

# Diagrama de control CUSUM FIR

highchart() |>
  hc_add_series(
    data = data_cusum_FIR,
    type = "column",
    hcaes(
      x = muestra, y = c_mas, color = color
    ),
    showInLegend = TRUE, lineWidth = 3, useHTML = TRUE,
    name = "<p><em><strong>C<sub>i</sub><sup>+</sup></strong></em></p>",
    marker = list(enabled = TRUE),
    tooltip = list(
      valueDecimals = 2, useHTML = TRUE, crosshairs = TRUE,
      shared = TRUE, headerFormat = "<b><i>{point.x}</i></b><br>",
      pointFormat = "<p><em><strong>C<sub>i</sub><sup>+</sup></strong></em></p>: {point.y}"
    ),
    min = min(data_cusum_FIR$c_mas),
    max = max(data_cusum_FIR$c_mas),
    plotOptions = list(stacking = "normal")
  ) |>
  hc_add_series(
    data = data_cusum_FIR,
    type = "column",
    hcaes(
      x = muestra, y = -c_menos
    ),
    color = "blue",
    showInLegend = TRUE, lineWidth = 3, useHTML = TRUE,
    name = "<p><em><strong>C<sub>i</sub><sup>-</sup></strong></em></p>",
    marker = list(enabled = TRUE),
    tooltip = list(
      valueDecimals = 2, useHTML = TRUE, crosshairs = TRUE,
      shared = TRUE, headerFormat = "<b><i>{point.x}</i></b><br>",
      pointFormat = "<p><em><strong>C<sub>i</sub><sup>-</sup></strong></em></p>: {point.y}"
    ),
    plotOptions = list(stacking = "normal")
  ) |>
  # hc_plotOptions(series = list(stacking = "normal")) |>
  hc_add_series(
    data = data_cusum_FIR,
    hcaes(x = muestra, low = 0, high = 5),
    marker = list(enabled = FALSE, visible = FALSE),
    type = "arearange", fillOpacity = 0.2, 
    name = "", color = "green",
    showInLegend = FALSE, tooltip = list(pointFormat = "{NULL}")
  ) |>
  hc_add_series(
    data = data_cusum_FIR,
    hcaes(x = muestra, low = -5, high = 0),
    marker = list(enabled = FALSE, visible = FALSE),
    type = "arearange", fillOpacity = 0.2, 
    name = "", color = "blue",
    showInLegend = FALSE, tooltip = list(pointFormat = "{NULL}")
  ) |>
  hc_add_series(
    data = data_cusum_FIR,
    type = "line", hcaes(x = muestra, y = 5), useHTML = TRUE,
    color = "red", dashStyle = "Dash", name = "<b><i>H</i></b>",
    step = "hv", marker = list(enabled = FALSE),
    tooltip = list(valueDecimals = 0)
  ) |>
  hc_add_series(
    data = data_cusum_FIR,
    type = "line", hcaes(x = muestra, y = -5), useHTML = TRUE,
    color = "red", dashStyle = "Dash", name = "<b><i>-H</i></b>",
    step = "hv", marker = list(enabled = FALSE),
    tooltip = list(valueDecimals = 0)
  ) |>
  hc_add_series(
    data = data_cusum_FIR,
    type = "line", hcaes(x = muestra, y = 2.5), useHTML = TRUE,
    color = "black", dashStyle = "Dash", name = "<b><i>K</i></b>",
    step = "hv", marker = list(enabled = FALSE),
    tooltip = list(valueDecimals = 1)
  ) |>
  hc_add_series(
    data = data_cusum_FIR,
    type = "line", hcaes(x = muestra, y = -2.5), 
    useHTML = TRUE, color = "black", dashStyle = "Dash", 
    name = "<b><i>-K</i></b>",
    step = "hv", marker = list(enabled = FALSE),
    tooltip = list(valueDecimals = 1)
  ) |>
  hc_yAxis_multiples(
    list(
      title = list(
        text = "<p><strong>⟵</strong><em><strong>C<sub>i</sub><sup>-</sup> &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; C<sub>i</sub><sup>+</sup></strong></em><strong>⟶</strong>  &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;</p>"#,
        # align = "center", y = 0, margin = 0
      ),
      style = list(fontWeight = "bold", fontSize = "auto"),
      labels = list(
        format = "{value}",
        style = list(fontSize = "auto", color = "black")
      ),
      opposite = FALSE,
      min = -8,
      max = 8,
      lineWidth = 2
    ),
    list(
      title = list(text = "<p><em><strong>x<sub>i</sub></strong></em></p>"),
      style = list(fontWeight = "bold", fontSize = "auto"),
      labels = list(
        format = "{value}"
      ),
      min = min(data_cusum_FIR$x_i),
      max = max(data_cusum_FIR$x_i),
      opposite = FALSE, 
      lineWidth = 2
    )
  ) |>
  hc_add_series(
    data = data_cusum_FIR,
    hcaes(x = muestra, y = x_i),
    showInLegend = TRUE, type = "scatter", useHTML = TRUE,
    name = "<p><em><strong>x<sub>i</sub></strong></em></p>",
    marker = list(
      enabled = TRUE, radius = 3, lineColor = "black"
    ),
    tooltip = list(
      valueDecimals = 2, useHTML = TRUE, crosshairs = TRUE,
      shared = TRUE, headerFormat = "<b><i>{point.x}</i></b><br>",
      pointFormat = "<p><em><strong>x<sub>i</sub></strong></em></p>: {point.y}"
    ),
    yAxis = 1, xAxis = 0,
    plotOptions = list(stacking = "normal")
  ) |>
  hc_tooltip(split = TRUE) |>
  hc_xAxis(
    title = list(
      useHTML = TRUE,
      text = "<b><i>Periodo</i></b><br>"
    ),
    plotLines = list(
      list(color = "blue", value = 10.5, width = 2, dashStyle = "Dash")
    )
  ) |> 
  hc_plotOptions(
    area = list(stacking = "normal"),
    column = list(stacking = "normal"),
    line = list(
      stacking = "normal",
      color = "blue",
      marker = list(
        stacking = "normal",
        fillColor = color,
        lineWidth = 2,
        lineColor = "blue"
      )
    )
  )
```
Figura 5.10: Diagrama CUSUM del status con la ventaja anticipada (FIR) con el paquete highcharter
Código 
```{r}
#| label: fig-cusum-basico
#| fig-cap: "Diagrama CUSUM del status básico con el paquete `highcharter`"

# Diagrama de control CUSUM básico
data_cusum2 <- tibble::tibble(
  `muestra` = 1:20,
  `x_i` = c(data_cusum_FIR_0$x_i, data_cusum_FIR_1$x_i)
)
data_cusum2 <- data_cusum2 |>
  dplyr::mutate(
    col_b = x_i - (10 + 1 / 2),
    c_mas = cusum_stats(x = x_i, mu_0 = 10, K = 1 / 2)[[1]],
    n_mas = cusum_stats(x = x_i, mu_0 = 10, K = 1 / 2)[[2]],
    col_c = (10 - 1 / 2) - x_i,
    c_menos = cusum_stats(
      x = x_i, mu_0 = 10, K = 1 / 2
    )[[3]],
    n_menos = cusum_stats(x = x_i, mu_0 = 10, K = 1 / 2)[[4]]
  ) |>
  data.table::data.table()

color <- ifelse(
  data_cusum2$c_mas > 5 | -data_cusum2$c_menos < -5,
  "#FF0000", "#0000FF"
)

highchart() |>
  hc_add_series(
    data = data_cusum2,
    type = "column",
    hcaes(
      x = muestra, y = c_mas, color = color
    ),
    showInLegend = TRUE, lineWidth = 3, useHTML = TRUE,
    name = "<p><em><strong>C<sub>i</sub><sup>+</sup></strong></em></p>",
    marker = list(enabled = TRUE),
    tooltip = list(
      valueDecimals = 2, useHTML = TRUE, crosshairs = TRUE,
      shared = TRUE, headerFormat = "<b><i>{point.x}</i></b><br>",
      pointFormat = "<p><em><strong>C<sub>i</sub><sup>+</sup></strong></em></p>: {point.y}"
    ),
    min = min(data_cusum2$c_mas),
    max = max(data_cusum2$c_mas),
    plotOptions = list(stacking = "normal")
  ) |>
  hc_add_series(
    data = data_cusum2,
    type = "column",
    hcaes(
      x = muestra, y = -c_menos
    ),
    color = "blue",
    showInLegend = TRUE, lineWidth = 3, useHTML = TRUE,
    name = "<p><em><strong>C<sub>i</sub><sup>-</sup></strong></em></p>",
    marker = list(enabled = TRUE),
    tooltip = list(
      valueDecimals = 2, useHTML = TRUE, crosshairs = TRUE,
      shared = TRUE, headerFormat = "<b><i>{point.x}</i></b><br>",
      pointFormat = "<p><em><strong>C<sub>i</sub><sup>-</sup></strong></em></p>: {point.y}"
    ),
    plotOptions = list(stacking = "normal")
  ) |>
  # hc_plotOptions(series = list(stacking = "normal")) |>
  hc_add_series(
    data = data_cusum2,
    hcaes(x = muestra, low = 0, high = 5),
    marker = list(enabled = FALSE, visible = FALSE),
    type = "arearange", fillOpacity = 0.2, 
    name = "", color = "green",
    showInLegend = FALSE, 
    tooltip = list(pointFormat = "{NULL}")
  ) |>
  hc_add_series(
    data = data_cusum2,
    hcaes(x = muestra, low = -5, high = 0),
    marker = list(enabled = FALSE, visible = FALSE),
    type = "arearange", fillOpacity = 0.2, 
    name = "", color = "blue",
    showInLegend = FALSE, tooltip = list(pointFormat = "{NULL}")
  ) |>
  hc_add_series(
    data = data_cusum2,
    type = "line", hcaes(x = muestra, y = 5), useHTML = TRUE,
    color = "red", dashStyle = "Dash", name = "<b><i>H</i></b>",
    step = "hv", marker = list(enabled = FALSE),
    tooltip = list(valueDecimals = 0)
  ) |>
  hc_add_series(
    data = data_cusum2,
    type = "line", hcaes(x = muestra, y = -5), 
    useHTML = TRUE,
    color = "red", dashStyle = "Dash", name = "<b><i>-H</i></b>",
    step = "hv", marker = list(enabled = FALSE),
    tooltip = list(valueDecimals = 0)
  ) |>
  hc_add_series(
    data = data_cusum2,
    type = "line", hcaes(x = muestra, y = 2.5), 
    useHTML = TRUE,
    color = "black", dashStyle = "Dash", name = "<b><i>K</i></b>",
    step = "hv", marker = list(enabled = FALSE),
    tooltip = list(valueDecimals = 1)
  ) |>
  hc_add_series(
    data = data_cusum2,
    type = "line", hcaes(x = muestra, y = -2.5), 
    useHTML = TRUE,
    color = "black", dashStyle = "Dash", name = "<b><i>-K</i></b>",
    step = "hv", marker = list(enabled = FALSE),
    tooltip = list(valueDecimals = 1)
  ) |>
  hc_yAxis_multiples(
    list(
      title = list(
        text = "<p><strong>⟵</strong><em><strong>C<sub>i</sub><sup>-</sup> &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; C<sub>i</sub><sup>+</sup></strong></em><strong>⟶</strong> &nbsp; &nbsp; &nbsp;</p>"#,
        # align = "center", y = 0, margin = 0
      ),
      style = list(fontWeight = "bold", fontSize = "auto"),
      labels = list(
        format = "{value}",
        style = list(fontSize = "auto", color = "black")
      ),
      opposite = FALSE,
      min = -8,
      max = 8,
      lineWidth = 2
    ),
    list(
      title = list(text = "<p><em><strong>x<sub>i</sub></strong></em></p>"),
      style = list(fontWeight = "bold", fontSize = "auto"),
      labels = list(
        format = "{value}"
      ),
      min = min(data_cusum2$x_i),
      max = max(data_cusum2$x_i),
      opposite = FALSE,
      lineWidth = 2
    )
  ) |>
  hc_add_series(
    data = data_cusum2,
    hcaes(x = muestra, y = x_i),
    showInLegend = TRUE, type = "scatter", useHTML = TRUE,
    name = "<p><em><strong>x<sub>i</sub></strong></em></p>",
    marker = list(
      enabled = TRUE, radius = 3, lineColor = "black"
    ),
    tooltip = list(
      valueDecimals = 2, useHTML = TRUE, crosshairs = TRUE,
      shared = TRUE, headerFormat = "<b><i>{point.x}</i></b><br>",
      pointFormat = "<p><em><strong>x<sub>i</sub></strong></em></p>: {point.y}"
    ),
    yAxis = 1, xAxis = 0,
    plotOptions = list(stacking = "normal")
  ) |>
  hc_tooltip(split = TRUE) |>
  hc_xAxis(
    title = list(
      useHTML = TRUE,
      text = "<b><i>Periodo</i></b><br>"
    ),
    plotLines = list(
      list(
        color = "blue", value = 10.5, 
        width = 2, dashStyle = "Dash"
      )
    )
  ) |> 
  hc_plotOptions(
    area = list(stacking = "normal"),
    column = list(stacking = "normal"),
    line = list(
      stacking = "normal",
      color = "blue",
      marker = list(
        stacking = "normal",
        fillColor = color,
        lineWidth = 2,
        lineColor = "blue"
      )
    )
  ) # |>
  # hc_add_dependency("plugins/multicolor_series.js")
```
Figura 5.11: Diagrama CUSUM del status básico con el paquete highcharter

5.1.7 Una CUSUM para Monitorear la Variabilidad del Proceso

También pueden construirse diagramas de control CUSUM para monitorear la variabilidad del proceso. Puesto que la CUSUM se emplea generalmente con observaciones individuales, el procedimiento creado por D. M. Hawkins (1981) es de gran utilidad. El procedimiento es el siguiente, sea \(X_i\) la medición de una característica de la calidad del proceso que tiene una distribución normal con media o valor objetivo \(\mu_0\) y desviación estándar \(\sigma\) y \(Z_i\) la estandarización de \(X_i\) tal como se definió en la ecuación 5.8. Hawkins(1981, 1993) sugiere crear la nueva cantidad estandarizada

CUSUM para Monitorear la Variabilidad

\[ \begin{equation} V_i=\frac{\sqrt{\left| Z_i \right| -0,822}}{0,349} \end{equation} \tag{5.10}\]

Hawkins sugiere que las \(V_i\) son más sensitivas a los corrimientos en la varianza que a los de la media. De hecho, el estadístico \(V_i\) es sensitivo a los corrimientos tanto de la media como de la varianza. Puesto que la distribución bajo control de \(V_i\) es aproximadamente \(N(0,1)\), es posible establecer dos CUSUMS unilaterales a escala estandarizada (es decir, la desviación estándar) como sigue.

La CUSUM a escala

\[ \begin{align*} S_i^+ & = \mathrm{máx}\left\{0, \, V_i - k + S_{i - 1}^+ \right\}\\ S_i^- & = \mathrm{máx}\left\{0, \, - k - V_i + C_{i - 1}^ - \right\}\\ \end{align*}\:. \tag{5.11}\]

Donde \(S_0^+ = S_i^- = 0\) (a menos que se use la característica FIR) y los valores de \(k\) y \(h\) se seleccionan como en la CUSUM para controlar la media del proceso.

La interpretación de la CUSUM a escala es similar a la interpretación de la CUSUM para la media. Si la desviación estándar del proceso se incrementa, los valores de \(S_i^+\) se incrementarán y terminarán por exceder \(h\), mientras que si la desviación estándar se reduce, los valores \(S_i^-\) se incrementarán y terminarán por exceder a \(h\).

Aun cuando podría establecerse diagramas de control CUSUM del status separadas para la media y la varianza D. M. Hawkins (1993) sugiere trazarla en la misma gráfica. También ofrece varios ejemplos excelentes y una discusión más amplia de este procedimiento. El estudio de sus ejemplos será valioso a fin de mejorar la habilidad del lector para determinar corrimientos en la variabilidad del proceso a partir de la CUSUM a escala. Si la CUSUM a escala envía una señal, se presumirá un cambio en la varianza, pero si ambas CUSUM envían una señal, se presumiría un corrimiento en la media.

El prosedimiento CUSUM máscara V

Barnard (1959) propuso un procedimiento alternativo al CUSUM tabular, denominado procedimiento máscara V, el cual no se tratará en este libro debido a las recomendaciones de D. C. Montgomery (2013, 430).

5.2 Diagrama de Control de Promedio Móvil Ponderado Exponencialmente

El diagrama de control del promedio móvil ponderado exponencialmente (o EWMA) fue propuesta por Roberts (1959). Este diagrama es una buena alternativa para la carta de control de Shewhart cuando se está interesado en detectar corrimientos pequeños. El desempeño del diagrama de control EWMA es aproximadamente equivalente al del diagrama de control CUSUM y, en cierto modo, es más sencilla de establecer y operar. Como con la CUSUM, el EWMA se usa de manera típica con observaciones individuales y, por consiguiente, se discute primero este caso. Se presentará también los resultados para subgrupos racionales de tamaño \(n > 1\).

En este diagrama el estadístico que se grafica es el promedio móvil ponderado exponencialmente (EWMA), el cual se define como

Estadístico EWMA recursivo

\[ \begin{equation} Z_{i} = \lambda X_{i} + \left( 1 - \lambda \right) Z_{i - 1} \end{equation} \tag{5.12}\]

donde \(0< \lambda ≤1\) es una constante y el valor inicial (requerido con la primera muestra en \(i = 1\)) es el objetivo del proceso, es decir, \(Z_0 = \mu_0\).

En ocasiones se usa el promedio de datos preliminares como valor inicial del EWMA, de tal modo que \(Z_0 = \bar{X}\).

Para demostrar que el EWMA \(Z_i\) es un promedio ponderado de todas las observaciones anteriores, puede sustituirse

\[ Z_{i - 1} = \lambda X_{i - 1} + \left( 1 - \lambda \right) Z_{i - 2} \]

en el segundo miembro de la ecuación 5.12 para obtener

\[ \begin{align*} Z_{i} & = \lambda X_{i} + \left( 1 - \lambda \right) \left[ \lambda X_{i - 1} + \left( 1 - \lambda \right) Z_{i - 2} \right] \\ & = \lambda X_{i} + \lambda \left( 1 - \lambda \right) X_{i - 1} + (1 - \lambda)^{2} Z_{i - 2} \end{align*} \]

Al continuar las sustituciones recursivamente para \(Z_{i - j}\), \(j = 2, 3, \dots, t\); se obtiene

\[ \begin{equation} Z_{i} = \lambda \sum_{j = 0}^{i - 1}\left( 1 - \lambda \right)^{j} X_{i-j} + \left( 1 - \lambda \right)^{i} Z_{0} \end{equation} \tag{5.13}\]

Las ponderaciones \(\lambda \left( 1 - \lambda \right)^{j}\) disminuyen geométricamente con la antigüedad de la media muestral. Es claro que el parámetro \(\lambda\) determina la profundidad de la memoria de la EWMA: mientras más cerca esté de cero es mayor el peso de los datos históricos, es decir, recuerda más el pasado. Mientras que, si está más cerca de uno, tiene más influencia la última observación y el pasado tiene menos peso. De tal forma que cuando \(\lambda = 1\) sería equivalente al diagrama de observaciones individuales tradicional, que no da ningún peso la información anterior a un punto dado. La experiencia ha mostrado que lo adecuado es que \(0,1\le \lambda \le 0,3\), y el valor 0.2 es el más típico. Además, la suma de las ponderaciones es la unidad, ya que

\[ \lambda \sum_{j = 0}^{i - 1}\left( 1 - \lambda \right)^{j} = \lambda\left[ \frac{1 - \left( 1 - \lambda \right)^{i}}{1 - \left( 1 - \lambda \right)} \right] = 1 -\left( 1 - \lambda \right)^{i} \]

Si \(\lambda = 0,2\), entonces las ponderaciones asignadas a la media muestral actual es 0,2 y las ponderaciones dadas a las medias anteriores son 0,16, 0,128, 0,1024, y así sucesivamente. En la Figura 5.12 se muestra una comparación de estas ponderaciones con las de un promedio móvil de cinco periodos. Debido a que estas ponderaciones declinan geométricamente cuando se unen con una curva suave, al EWMA se le llama ocasionalmente promedio móvil geométrico (GMA, por sus siglas en inglés).

Código 
```{r}
#| label: fig-ewma-peso
#| fig-cap: "Ponderación de las muestras en el tiempo para $\\lambda = 0,2$"

pesos <- data.frame(
  i = 0:10,
  peso = 0.2 * (1 - 0.2)^(0:10)
)

highchart() |>
  hc_add_series(
    data = pesos,
    type = "line", color = "black",
    hcaes(x = i, y = peso),
    marker = list(radius = 3),
    useHTML = TRUE, name = "<b><i>Peso</i></b>",
    dashStyle = "Dash",
    tooltip = list(valueDecimals = 4)
  ) |>
  hc_add_series(
    data = data.frame(
      i = c(0:4, rep(NA, 6)),
      peso = c(rep(0.2, 5), rep(NA, 6))
    ),
    hcaes(x = i, y = peso),
    marker = list(enabled = FALSE, visible = FALSE),
    type = "area", fillOpacity = 0.2, name = "", 
    color = "blue",
    showInLegend = FALSE, 
    tooltip = list(pointFormat = "{NULL}")
  ) |>
  hc_xAxis(
    title = list(
      useHTML = TRUE,
      text = "<b><i>Antiguedad de la media muestral</i></b><br>"
    )
  ) |>
  hc_yAxis(
    title = list(
      useHTML = TRUE,
      text = "<b><i>Ponderación</i></b><br>"
    )
  ) |>
  hc_tooltip(crosshairs = TRUE, shared = TRUE) |> 
  hc_plotOptions(series = list(animation = FALSE) )
```
Figura 5.12: Ponderación de las muestras en el tiempo para \(\lambda = 0,2\)

Puesto que el EWMA pude considerarse como un promedio ponderado de todas las observaciones pasadas y actuales, es muy insensible al supuesto de normalidad. Por lo tanto, es un diagrama de control ideal para usarse con observaciones individuales.

Si las observaciones \(X_i\) son variables aleatorias independientes con media \(\mu_0\) y varianza \(\sigma^{2}\), entonces la media de \(Z_i\) viene dado por

\[ \begin{align*} \mu_{_{Z_i}} & = E\left( Z_i \right) = E\left[ \lambda \sum_{j = 0}^{i - 1}\left( 1 - \lambda \right)^{j} X_{i-j} + \left( 1 - \lambda \right)^{i} Z_{0} \right] \\ & = \lambda \sum_{j = 0}^{i - 1} \left( 1 - \lambda \right)^{j} E\left( X_{i- j}\right) + \left( 1 - \lambda \right)^{i} E\left( \mu_{0} \right)\\ & = \lambda \sum_{j = 0}^{i - 1}\left( 1 - \lambda \right)^{j} \mu_{0} + \left( 1 - \lambda \right)^{i} \mu_{0} \\ & = \mu_{0}\left[ \lambda \sum_{j = 0}^{i - 1} \left( 1 - \lambda \right)^{j} + \left( 1 - \lambda \right)^{i}\right] \end{align*} \:. \]

Como se demostró anteriormente que

\[ \lambda\sum_{j = 0}^{i - 1}\left( 1 - \lambda \right)^{j} = 1 - \left( 1 - \lambda \right)^{i} \:. \]

Entonces,

\[ \mu_{_{Z_i}} = \mu_{0}\left[ 1 - \left( 1 - \lambda \right)^{i} + \left( 1 - \lambda \right)^{i}\right] = \mu_{0} \:. \]

Mientras que la varianza de \(Z_i\) se obtiene de la siguiente manera

\[ \begin{align*} \sigma_{_{\!Z_i}}^{2} & = V\left( Z_i \right) = V\left[ \lambda\sum_{j = 0}^{i - 1}\left( 1 - \lambda \right)^{j} X_{i-j} + \left( 1 - \lambda \right)^{i} Z_{0}\right] \\ & = \lambda^{2} \sum_{j = 0}^{i - 1}\left( 1 - \lambda \right)^{2j} V\left( X_{i - j} \right) + \left( 1 - \lambda \right)^{2i}V\left( Z_{0} \right) \\ & = \sigma^{2}\left[ \lambda^{2}\sum_{j = 0}^{i - 1} \left( 1 - \lambda \right)^{2j}\right] \end{align*} \:. \]

Dado que la suma anterior es una progresión geométrica que tiene como primer término uno y razón \((1 - \lambda)^2\), entonces la varianza puede ser reescrita como

\[ \begin{align*} \sigma_{_{\!Z_i}}^{2} & = \sigma^{2}\left[ \lambda^{2} \left( \frac{1 - \left( 1 - \lambda \right)^{2j}}{1 - \left( 1 -\lambda \right)^{2}} \right) \right] \\ & = \sigma^{2}\left[ \lambda^{2} \left( \frac{1 - \left( 1 - \lambda \right)^{2j}}{1 - \left( 1 -2\lambda + \lambda^{2} \right)} \right) \right] \\ & = \sigma^{2}\left[ \lambda^{2} \left( \frac{1 - \left( 1 - \lambda \right)^{2j}}{2\lambda - \lambda^{2} } \right) \right] \\ & = \sigma^{2}\left[ \lambda^{2} \left( \frac{1 - \left( 1 - \lambda \right)^{2j}}{\lambda\left( 2 - \lambda \right) } \right) \right] \\ & = \sigma^{2}\left( \frac{\lambda}{2 - \lambda} \right)\left[ 1 - \left( 1 - \lambda \right)^{2i} \right] \end{align*} \:. \]

En resumen, la media y la varianza de \(Z_i\) vienen dadas por

Media y varianza de \(Z_i\)

\[ \begin{align*} \mu_{_{Z_i}} & = \mu_{0} \\ \sigma_{_{\!Z_i}}^{2} & = \sigma^{2}\left( \frac{\lambda}{2 - \lambda} \right)\left[ 1 - \left( 1 - \lambda \right)^{2i} \right] \end{align*} \:. \tag{5.14}\]

Por lo tanto, el diagrama de control EWMA se construye graficando \(Z_i\) contra el número de muestra \(i\) (o el tiempo). La línea central y los límites de control, tipo Shewhart, para el diagrama de control EWMA son los siguientes

Límites de control estándar del diagrama EWMA

\[ \begin{align*} LCS & = \mu_{0} + L\sigma\sqrt{\left( \frac{\lambda}{2 - \lambda} \right)\left[ 1 - \left( 1 - \lambda \right)^{2i} \right]} \\ LC & = \mu_{0} \\ LCI & = \mu_{0} - L\sigma\sqrt{\left( \frac{\lambda}{2 - \lambda} \right)\left[ 1 - \left( 1 - \lambda \right)^{2i} \right]} \\ \end{align*} \:. \tag{5.15}\]

En la ecuación 5.15, el factor \(L\) es la anchura de los límites de control. Se discutirá brevemente la elección de los parámetros \(L\) y \(\lambda\).

Obsérvese que el término \(1 - \left( 1 - \lambda \right)^{2i}\) en la ecuación 5.15 tiende a la unidad cuando \(i\) se hace grande. Esto significa que después que el diagrama de control EWMA ha estado operando durante varios periodos de tiempo, los límites de control tenderán a valores de estado estable o régimen permanente dados por

\[ \begin{split} LCS & = \mu_{0} + L\sigma\sqrt{\left( \frac{\lambda}{2 - \lambda} \right)} \\ LC & = \mu_{0} \\ LIC & = \mu_{0} - L\sigma\sqrt{\left( \frac{\lambda}{2 - \lambda} \right)} \\ \end{split} \:. \tag{5.16}\]

Sin embargo se recomienda enfáticamente usar los límites de control exactos dados en la ecuación ecuación 5.15 para valores pequeños de \(i\). Con esto se mejorará de manera significativa el desempeño de la carta de control para detectar un proceso fuera del objetivo inmediatamente después de que se ha iniciado el EWMA

Ejemplo 5.4 (Aplicación del Diagrama EWMA para Observaciones Individuales)
Se aplicará el diagrama de control EWMA con \(\lambda = 0,\!1\) y \(L = 2,\!7\) a los datos de la Tabla 5.1. Recuerde que el valor objetivo de la media es \(\mu_0 = 100\) y que la desviación estándar es \(\sigma = 5\). En la Tabla 5.9 se muestran los cálculos de los EWMA y los límites de control del diagrama. Mientras que la Figura 5.15 muestra el diagrama de control para estos datos.

Para ilustrar los cálculos, considérese la primera observación, \(x_1 = 96,9663\). El primer valor del EWMA es

\[ z_{1} = \lambda x_{1} + \left( 1 - \lambda \right) z_{0} = 0,1 * \left(96,9663\right) + \left( 1 - 0,1 \right) * \left( 100 \right) = 99,6966 \:. \]

Por lo tanto, \(z_1 = 99,6966\) es el primer punto graficado en la Figura 5.15. El segundo valor del EWMA es

\[ z_{2} = \lambda x_{2} + \left( 1 - \lambda \right) z_{1} = 0,1 * \left(97,2494\right) + \left( 1 - 0,1 \right) * \left( 99,6966 \right) = 99,4519 \:. \]

Los demás valores del estadístico EWMA se calculan de mantera similar.

Código 
```{r}
#| label: tbl-ewma
#| tbl-cap: "Cálculo del EWMA para los datos de la @tbl-datos-cusum"
#| table-caption: left
#| class: output

set.seed(347)
data_ewma <- tibble::tibble(
  `muestra` = 1:30,
  `x_i` = c(
    sample(
      rnorm(n = 10000, mean = 100, sd = 5),
      size = 20, replace = TRUE
    ),
    sample(
      rnorm(n = 10000, mean = 105, sd = 5),
      size = 10, replace = TRUE
    )
  )
)

data_ewma <- data_ewma |>
  dplyr::mutate(
    ewma = qcc::ewma(
      data = data_ewma$x_i, sizes = 1, center = 100, std.dev = 5, 
      lambda = 0.1, nsigmas = 2.7, plot = FALSE
    )$y
  ) |>
  dplyr::bind_cols(
    qcc::ewma(
      data = data_ewma$x_i, sizes = 1, center = 100, std.dev = 5, 
      lambda = 0.1, nsigmas = 2.7, plot = FALSE
    )$limits
  ) |> 
  data.table::data.table()

knitr::kable(
  data_ewma,
  format = "markdown",
  digits = 4,
  format.args = list(decimal.mark = ",", big.mark = "."),
  col.names = c(
    "$\\textbf{Periodo} \\: \\left( i \\right)$", "$x_i$", 
    "$\\textbf{EWMA} \\: \\left( z_i \\right)$",
    "$LCI$", "$LCS$"
  ),
  align = "ccccc",
  escape = FALSE
)
```
Tabla 5.9: Cálculo del EWMA para los datos de la Tabla 5.1
\(\textbf{Periodo} \: \left( i \right)\) \(x_i\) \(\textbf{EWMA} \: \left( z_i \right)\) \(LCI\) \(LCS\)
1 96,9663 99,6966 98,6500 101,3500
2 97,2494 99,4519 98,1838 101,8162
3 103,6353 99,8702 97,8800 102,1200
4 96,7796 99,5612 97,6627 102,3373
5 110,3248 100,6375 97,5005 102,4995
6 97,7647 100,3503 97,3765 102,6235
7 97,4743 100,0627 97,2801 102,7199
8 100,8277 100,1392 97,2045 102,7955
9 91,8081 99,3061 97,1448 102,8552
10 107,2495 100,1004 97,0973 102,9027
11 94,6188 99,5522 97,0593 102,9407
12 105,6507 100,1621 97,0290 102,9710
13 94,9607 99,6420 97,0046 102,9954
14 107,0322 100,3810 96,9850 103,0150
15 102,3922 100,5821 96,9692 103,0308
16 100,0102 100,5249 96,9565 103,0435
17 96,8418 100,1566 96,9463 103,0537
18 97,6166 99,9026 96,9380 103,0620
19 97,9844 99,7108 96,9313 103,0687
20 99,4052 99,6802 96,9259 103,0741
21 110,7501 100,7872 96,9215 103,0785
22 102,4706 100,9555 96,9179 103,0821
23 108,5581 101,7158 96,9151 103,0849
24 113,1497 102,8592 96,9128 103,0872
25 112,3129 103,8046 96,9109 103,0891
26 107,0480 104,1289 96,9094 103,0906
27 110,5134 104,7673 96,9081 103,0919
28 95,7660 103,8672 96,9071 103,0929
29 107,5611 104,2366 96,9063 103,0937
30 98,1360 103,6265 96,9057 103,0943


Sustituyendo \(\mu_0 = 100\), \(\sigma = 5\) y \(L = 2,7\) en la ecuación 5.15 se obstinen los límites de control y la línea central del diagrama, como se muestra a continuación

\[ \begin{align*} LCS & = 100 + 2,7\left(5 \right)\sqrt{\left( \frac{0,1}{2 - 0,1} \right)\left[ 1 - \left( 1 - 0,1\right)^{2i} \right]} \\ LC & = 100 \\ LCI & = 100 - 2,7\left(5 \right)\sqrt{\left( \frac{0,1}{2 - 0,1} \right)\left[ 1 - \left( 1 - 0,1\right)^{2i} \right]} \\ \end{align*} \:. \]

Haciendo \(i = 1\) en las ecuaciones anteriores se obtienen los límites de control para el primer periodo, como sigue

\[ \begin{align*} {LCS}_1 & = 100 + 2,7\left(5 \right)\sqrt{\left( \frac{0,1}{2 - 0,1} \right)\left[ 1 - \left( 1 - 0,1\right)^{2\left( 1\right)} \right]} = 101,35\\ {LC}_1 & = 100 \\ {LCI}_1 & = 100 - 2,7\left(5 \right)\sqrt{\left( \frac{0,1}{2 - 0,1} \right)\left[ 1 - \left( 1 - 0,1\right)^{2\left( 1\right)} \right]} = 98,65\\ \end{align*} \:. \]

Código 
```{r}
#| label: fig-ewma-ggplot2
#| fig-cap: "Diagrama EWMA para los datos de @tbl-datos-cusum con el paquete `ggplot2`"

q <-  qcc::ewma(
  data = data_ewma$x_i, sizes = 1, center = 100, std.dev = 5, 
  lambda = 0.1, nsigmas = 2.7, plot = FALSE
) 

# summary(q)

data_ewma$violacion <- ifelse(
  data_ewma$muestra == unique(q$violations), "Si", "No"
)

ggplot(data = data_ewma, aes(x = muestra, y = ewma), show.legend = FALSE) +
  geom_step(
    mapping = aes(x = muestra, y = LCL),
    direction = "hv", color = "red", lty = "dashed"
  ) +
  geom_step(
    mapping = aes(x = muestra, y = 100),
    direction = "hv", color = "blue"
  ) +
  geom_step(
    mapping = aes(x = muestra, y = UCL),
    direction = "hv", color = "red", lty = "dashed"
  ) +
  geom_point(
    mapping = aes(
      x = muestra, y = ewma, color = factor(violacion), 
      shape = factor(violacion)
    ),
    size = 2
  ) +
  scale_color_manual(
    values = c("No" = "black", "Si" = "red")
  ) +
  scale_shape_manual(values = c(16, 8)) + 
  geom_line(
    mapping = aes(x = muestra, y = ewma), color = "black"
  ) +
  pammtools::geom_stepribbon(
    aes(ymin = LCL, ymax = 100), fill = "blue", alpha = 0.2
  ) +
  pammtools::geom_stepribbon(
    aes(ymin = 100, ymax = UCL), fill = "green", alpha = 0.2
  ) +
  annotate(
    "text",
    x = nrow(data_ewma) + 1, y = last(data_ewma$LCL), 
    label = "LCI", color = "red"
  ) +
  annotate(
    "text",
    x = nrow(data_ewma) + 1, y = last(data_ewma$UCL), 
    label = "LCS", color = "red"
  ) +
  annotate(
    "text",
    x = nrow(data_ewma) + 1, y = 100, label = "LC", 
    color = "blue"
  ) +
  labs(
    x = latex2exp::TeX(
      "$Periodo \\, \\left(i\\right)$", 
      bold = TRUE, italic = TRUE
    ),
    y = latex2exp::TeX(
      "$\\textbf{EWMA} \\, \\left( z_i \\right)$",
      bold = TRUE, italic = TRUE
    ),
    shape = "Fuera de control", color = "Fuera de control"
  ) +
  theme_bw() +
  theme(
    legend.position = "bottom",
    legend.text = element_text(face = "bold"),
    legend.title = element_text(face = "bold"),
    legend.title.position = "top"
  ) 
```
Figura 5.13: Diagrama EWMA para los datos de Tabla 5.1 con el paquete ggplot2
Código 
```{r}
#| label: fig-ewma-qcc
#| fig-cap: "Diagrama EWMA para los datos de @tbl-datos-cusum con el paquete `ggplot2`"

plot(
  x = q, add.stats = FALSE, label.limits = c("LCI", "LCS"),
  title = "",xlab = "Periodo", ylab = "EWMA", axes.las = 0, 
  digits = getOption("digits"), restore.par = TRUE
)
```
Figura 5.14: Diagrama EWMA para los datos de Tabla 5.1 con el paquete ggplot2
Código 
```{r}
#| label: fig-ewma-highcharter
#| fig-cap: "Diagrama EWMA con el paquete `highcharter`"

q <-  ewma(
  data = data_ewma$x_i, sizes = 1, center = 100, std.dev = 5, 
  lambda = 0.1, nsigmas = 2.7, plot = FALSE
) 

color <- ifelse(data_ewma$muestra == unique(q$violations), "#FF0000", "#000000")
symbol <- ifelse(data_ewma$muestra == unique(q$violations), "asterisk-open", "circle")

highchart() |>
  hc_add_series(
    data = data_ewma, hcaes(x = muestra, low = UCL, high = 100),
    step = "hv", marker = list(enabled = FALSE, visible = FALSE),
    type = "arearange", fillOpacity = 0.2, 
    name = "", color = "green",
    showInLegend = FALSE, tooltip = list(pointFormat = "{NULL}")
  ) |>
  hc_add_series(
    data = data_ewma, hcaes(x = muestra, low = LCL, high = 100),
    step = "hv", marker = list(enabled = FALSE, visible = FALSE),
    type = "arearange", fillOpacity = 0.2, 
    name = "", color = "blue",
    showInLegend = FALSE, tooltip = list(pointFormat = "{NULL}")
  ) |>
  hc_add_series(
    data_ewma, type = "line", 
    hcaes(x = muestra, y = UCL), useHTML = TRUE,
    color = "red", dashStyle = "Dash", name = "<b><i>LCS</i></b>",
    step = "hv", marker = list(enabled = FALSE),
    tooltip = list(valueDecimals = 4)
  ) |>
  hc_add_series(
    data_ewma,
    type = "line", hcaes(x = muestra, y = 100), useHTML = TRUE,
    color = "blue", name = "<b><i>LC</i></b>",
    step = "hv", marker = list(enabled = FALSE),
    tooltip = list(valueDecimals = 0)
  ) |>
  hc_add_series(
    data_ewma, type = "line", hcaes(x = muestra, y = LCL), 
    useHTML = TRUE, color = "red", dashStyle = "Dash", 
    name = "<b><i>LCI</i></b>",
    step = "hv", marker = list(enabled = FALSE),
    tooltip = list( valueDecimals = 4)
  ) |>
  hc_add_series(
    data_ewma, type = "line", color = "black",
    hcaes(
      x = muestra, y = ewma, segmentColor = color, 
      color = color, symbol = symbol
    ),
    marker = list(radius = 3, color = color, symbol = symbol),
    useHTML = TRUE, name = "<b><i>EWMA</i></b>",
    dashStyle = "Dash",
    tooltip = list( valueDecimals = 4)
  ) |>
  hc_xAxis(
    title = list(
      useHTML = TRUE,
      text = "<b><i>Periodo (i)</i></b><br>"
    )
  ) |>
  hc_yAxis(
    title = list(
      useHTML = TRUE,
      text = "<b><i>EWMA",
      # max = 0.5,
      # min = 0,
      startOnTick = FALSE
    )
  ) |>
  hc_tooltip(crosshairs = TRUE, shared = TRUE) |> 
  hc_plotOptions( series = list(animation = TRUE)) # |> 
  # hc_add_dependency("plugins/multicolor_series.js")
```
Figura 5.15: Diagrama EWMA con el paquete highcharter
Código 
```{r}
#| label: fig-ewma-plotly
#| fig-cap: "Diagrama EWMA con el paquete `plotly`"

plot_ly(data_ewma, x = ~muestra) |>
  add_lines(
    y = ~UCL, name = "<b><i>LCS</i></b>", 
    line = list(
      shape = "hvh", color = "red", dash = "dash"
    ),
    hovertemplate = paste(
      "<b><i>i = %{x}</i></b>",
      "<br><b><i>LCS<b><i> = %{y:.4f}<br><extra></extra>"
    )
  ) |>
  add_lines(
    y = 100, name = "<b><i>LC</i></b>", fill = "tonexty", 
    fillcolor = 'green',
    hovertemplate = paste(
      "<b><i>i = %{x}</i></b>",
      "<br><b><i>LC<b><i> = %{y:.0f}<br><extra></extra>"
    ),
    line = list(
      shape = "linear", color = "blue"
    )
  ) |>
  add_lines(
    y = ~LCL, name = "<b><i>LCI</i></b>", fill = "tonexty",
    fillcolor = 'blue',
    hovertemplate = paste(
      "<b><i>i = %{x}</i></b>",
      "<br><b><i>LCI<b><i> = %{y:.4f}<br><extra></extra>"
    ),
    line = list(
      shape = "hvh", color = "red", dash = "dash"
    )
  ) |>
  add_trace(
    y = ~ewma,
    name = "<b><i>EWMA<b><i>", type = "scatter",
    mode = "lines+markers",
    marker = list(color = color, symbol = symbol),
    hovertemplate = paste(
      "<b><i>i = %{x}</i></b>",
      "<br><b><i>z<sub>i</sub><b><i> = %{y:.4f}<br><extra></extra>"
    ),
    line = list(color = "black", dash = "dash")
  ) |>
  layout(
    legend = list(
      x = 0.5, y = -0.23, xref = 'paper', yref = 'paper',
      orientation = "h", xanchor = 'center', showarrow = F
    ),
    xaxis = list(
      title = "<b><i>Periodo</i></b><br>",
      zeroline = F
    ),
    yaxis = list(
      title = "<b><i>EWMA (z<sub>i</sub>)</i></b><br>",
      zeroline = FALSE
    )
  ) |>
  plotly::config(locale = "es", mathjax = "cdn") |>
  fillOpacity(alpha = 0.3)
```
Figura 5.16: Diagrama EWMA con el paquete plotly

5.2.1 Diseño del Diagramas de Control EWMA

El diagrama de control EWMA es muy efectivo para detectar corrimientos pequeños en el proceso. Los parámetros de diseño de este diagrama son el múltiplo de sigma usado en los límites de control (\(L\)) y el valor de \(\lambda\). Es posible elegir estos parámetros para obtener un desempeño de la \(\textit{LMC}\) del diagrama EWMA muy cercano al desempeño de la \(\textit{LMC}\) del diagrama CUSUM para detectar corrimientos pequeños.

Se han realizado diversos estudios teóricos de las propiedades de la longitud media de corrida del diagrama de control EWMA. Por ejemplo, las publicaciones de Crowder(1987, 1989) y Lucas y Saccucci (1990). En estos estudios se presentan tablas o gráficas de la longitud media de corrida para un rango de valores de \(\lambda\) y \(L\). En la Tabla 5.10 se muestra el desempeño de la longitud media de corrida para varios esquemas de control EWMA. El procedimiento de diseño óptimo consistirá en especificar las longitudes medias de corridas bajo control (\({\textit{LMC}}_0\)) y fuera de control (\({\textit{LMC}}_1\)), y la magnitud del corrimiento que se desea detectar, para después seleccionar los parámetros \(\lambda\) y \(L\) que proporcione el desempeño de la \(\textit{LMC}\) deseado.

Código 
```{r}
#| label: tbl-lmc-ewma
#| tbl-cap: "Longitud media de corrida para varios esquemas de control EWMA $[$adaptación de @Lucas1990$]$"
#| table-caption: left

LMC_Lucas <- data.frame(
  D = c(0, 0.25, 0.5, 0.75, 1, 1.5, 2, 2.5, 3, 4),
  L1 = c(500, 224, 71.2, 28.4, 14.3, 5.9, 3.5, 2.5, 2, 1.4),
  L2 = c(500, 170, 48.2, 20.1, 11.1, 5.5, 3.6, 2.7, 2.3, 1.7),
  L3 = c(500, 150, 41.8, 18.2, 10.5, 5.5, 3.7, 2.9, 2.4, 1.9),
  L4 = c(500, 106, 31.3, 15.9, 10.3, 6.1, 4.4, 3.4, 2.9, 2.2),
  L5 = c(500, 84.1, 28.8, 16.4, 11.4, 7.1, 5.2, 4.2, 3.5, 2.7)
) |> 
  data.table::as.data.table()

# options(
#   reactable.language = reactable::reactableLang(
#     searchPlaceholder = "buscar...",
#     pageNext = "Siguiente",
#     pagePrevious = "Anterior",
#     noData = "No entries found",
#     pageInfo = "{rowStart} de {rowEnd} de {rows} filas"
#   )
# )

# reactable::reactable(
#   LMC_Lucas,
#   # defaultColDef = colDef(
#   #   headerStyle = list(background = "#f7f7f8")
#   # ),
#   columns = list(
#     D = colDef(
#       name = "D", align = "right",
#       format = colFormat(
#         separators = TRUE, digits = 2, locales = "es-ES"
#       )
#     ),
#     L1 = colDef(
#       name = "L = 3,054", align = "right",
#       format = colFormat(
#         separators = TRUE, digits = 2, locales = "es-ES"
#       )
#     ),
#     L2 = colDef(
#       name = "L = 2,998", align = "right",
#       format = colFormat(
#         separators = TRUE, digits = 2, locales = "es-ES"
#       )
#     ),
#     L3 = colDef(
#       name = "L = 2,962", align = "right",
#       format = colFormat(
#         separators = TRUE, digits = 2, locales = "es-ES"
#       )
#     ),
#     L4 = colDef(
#       name = "L = 2,814", align = "right",
#       format = colFormat(
#         separators = TRUE, digits = 2, locales = "es-ES"
#       )
#     ),
#     L5 = colDef(
#       name = "L = 2,615", align = "right",
#       format = colFormat(
#         separators = TRUE, digits = 2, locales = "es-ES"
#       )
#     )
#   ),
#   striped = TRUE,
#   bordered = TRUE,
#   highlight = TRUE,
#   filterable = FALSE,
#   minRows = 10,
#   columnGroups = list(
#     colGroup(name = "Lambda = 0,40", columns = "L1"),
#     colGroup(name = "Lambda = 0,25", columns = "L2"),
#     colGroup(name = "Lambda = 0,20", columns = "L3"),
#     colGroup(name = "Lambda = 0,10", columns = "L4"),
#     colGroup(name = "Lambda = 0,05", columns = "L5")
#   )
# ) 

LMC_Lucas |> 
  gt() |> 
   cols_label(
    D = html(
      "<p><em><strong>Corrimiento en la media (m&uacute;ltiplo de &sigma;)</strong></em></p>"
    ),
    L1 = html("<p><em><strong>L = 3,054</strong></em></p>"),
    L2 = html("<p><em><strong>L = 2,998</strong></em></p>"),
    L3 = html("<p><em><strong>L = 2,962</strong></em></p>"),
    L4 = html("<p><em><strong>L = 2,814</strong></em></p>"),
    L5 = html("<p><em><strong>L = 2,615</strong></em></p>"),
  ) |> 
  tab_spanner(
    label = html(
      "<p><em><strong>&lambda; = 0,40</strong></em></p>"
    ),
    columns = c(L1)
  ) |> 
  tab_spanner(
    label =  html(
      "<p><em><strong>&lambda; = 0,25</strong></em></p>"
    ),
    columns = c(L2)
  ) |> 
  tab_spanner(
    label =  html(
      "<p><em><strong>&lambda; = 0,20</strong></em></p>"
    ),
    columns = c(L3)
  ) |> 
  tab_spanner(
    label =  html(
      "<p><em><strong>&lambda; = 0,10</strong></em></p>"
    ),
    columns = c(L4)
  ) |> 
  tab_spanner(
    label =  html(
      "<p><em><strong>&lambda; = 0,05</strong></em></p>"
    ),
    columns = c(L5)
  ) |> 
   tab_style(
    style = list(
      cell_text(weight = "bold")
      ),
    locations = cells_body(
      columns = D,
    )
  ) |> 
  fmt_number(decimals = 1, sep_mark = ".", dec_mark = ",") |> 
  opt_align_table_header(align = "left") |> 
  cols_align(
    align = "center",
    columns = everything()
  )
```
Tabla 5.10: Longitud media de corrida para varios esquemas de control EWMA \([\)adaptación de Lucas y Saccucci (1990)\(]\)

Corrimiento en la media (múltiplo de σ)

λ = 0,40

λ = 0,25

λ = 0,20

λ = 0,10

λ = 0,05

L = 3,054

L = 2,998

L = 2,962

L = 2,814

L = 2,615
0,0 500,0 500,0 500,0 500,0 500,0
0,2 224,0 170,0 150,0 106,0 84,1
0,5 71,2 48,2 41,8 31,3 28,8
0,8 28,4 20,1 18,2 15,9 16,4
1,0 14,3 11,1 10,5 10,3 11,4
1,5 5,9 5,5 5,5 6,1 7,1
2,0 3,5 3,6 3,7 4,4 5,2
2,5 2,5 2,7 2,9 3,4 4,2
3,0 2,0 2,3 2,4 2,9 3,5
4,0 1,4 1,7 1,9 2,2 2,7

Por otro lado, Bracho (2000), basado en el esquema de Crowder (1989), construyo unas tablas en las cuales se puede seleccionar los pares de valores \(\lambda\) y \(L\) para diferentes longitudes medias de corridas bajo control, específicamente para \({\textit{LMC}}_0 = 100\), \({\textit{LMC}}_0 = 250\), \({\textit{LMC}}_0 = 350\) y \({\textit{LMC}}_0 = 500\), que minimicen la longitud media de corrida fuera de control. Estas tablas se muestran en el Apéndice B.

Hunter (1989) ha estudiado el EWMA y ha sugerido escoger \(\lambda\) de tal modo que la ponderación dada a las observaciones actuales y anteriores coincida en la medida de lo posible con las ponderaciones dadas a estas observaciones en una carta Shewhart con las reglas de la Western Electric. Se llega así a un valor recomendado de \(\lambda=0,\!4\). Si \(L=3,\!054\), entonces la Tabla 5.10 indica que esta carta tendría el valor \({\textit{LMC}}_0 \simeq 500\) y, para detectar un corrimiento de una desviación estándar en la media del proceso, \({\textit{LMC}}_1 \simeq 14,\!3\).

En general, es práctica habitual tomar valores de \(\lambda\) en el intervalo \(0,\!05 ≤ \lambda ≤ 0,\!25\), dado que funcionan bien en la práctica, siendo \(\lambda = 0,05\), \(\lambda = 0,\!10\) y \(\lambda = 0,\!20\) elecciones populares. Una buena regla empírica es usar valores pequeños de \(\lambda\) para detectar los corrimientos más pequeños. Se ha encontrado también que \(L = 3\) (los límites tres sigma usuales) funciona razonablemente bien, en particular con el valor más grande de \(\lambda\), aunque cuando \(\lambda\) es pequeña (por ejemplo, \(\lambda ≤ 0,\!10\)), se consigue una ventaja al reducir la anchura de los límites de control utilizando un valor de \(L\) entre 2,6 y 2,8. Recuerde que en el ejemplo del diagrama EWMA se usó \(\lambda = 0,\!10\) y \(L = 2,\!7\). Se esperaría que esta elección de los parámetros resultara en una \(\textit{LMC}\) bajo control de \({\textit{LMC}}_0 = 500\) y en un \(\textit{LMC}\) para detectar un corrimiento de una desviación estándar en la media de \({\textit{LMC}}_1 = 1010,\!30\). Por tanto, este diseño es el equivalente aproximado de la CUSUM con \(h = 5\) y \(k = 1/2\).

Al igual que la CUSUM, el EWMA tiene un buen desempeño contra los corrimientos pequeños, pero no reacciona con tanta rapidez como el diagrama de Shewhart a los corrimientos grandes. Sin embargo, el EWMA con frecuencia es superior a la CUSUM para los corrimientos grandes, en particular si \(λ > 0,\!10\). Una buena manera de mejorar adicionalmente la sensibilidad del procedimiento de control a los corrimientos grandes sin sacrificar la habilidad para detectar corrimientos pequeños es combinar el diagrama de Shewhart con el EWMA. Estos procedimientos de control combinados Shewhart-EWMA son efectivos tanto contra los corrimientos pequeños como contra los corrimientos grandes. Al usar estos esquemas, D. C. Montgomery (2005), ha encontrado útil usar límites con anchura un poco mayor que lo usual en el diagrama de Shewhart (por ejemplo, 3,25 sigma, o incluso 3,5 sigma). También es posible graficar tanto \(X_i\) (o \(\bar{X}_i\)) como el estadístico EWMA \(Z_i\) en la misma gráfica de control junto con los límites de Shewhart y EWMA. Se produce así una carta para el procedimiento de control combinado en cuya interpretación los operadores pronto se vuelven expertos. Cuando las gráficas son generadas por computadora, pueden usarse diferentes colores o símbolos gráficos para los dos conjuntos de límites de control y estadísticos.

5.2.2 Diagrama de Control EWMA para Subgrupos Racionales

El diagrama de control EWMA suele usarse con mediciones individuales. Sin embargo, si se toman subgrupos racionales de tamaño \(n > 1\), entonces el estadístico EWMA dado en la ecuación 5.12 se redefine sustituyendo \(X_i\) con \(\overline{X}_i\), como se muestra a continuación.

Estadístico EWMA recursivo para subgrupos racionales

\[ \begin{equation} Z_{i} = \lambda \overline{X}_i + \left( 1 - \lambda \right) Z_{i - 1} \end{equation} \:. \tag{5.17}\]

Siguiendo un procedimiento análogo al realizado para obtener la media y la varianza del estadístico EWMA recursivo para observaciones individuales dado en la ecuación 5.12, se puede verificar que la media y la varianza del estadístico EWMA recursivo para subgrupos racionales están determinadas por:

Media y varianza de \(Z_i\) para subgrupos racionales

\[ \begin{align*} \mu_{_{Z_i}} & = \mu_{0} \\ \sigma_{_{\!Z_i}}^{2} & = \frac{\sigma^{2}}{n}\left( \frac{\lambda}{2 - \lambda} \right)\left[ 1 - \left( 1 - \lambda \right)^{2i} \right] \end{align*}\:. \tag{5.18}\]

De aquí que los límites de control para el diagrama EWMA para subgrupos racionales vienen dados por

Límites de control estándar del diagrama de control EWMA para subgrupos racionales

\[ \begin{align*} LCS & = \mu_{0} + L\sigma\sqrt{\left[ \frac{\lambda}{n\left( 2 - \lambda \right)} \right]\left[ 1 - \left( 1 - \lambda \right)^{2i} \right]} \\ LC & = \mu_{0} \\ LCI & = \mu_{0} - L\sigma\sqrt{\left[ \frac{\lambda}{n\left( 2 - \lambda \right)} \right]\left[ 1 - \left( 1 - \lambda \right)^{2i} \right]} \\ \end{align*} \:. \tag{5.19}\]

Note en la ecuación 5.19, si tomamos \(L = 3\) y \(\lambda = 1\) se observa que los límites dados en esta ecuación son iguales a los límites de control del diagrama \(\bar{x}\) mostrados en la ecuación 4.4. Lo que implica que, bajo estas condiciones, el diagrama \(\bar{x}\) es equivalente al diagrama EWMA.

Por otro lado, como ya se dijo anteriormente, el término \(1 - \left( 1 - \lambda \right)^{2i}\) en la ecuación 5.19 tiende a uno cuando \(i\) se hace grande, de manera que la varianza de \(Z_i\) se va incrementando y con ello la longitud de los límites de control del diagrama EWMA se incrementan en los primeros puntos hasta estabilizarse en:

\[ \begin{align*} LCS & = \mu_{0} + L\sigma\sqrt{\left[ \frac{\lambda}{n\left( 2 - \lambda \right)} \right]} \\ LC & = \mu_{0} \\ LCI & = \mu_{0} - L\sigma\sqrt{\left[ \frac{\lambda}{n\left( 2 - \lambda \right)} \right]} \\ \end{align*} \:. \tag{5.20}\]

5.2.3 Robustez del Diagrama de Control EWMA para la Normalidad

En la discusión sobre los diagramas de control para variables tipo Shewhart planteada en el Capítulo 4, se destacó que el diagrama para observaciones individuales es muy sensible a la falta de normalidad de la característica de la calidad debido a que provoca que la longitud media de corrida bajo control (\({\textit{LMC}}_0\)) sea notablemente menor que el valor esperado bajo el supuesto de normalidad. Este problema es abordado por Borror, Montgomery, y Runger (1999) comparando el desempeño de la \(\textit{LMC}\) de los diagramas Shewhart y EWMA para observaciones individuales para el caso de distribuciones no normales. Específicamente, usaron la distribución gamma para representar el caso de distribuciones sesgadas y la distribución \(t\) para representar distribuciones simétricas con colas más altas que la distribución normal.

En la Tabla 5.10 y Tabla 5.12 se presenta la \({\textit{LMC}}_0\) del diagrama Shewhart para mediciones individuales y varios diagramas de control EWMA para distribuciones no normales. Hay dos aspectos de la información contenida en estas tablas que llaman mucho la atención. Primero, incluso las distribuciones que se apartan moderadamente de la normalidad tienen el efecto de reducir de manera considerable la \({\textit{LMC}}_0\) del diagrama Shewhart para mediciones individuales. Esto, desde luego, produce un incremento radical de la tasa de falsa alarma. Segundo, un EWMA con \(\lambda = 0,05\) o \(\lambda = 0,10\) y un límite de control elegido apropiadamente tendrá un desempeño muy bueno para distribuciones tanto normales como no normales. Cuando \(\lambda = 0,05\) y \(L=2,492\), la \({\textit{LMC}}_0\) de la EWMA está aproximadamente dentro del 8% de la \({\textit{LMC}}_0\) determinada por la teoría normal de 370, excepto en los casos más extremos. Además, las propiedades de detección de corrimiento del EWMA son todas uniformemente superiores a las del diagrama Shewhart para mediciones individuales.

Código 
```{r}
#| label: tbl-lmc-ewma-gamma
#| tbl-cap: "${\\textit{LMC}}_0$ para los diagramas de control EWMA y de mediciones individuales para varias distribuciones gamma"
#| table-caption: left

comp_gamma <- tibble::tibble(
  dist = c(
    "<span style='font-size: 1em;'><strong>Normal</strong></span>",
    "<span style='font-size: 1em;'><strong>Gamma(4, 1)</strong></span>",
    "<span style='font-size: 1em;'><strong>Gamma(3, 1)</strong></span>",
    "<span style='font-size: 1em;'><strong>Gamma(2, 1)</strong></span>",
    "<span style='font-size: 1em;'><strong>Gamma(1, 1)</strong></span>",
    "<span style='font-size: 1em;'><strong>Gamma(0,5, 1)</strong></span>"
  ),
  EWMA1 = c(370.4, 372, 372, 372, 369, 357),
  EWMA2 = c(370.8, 341, 332, 315, 274, 229),
  EWMA3 = c(370.5, 259, 238, 208, 163, 131),
  Shewhart = c(370.4, 97, 85, 71, 55, 45)
)

comp_gamma |>
  gt::gt() |>
   cols_label(
    dist = html(
      "<p><em><strong>Distribuci&oacute;n</strong></em></p>"
    ),
    EWMA1 = html("<p><em><strong>L = 2,492</strong></em></p>"),
    EWMA2 = html("<p><em><strong>L = 2,703</strong></em></p>"),
    EWMA3 = html("<p><em><strong>L = 2,860</strong></em></p>"),
    Shewhart = html("<p><em><strong>L = 3</strong></em></p>"),
  ) |>
  tab_spanner(
    label = html(
      "<p><em><strong>&lambda; = 0,05</strong></em></p>"
    ),
    columns = c(EWMA1)
  ) |>
  tab_spanner(
    label =  html(
      "<p><em><strong>&lambda; = 0,10</strong></em></p>"
    ),
    columns = c(EWMA2)
  ) |>
  tab_spanner(
    label =  html(
      "<p><em><strong>&lambda; = 0,20</strong></em></p>"
    ),
    columns = c(EWMA3)
  ) |>
  tab_spanner(
    label =  html(
      "<p><em><strong>&lambda; = 1</strong></em></p>"
    ),
    columns = c(Shewhart)
  ) |>
  tab_spanner(
    label =  html(
      "<p><em><strong>EWMA</strong></em></p>"
    ),
    columns = c("EWMA1", "EWMA2", "EWMA3")
  ) |>
  tab_spanner(
    label =  html(
      "<p><em><strong>Shewhart</strong></em></p>"
    ),
    columns = c("Shewhart")
  ) |>
  fmt_number(decimals = 1, sep_mark = ".", dec_mark = ",") |>
  opt_align_table_header(align = "left") |>
  cols_align(
    align = "center",
    columns = everything()
  ) |> 
  fmt_markdown(
    columns = dist
  ) |> 
  cols_align(
    align = "center",
    columns = everything()
  ) |> 
  tab_style(
    style = list(
      cell_text(align = "center")
    ),
    locations = cells_body(columns = everything())
  )
```
Tabla 5.11: \({\textit{LMC}}_0\) para los diagramas de control EWMA y de mediciones individuales para varias distribuciones gamma

EWMA

Shewhart

Distribución

λ = 0,05

λ = 0,10

λ = 0,20

λ = 1

L = 2,492

L = 2,703

L = 2,860

L = 3
Normal 370,4 370,8 370,5 370,4
Gamma(4, 1) 372,0 341,0 259,0 97,0
Gamma(3, 1) 372,0 332,0 238,0 85,0
Gamma(2, 1) 372,0 315,0 208,0 71,0
Gamma(1, 1) 369,0 274,0 163,0 55,0
Gamma(0,5, 1) 357,0 229,0 131,0 45,0

Basados en la información proporcionada por la Tabla 5.10 y Tabla 5.12, se recomienda usar el diagrama de control EWMA adecuadamente diseñado como gráfico de control para mediciones individuales en una amplia gama de aplicaciones, especialmente en la monitorización de procesos de la fase II. Es casi un procedimiento perfectamente no paramétrico (libre de distribución).

Código 
```{r}
#| label: tbl-lmc-ewma-t
#| tbl-cap: "${\\textit{LMC}}_0$ para los diagramas de control EWMA y de mediciones individuales para varias distribuciones $t$"
#| table-caption: left

comp_t <- tibble::tibble(
  dist = c(
    "<span style='font-size: 1em;'><strong>Normal</strong></span>", 
    "<span style='font-size: 1em;'><em><strong>t<sub>50</sub></strong></em></span>", 
    "<span style='font-size: 1em;'><em><strong>t<sub>40</sub></strong></em></span>", 
    "<span style='font-size: 1em;'><em><strong>t<sub>30</sub></strong></em></span>",
    "<span style='font-size: 1em;'><em><strong>t<sub>20</sub></strong></em></span>", 
    "<span style='font-size: 1em;'><em><strong>t<sub>15</sub></strong></em></span>", 
    "<span style='font-size: 1em;'><em><strong>t<sub>10</sub></strong></em></span>", 
    "<span style='font-size: 1em;'><em><strong>t<sub>8</sub></strong></em></span>", 
    "<span style='font-size: 1em;'><em><strong>t<sub>6</sub></strong></em></span>", 
    "<span style='font-size: 1em;'><em><strong>t<sub>4</sub></strong></em></span>"
  ),
  EWMA1 = c(
    370.4, 369, 369, 368, 367, 365,
    361, 358, 351, 343
  ),
  EWMA2 = c(
    370.8, 365, 363, 361, 355, 349, 335, 324, 305, 274
  ),
  EWMA3 = c(
    370.5, 353, 348, 341, 325, 310, 280, 259, 229, 188
  ),
  Shewhart = c(
    370.4, 283, 266, 242, 204, 176, 137, 117, 96, 76
  )
)

# knitr::kable(
#   comp_t,
#   format = "markdown",
#   digits = 3,
#   format.args = list(decimal.mark = ",", big.mark = "."),
#   col.names = c(
#     "Distribución",
#     "$\\textbf{EWMA}\\left(\\lambda = 0,05; L = 2, 492 \\right)$",
#     "$\\textbf{EWMA}\\left(\\lambda = 0,1; L = 2, 703 \\right)$",
#     "$\\textbf{EWMA}\\left(\\lambda = 0,2; L = 2, 86 \\right)$",
#     "$\\textbf{Shewhart}\\left(\\lambda = 12; L = 3 \\right)$"
#   ),
#   align = "lrrrr",
#   escape = FALSE
# )

comp_t |>
  gt::gt() |>
   cols_label(
    dist = html(
      "<p><em><strong>Distribuci&oacute;n</strong></em></p>"
    ),
    EWMA1 = html("<p><em><strong>L = 2,492</strong></em></p>"),
    EWMA2 = html("<p><em><strong>L = 2,703</strong></em></p>"),
    EWMA3 = html("<p><em><strong>L = 2,860</strong></em></p>"),
    Shewhart = html("<p><em><strong>L = 3</strong></em></p>"),
  ) |>
  tab_spanner(
    label = html(
      "<p><em><strong>&lambda; = 0,05</strong></em></p>"
    ),
    columns = c(EWMA1)
  ) |>
  tab_spanner(
    label =  html(
      "<p><em><strong>&lambda; = 0,10</strong></em></p>"
    ),
    columns = c(EWMA2)
  ) |>
  tab_spanner(
    label =  html(
      "<p><em><strong>&lambda; = 0,20</strong></em></p>"
    ),
    columns = c(EWMA3)
  ) |>
  tab_spanner(
    label =  html(
      "<p><em><strong>&lambda; = 1</strong></em></p>"
    ),
    columns = c(Shewhart)
  ) |>
  tab_spanner(
    label =  html(
      "<p><em><strong>EWMA</strong></em></p>"
    ),
    columns = c("EWMA1", "EWMA2", "EWMA3")
  ) |>
  tab_spanner(
    label =  html(
      "<p><em><strong>Shewhart</strong></em></p>"
    ),
    columns = c("Shewhart")
  ) |>
  fmt_number(decimals = 1, sep_mark = ".", dec_mark = ",") |>
  opt_align_table_header(align = "left") |>
  cols_align(
    align = "center",
    columns = everything()
  ) |> 
  cols_label(
    dist = html(
      "<p><em><strong>Distribuci&oacute;n</strong></em></p>"
    )
  ) |> 
  fmt_markdown(
    columns = dist
  ) |> 
  cols_align(
    align = "center",
    columns = everything()
  ) |> 
  tab_style(
    style = list(
      cell_text(align = "center")
    ),
    locations = cells_body(columns = everything())
  )
```
Tabla 5.12: \({\textit{LMC}}_0\) para los diagramas de control EWMA y de mediciones individuales para varias distribuciones \(t\)

EWMA

Shewhart

Distribución

λ = 0,05

λ = 0,10

λ = 0,20

λ = 1

L = 2,492

L = 2,703

L = 2,860

L = 3
Normal 370,4 370,8 370,5 370,4
t50 369,0 365,0 353,0 283,0
t40 369,0 363,0 348,0 266,0
t30 368,0 361,0 341,0 242,0
t20 367,0 355,0 325,0 204,0
t15 365,0 349,0 310,0 176,0
t10 361,0 335,0 280,0 137,0
t8 358,0 324,0 259,0 117,0
t6 351,0 305,0 229,0 96,0
t4 343,0 274,0 188,0 76,0

5.2.4 Diagrama de Control EWMA con la Cararística de Respuesta Inicial Rápida

Es posible incorporar al diagrama de control EWMA la característica de respuesta inicial rápida (FIR) o ventaja anticipada. Al igual que en el diagrama de control CUSUM, la ventaja de este procedimiento es detectar más rápidamente un proceso que está fuera del objetivo desde el inicio.

En este libro se aconsejan dos enfoques. El primero, desarrollado por Rhoads, Montgomery, y Mastrangelo (1996) los cuales establecen dos diagramas EWMA unilaterales y los inician a la mitad entre el objetivo y el límite de control. Se supone que ambos diagramas unilaterales tienen límites que varían en el tiempo.

El segundo enfoque fue propuesto por Steiner (1999) y este utiliza una sola carta de control, pero estrecha aún más la amplitud de los límites que varían con el tiempo para los primeros puntos de la muestra. Emplea un ajuste exponencialmente decreciente para estrechar aún más los límites, de modo que los límites de control se ubiquen a una distancia

\[ \pm \, L\sigma \left\{\left [ 1 - \left ( 1 -f \right )^{1 + a\left ( t -1 \right )} \right ] \sqrt{\frac{\lambda }{2 -\lambda }}\left [ 1 - \left ( 1 - \lambda \right )^{2t} \right ]\right\} \]

alrededor del objetivo.

Las constantes \(f\) y \(a\) deben determinarse. Steiner sugiere elegir \(a\) de manera que la FIR tenga poco efecto después de unas 20 observaciones. Esto lleva a elegir \(a = \left[ −2/\log\left( 1 − f \right) − 1 \right]/19\). Por ejemplo, si \(f=0,5\), entonces \(a = 0,3\). La elección de \(f=0,5\) es atractiva porque simula la ventaja anticipada del 50% que a menudo se usa en el diagrama de control CUSUM, discutida en Sección 5.1.6.

Ambos métodos son efectivos para disminuir la \(\textit{LMC}\) y detectar rápidamente un proceso que está fuera del objetivo al arranque. En la práctica, el procedimiento de Steiner es más sencillo de aplicar.

5.2.5 Diagrama de Control EWMA para Monitorear la Variabilidad del Proceso

MacGregor y Harris (1993) presentan el siguiente estadístico basado en el error cuadrático medio exponencialmente ponderado (EWMS, por su siglas en inglés):

Estadístico EWMS Recursivo

\[ \begin{align*} S_{i}^2 = \lambda\left(X_i - \mu\right)^2 + \left(1 - \lambda \right)S_{i-1}^2 \end{align*} \:. \tag{5.21}\]

Estos autores demuestran que bajo el supuesto de que los \(X_i\) son variables aleatorias normalmente distribuidas con media \(\mu\) y desviación estándar \(\sigma\), entonces \(E\left(S_{i}^2\right) = \sigma^2\) (para \(i\) grande) y si las observaciones \(X_i\) son independientes \(S_{i}^2/\sigma^2\) tiene una distribución aproximada chi-cuadrado con \(v=\left( 2 - \lambda \right)/\lambda\) grados de libertad. En consecuencia, si \(\sigma_{_0}\) representa el valor bajo control u objetivo de la desviación estándar del proceso, \(\sqrt{S_{i}^2}\) podría graficarse en un diagrama de control de la raíz cuadrada media exponencialmente ponderada o EWRMS (por sus siglas en inglés) con límites de control:

Límites de control del Diagrama EWMA para la Variabilidad

\[ \begin{align*} LCS &= \sigma_{_0}\sqrt{\frac{\chi_{v,\, 1 - \left( \alpha / 2 \right)}^2}{v}} \\ LC & = \sigma_{_0} \\ LCS &= \sigma_{_0}\sqrt{\frac{\chi_{v,\, \alpha / 2 }^2}{v}} \\ \end{align*} \:. \tag{5.22}\]

MacGregor y Harris (1993) señalan que la estadística EWMS puede ser sensible a cambios tanto en la media como en la desviación estándar del proceso. Por lo tanto, sugieren reemplazar \(\mu\) en la ecuación 5.21 con una estimación en cada punto en el tiempo. Una estimación lógica de \(\mu\) resulta ser el EWMA \(Z_i\) ordinaria. Estos autores derivan límites de control para la varianza móvil exponencialmente ponderada exponencialmente resultante (EWMV, por sus siglas en inglés).

5.2.6 Diagrama de Control EWMA para Datos de Poisson

Así como la CUSUM puede usarse como base de una carta de control efectiva para conteo de Poisson, lo mismo puede hacerse con un EWMA diseñado apropiadamente. Borror, Champ, y Rigdon (1998) describen el procedimiento, indican la forma de diseñar el diagrama de control y ofrecen un ejemplo. Si \(X_i\) es un conteo, entonces la recursión básica del EWMA se mantiene sin cambio:

\[ \begin{equation} Z_{i} = \lambda X_i + \left( 1 - \lambda \right) Z_{i - 1} \end{equation} \tag{5.23}\]

con \(z_0 = \mu_0\) la tasa de conteo bajo control u objetivo. Los parámetros de la carta de control son los siguientes:

Límites de control del Diagrama EWMA para Datos de Poisson

\[ \begin{align*} LCS & = \mu_{0} + A_S\sqrt{\left( \frac{\lambda\mu_0}{2 - \lambda} \right)\left[ 1 - \left( 1 - \lambda \right)^{2i} \right]} \\ LC & = \mu_{0} \\ LCI & = \mu_{0} - A_I\sqrt{\left( \frac{\lambda \mu_0}{2 - \lambda} \right)\left[ 1 - \left( 1 - \lambda \right)^{2i} \right]} \\ \end{align*} \:. \tag{5.24}\]

donde \(A_S\) y \(A_I\) son los factores de los límites de control superior e inferior, respectivamente. En muchas aplicaciones se escoge \(A_S = A_I = A\). Borror, Champ, y Rigdon (1998) presentan gráficas del desempeño de la \(\textit{LMC}\) del diagrama de control EWMA de Poisson como una función de \(\lambda\) y \(A\) y para varias tasas de conteo \(\mu_0\) bajo control u objetivo. Una vez que se determina \(\mu_0\) y que se especifica un valor para \(\lambda\), estos diagramas pueden usarse para seleccionar el valor de \(A\) que produce la \({\textit{LMC}}_0\) deseada. Los autores indican asimismo que la habilidad de este diagrama de control para detectar causas asignables es considerablemente mejor que la del diagrama \(c\) de Shewhart.

5.3 Diagrama de Control del Promedio Móvil

Tanto en el diagrama de control CUSUM como el EWMA los estadísticos graficados son ponderados en el tiempo. El diagrama EWMA utiliza un promedio ponderado como el estadístico del diagrama. Ocasionalmente, puede ser de interés otro tipo de diagrama de control ponderado con el tiempo basada en un promedio móvil no ponderado simple.

Supóngase que se tiene un proceso de producción en el cual se está monitoreando una característica de la calidad \(X\), y que \(X_1,X_2, \dotsc, X_m\) denotan cualquier muestra de tamaño \(m\) tomada del proceso. Entonces, el promedio o media móvil de extensión \(w\) en el tiempo \(i\) se define como

Estadístico media móvil de extensión \(w\)

\[ \begin{equation} M_i = \frac{X_{i - w + 1} + X_{i - w + 2}+ \cdots + X_{i}}{w} = \frac{\sum_{j = i - w +1}^iX_j}{w} \end{equation} \:. \tag{5.25}\]

Es decir, en el periodo de tiempo \(i\), se sustituye la observación más antigua del conjunto del promedio móvil y se agrega al conjunto la más reciente.

Note de la ecuación anterior (ecuación 5.25) que el rango móvil de extensión \(w\) solo puede ser calculado a partir de los periodos \(i \ge w\). Por lo que, para los periodos \(i < w\) se usa el estadístico media móvil de extensión \(i\) . Es decir,

\[ \begin{equation} M_i = \frac{\sum_{j = 1}^{i}X_j}{i} \end{equation} \:. \tag{5.26}\]

De manera general, se define el estadístico promedio o media móvil de la siguiente manera,

Estadístico media móvil

\[ M_i = \begin{cases} \frac{\sum_{j = 1}^{i}X_j}{i}, & \text{ si }\:i < w \\ \frac{\sum_{j = i - w +1}^iX_j}{w}, & \text{ si }\: i \ge w \end{cases} \tag{5.27}\]

Si las observaciones \(X_i\) son variables aleatorias independientes con media \(\mu_0\) y varianza \(\sigma^2\), entonces la media de \(M_i\) viene dado por:

  • Media de \(M_i\), para \(i < 5\) \[ \begin{align*} \mu_{_{M_i}} &= E\left( M_i \right) = E\left( \frac{\sum_{j = 1}^{i}X_j}{i} \right) = \frac{1}{i}E\left( \sum_{j = 1}^{i}X_j\right) \\ &= \frac{1}{i}\sum_{j = 1}^{i}E\left( X_j \right)= \frac{1}{i}\sum_{j = 1}^{i}\mu_0 = \frac{1}{i}i\mu_0 = \mu_0\\ \end{align*} \]
  • Media de \(M_i\), para \(i \ge 5\) \[ \begin{align*} \mu_{_{M_i}} &= E\left( M_i \right) = E\left( \frac{\sum_{j = i - w +1}^{i}X_j}{w} \right) = \frac{1}{w}E\left( \sum_{j = i - w +1}^{i}X_j\right) \\ &= \frac{1}{w}\sum_{j = i - w + 1}^{i}E\left( X_j \right)= \frac{1}{w}\sum_{j = i - w + 1}^{i}\mu_0 = \frac{1}{w}w\mu_0 = \mu_0\\ \end{align*} \]

Mientras que la varianza del promedio móvil \(M_i\) es:

  • Varianza de la media móvil para \(i < w\) \[ \begin{align*} \sigma_{_{\!M_i}}^{2} &= V\left( M_i \right)= V\left( \frac{\sum_{j = i - w +1}^{i}X_j}{i} \right) = \frac{1}{i^2}V\left( \sum_{j = i - w +1}^{i}X_j\right) \\ &= \frac{1}{i^2}\sum_{j = i - w + 1}^{i}V\left( X_j \right)= \frac{1}{i^2}\sum_{j = i - w + 1}^{i}\sigma^2 = \frac{1}{i^2}i\sigma^2 = \frac{\sigma^2}{i}\\ \end{align*} \]
  • Varianza de la media móvil para \(i \ge w\) \[ \begin{align*} \sigma_{_{\!M_i}}^{2} &= V\left( M_i \right)= V\left( \frac{\sum_{j = i - w +1}^{i}X_j}{w} \right) = \frac{1}{w^2}V\left( \sum_{j = i - w +1}^{i}X_j\right) \\ &= \frac{1}{w^2}\sum_{j = i - w + 1}^{i}V\left( X_j \right)= \frac{1}{w^2}\sum_{j = i - w + 1}^{i}\sigma^2 = \frac{1}{w^2}w\sigma^2 = \frac{\sigma^2}{w}\\ \end{align*} \]

En resumen, la media y la varianza del estadístico \(M_i\) vienen dadas por

Media y varianza de \(M_i\)

\[ \begin{align*} \mu_{_{W_i}} & = \mu_{0} \\ \sigma_{_{\!W_i}}^{2} & = \begin{cases} \frac{\sigma^2}{i}, & \text{ si }\:i < w \\ \frac{\sigma^2}{w}, & \text{ si }\: i \ge w \end{cases} \end{align*} \:. \tag{5.28}\]

Por lo tanto, si \(μ_0\) denota el valor objetivo del proceso, entonces la línea central y los límites de control tres sigma para el diagrama del promedio móvil \(M_i\) vienen dado por

Límites de control del diagrama Media Móvil

\[ \begin{align*} LIC & = \begin{cases} \mu_0 - 3\frac{\sigma}{\sqrt{i}}, & \text{ si }\:i < w \\ \mu_0 - 3\frac{\sigma}{\sqrt{w}}, & \text{ si }\: i \ge w \end{cases} \\ LSC & = \begin{cases} \mu_0 + 3\frac{\sigma}{\sqrt{i}}, & \text{ si }\:i < w \\ \mu_0 + 3\frac{\sigma}{\sqrt{w}}, & \text{ si }\: i \ge w \end{cases} \\ \end{align*} \:. \tag{5.29}\]

Nota

Como se puede ver de la ecuación anterior (ecuación 5.29), los límites de control del diagrama promedio móvil para los periodos \(i < w\) son variables, mientras que para los periodos \(i \ge w\) son constantes (ver figura Figura 5.17)

El procedimiento de control consistiría en calcular el nuevo promedio móvil \(M_i\) cada vez que se disponga de una nueva observación \(X_i\), graficar \(M_i\) en un diagrama de control con línea central y límites de control superior e inferior dados por la ecuación ecuación 5.29, y concluir que el proceso está fuera de control si \(M_i\) excede los límites de control. En general, la magnitud del corrimiento de interés y \(w\) están relacionados de manera inversa; se contaría con una protección más efectiva para los corrimientos más pequeños utilizando promedios móviles de extensión más grande, a costa de respuesta rápida a los corrimientos grandes.

Ejemplo 5.5 (Ejemplo de Aplicación del Diagrama de Control del Promedios Móviles para Observaciones Individuales)
Establezca un diagrama de control para los datos de la Tabla 5.1, utilizando \(w = 5\). Recuerde que para estos datos el valor objetivo de la media es \(\mu_0 = 100\) y que la desviación estándar es \(\sigma = 5\).

En la tabla Tabla 5.13 se muestran las observaciones \(x_i\), los cálculos de los \(M_i\) y los límites de control, para los periodos \(1 \le i \ge 30\).

Código 
```{r}
#| label: tbl-datos-promedio-movil
#| tbl-cap: "Diagrama de control del promedio móvil para los datos de la @tbl-datos-cusum"
#| table-caption: left
#| class: output

# Función para calcular los promedios móviles
media_movil <- function(data, w) {
  m_i <- vector(mode = "double", length =  length(data))
  for (i in 1:length(data)) {
    ifelse(
      i < w,
      m_i[i] <- mean(data[1:i]),
      m_i[i] <- mean(data[(i - w + 1):i])
    )
  }
  return(m_i)
}

# Función para calcular los límites de control del diagrama de control promedio móvil

lc_dcmm <- function(data, w, media, ds) {
  LSC <- vector(mode = "double", length =  length(data))
  LIC <- vector(mode = "double", length =  length(data))
  for (i in 1:length(data)) {
    ifelse(
      i < w,
      {
        LSC[i] <- media + 3 * ds / sqrt(i)
        LIC[i] <- media - 3 * ds / sqrt(i)
      },
      {
        LSC[i] <- media + 3 * ds / sqrt(w)
        LIC[i] <- media - 3 * ds / sqrt(w)
      }
    )
  }
  return(list(LIC, LSC))
}

data_mm <- data_cusum |>
  dplyr::select(muestra, x_i) |>
  dplyr::mutate(
    m_i = media_movil(data = x_i, w = 5),
    LIC = lc_dcmm(
      data = data_cusum$x_i, w = 5, media = 100, ds = 5
    )[[1]],
    LSC = lc_dcmm(
      data = data_cusum$x_i, w = 5, media = 100, ds = 5
    )[[2]]
  )

knitr::kable(
  data_mm,
  format = "markdown",
  digits = 4,
  format.args = list(decimal.mark = ",", big.mark = "."),
  col.names = c(
    "$Muestra \\: (i)$", "$x_i$", "$M_i$", "$LIC$", "$LSC$"
  ),
  align = "ccccc",
  escape = FALSE
  )
```
Tabla 5.13: Diagrama de control del promedio móvil para los datos de la Tabla 5.1
\(Muestra \: (i)\) \(x_i\) \(M_i\) \(LIC\) \(LSC\)
1 96,9663 96,9663 85,0000 115,0000
2 97,2494 97,1078 89,3934 110,6066
3 103,6353 99,2837 91,3397 108,6603
4 96,7796 98,6577 92,5000 107,5000
5 110,3248 100,9911 93,2918 106,7082
6 97,7647 101,1508 93,2918 106,7082
7 97,4743 101,1958 93,2918 106,7082
8 100,8277 100,6342 93,2918 106,7082
9 91,8081 99,6399 93,2918 106,7082
10 107,2495 99,0249 93,2918 106,7082
11 94,6188 98,3957 93,2918 106,7082
12 105,6507 100,0310 93,2918 106,7082
13 94,9607 98,8576 93,2918 106,7082
14 107,0322 101,9024 93,2918 106,7082
15 102,3922 100,9309 93,2918 106,7082
16 100,0102 102,0092 93,2918 106,7082
17 96,8418 100,2474 93,2918 106,7082
18 97,6166 100,7786 93,2918 106,7082
19 97,9844 98,9690 93,2918 106,7082
20 99,4052 98,3716 93,2918 106,7082
21 110,7501 100,5196 93,2918 106,7082
22 102,4706 101,6454 93,2918 106,7082
23 108,5581 103,8337 93,2918 106,7082
24 113,1497 106,8667 93,2918 106,7082
25 112,3129 109,4483 93,2918 106,7082
26 107,0480 108,7078 93,2918 106,7082
27 110,5134 110,3164 93,2918 106,7082
28 95,7660 107,7580 93,2918 106,7082
29 107,5611 106,6403 93,2918 106,7082
30 98,1360 103,8049 93,2918 106,7082


El estadístico graficado en el diagrama de control del promedio móvil es

\[ \begin{equation} M_i = \frac{X_{i - 4} + X_{i - 3} + \cdots + X_{i }}{5} \end{equation}\:, \]

para los periodos \(i \ge 5\). Para los periodos de tiempo \(i < 5\), se grafica el promedio de las observaciones para los periodos \(1, 2, \dotsc, i\), es decir,

\[ \begin{equation} M_i = \frac{\sum_{j = 1}^{i} X_{i}}{i} \end{equation}\:. \]

A continuación, se muestra el cálculo de los primeras 7 medias móviles:

  • Media móvil para el periodo 1 \[ m_1 = \frac{x_1}{1} = \frac{96,9663}{1} = 96,9663 \]
  • Media móvil para el periodo 2 \[ m_2 = \frac{x_1 + x_2}{2} = \frac{96,9663 + 97,2494}{2} = 97,1078 \]
  • Media móvil para el periodo 3 \[ m_3 = \frac{x_1 + x_2 + x_3}{3} = \frac{96,9663 + 97,2494+ 103,6353}{3} = 99,2837 \]
  • Media móvil para el periodo 4 \[ \begin{align*} m_4 &= \frac{x_1 + x_2 + x_3 + x_4}{4} = \frac{96,9663 + 97,2494 + 103,6353+ 96,7796}{4}\\ &= 98,6577 \end{align*} \]
  • Media móvil para el periodo 5 \[ \begin{align*} m_5 &= \frac{x_1 + x_2 + x_3 + x_4 + X_5}{4}\\ &= \frac{96,9663 + 97,2494 + 103,6353+ 96,7796+ 110,3248}{5}\\ &= 100,9911 \end{align*} \]
  • Media móvil para el periodo 6 \[ \begin{align*} m_6 &= \frac{x_2 + x_3 + x_4 + x_5 + X_6}{4}\\ &= \frac{97,2494 + 103,6353 + 96,7796+ 110,3248+ 97,7647}{5}\\ &= 101,1508 \end{align*} \]
  • Media móvil para el periodo 7 \[ \begin{align*} m_7 &= \frac{x_3 + x_4 + x_5 + x_6 + X_7}{4}\\ &= \frac{103,6353 + 96,7796 + 110,3248+ 97,7647+ 97,4743}{5}\\ &= 101,1958 \end{align*} \]

Los límites de control para los primeros 5 periodos se muestran a continuación:

  • Periodo 1 \[ \begin{align*} {LIC}_1 &= {\mu}_0 - 3\frac{\sigma}{\sqrt{1}} = 100 - 3\frac{5}{\sqrt{1}} = 85\\ {LSC}_1 &= {\mu}_0 + 3\frac{\sigma}{\sqrt{1}} = 100 - 3\frac{5}{\sqrt{1}} = 115 \end{align*} \]
  • Periodo 2 \[ \begin{align*} {LIC}_2 &= {\mu}_0 - 3\frac{\sigma}{\sqrt{2}} = 100 - 3\frac{5}{\sqrt{2}} = 89,3934\\ {LSC}_2 &= {\mu}_0 + 3\frac{\sigma}{\sqrt{2}} = 100 + 3\frac{5}{\sqrt{2}} = 110,6066 \end{align*} \]
  • Periodo 3 \[ \begin{align*} {LIC}_3 &= {\mu}_0 - 3\frac{\sigma}{\sqrt{3}} = 100 - 3\frac{5}{\sqrt{3}} = 91,3397\\ {LSC}_3 &= {\mu}_0 + 3\frac{\sigma}{\sqrt{3}} = 100 + 3\frac{5}{\sqrt{3}} = 108,6603 \end{align*} \]
  • Periodo 4 \[ \begin{align*} {LIC}_4 &= {\mu}_0 - 3\frac{\sigma}{\sqrt{4}} = 100 - 3\frac{5}{\sqrt{4}} = 92,5\\ {LSC}_4 &= {\mu}_0 + 3\frac{\sigma}{\sqrt{4}} = 100 + 3\frac{5}{\sqrt{4}} = 107,5 \end{align*} \]
  • Periodo 5 \[ \begin{align*} {LIC}_5 &= {\mu}_0 - 3\frac{\sigma}{\sqrt{5}} = 100 - 3\frac{5}{\sqrt{5}} = 93,2918\\ {LSC}_5 &= {\mu}_0 + 3\frac{\sigma}{\sqrt{5}} = 100 + 3\frac{5}{\sqrt{5}} = 106,7082 \end{align*} \]

A continuación se muestra el diagrama de control promedio móvil generado a través de los paquetes ggplot2, highcharter y plotly.

Código 
```{r}
#| label: fig-dcmm-ggplot2
#| fig-cap: "Diagrama Promedio móvil para los datos de @tbl-datos-cusum con el paquete `ggplot2`"

q <-  qcc(
  data_mm$m_i[1:20], sizes = 5, type = "xbar", center = 100,
  std.dev = 5, newdata = data_mm$m_i[21:30], newsizes = 5,
  plot = FALSE
)

violacion <- union(
  q$violations$beyond.limits,
  q$violations$violating.runs
) |> 
  sort()

# summary(q)

violacion <- ifelse(
  data_mm$muestra %in% violacion, "Si", "No"
)

ggplot(
  data = data_mm, aes(x = muestra, y = m_i), show.legend = FALSE
) +
  geom_step(
    mapping = aes(x = muestra, y = LIC),
    direction = "hv", color = "red", lty = "dashed"
  ) +
  geom_step(
    mapping = aes(x = muestra, y = 100),
    direction = "hv", color = "blue"
  ) +
  geom_step(
    mapping = aes(x = muestra, y = LSC),
    direction = "hv", color = "red", lty = "dashed"
  ) +
  geom_point(
    mapping = aes(
      x = muestra, y = m_i, color = factor(violacion), 
      shape = factor(violacion)
    ),
    size = 2
  ) +
  scale_color_manual(
    values = c("No" = "black", "Si" = "red")
  ) +
  scale_shape_manual(values = c(16, 8)) + 
  geom_line(
    mapping = aes(x = muestra, y = m_i), color = "black"
  ) +
  pammtools::geom_stepribbon(
    aes(ymin = LIC, ymax = 100), fill = "blue", alpha = 0.2
  ) +
  pammtools::geom_stepribbon(
    aes(ymin = 100, ymax = LSC), fill = "green", 
    alpha = 0.2
  ) +
  geom_vline(
    xintercept = (20 + 21) / 2, linetype = 2, 
    color = "blue"
  ) +
  annotate(
    "text",
    x = nrow(data_mm) + 1, y = last(data_mm$LIC), 
    label = "LCI", color = "red"
  ) +
  annotate(
    "text",
    x = nrow(data_mm) + 1, y = last(data_mm$LSC), 
    label = "LCS", color = "red"
  ) +
  annotate(
    "text",
    x = nrow(data_mm) + 1, y = 100, label = "LC",
    color = "blue"
  ) +
  labs(
    x = latex2exp::TeX(
      "$Periodo \\, \\left(i\\right)$", 
      bold = TRUE, italic = TRUE
    ),
    y = latex2exp::TeX(
      "$\\textbf{Promedio \\, móvil}\\,\\left( M_i \\right)$",
      bold = TRUE, italic = TRUE
    ),
    shape = "Fuera de control", color = "Fuera de control"
  ) +
  theme_bw() +
  theme(
    legend.position = "bottom",
    legend.text = element_text(face = "bold"),
    legend.title.position = "top"
  ) 
```
Figura 5.17: Diagrama Promedio móvil para los datos de Tabla 5.1 con el paquete ggplot2
Código 
```{r}
#| label: fig-dcmm-highcharter
#| fig-cap: "Diagrama de contol promedio móvil con el paquete `highcharter`"


color <- ifelse(violacion == "Si", "#FF0000", "#000000")
symbol <- ifelse(violacion == "Si", "asterisk-open", "circle")

highchart() |>
  hc_add_series(
    data = data_mm, 
    hcaes(x = muestra, low = LSC, high = 100),
    step = "hv", 
    marker = list(enabled = FALSE, visible = FALSE),
    type = "arearange", fillOpacity = 0.2, 
    name = "", color = "green",
    showInLegend = FALSE, 
    tooltip = list(pointFormat = "{NULL}")
  ) |>
  hc_add_series(
    data = data_mm, 
    hcaes(x = muestra, low = LIC, high = 100),
    step = "hv", 
    marker = list(enabled = FALSE, visible = FALSE),
    type = "arearange", fillOpacity = 0.2, 
    name = "", color = "blue",
    showInLegend = FALSE, 
    tooltip = list(pointFormat = "{NULL}")
  ) |>
  hc_add_series(
    data_mm, type = "line", 
    hcaes(x = muestra, y = LSC), useHTML = TRUE,
    color = "red", dashStyle = "Dash", 
    name = "<b><i>LSC</i></b>",
    step = "hv", marker = list(enabled = FALSE),
    tooltip = list(valueDecimals = 4)
  ) |>
  hc_add_series(
    data_mm,
    type = "line", hcaes(x = muestra, y = 100), 
    useHTML = TRUE,
    color = "blue", name = "<b><i>LC</i></b>",
    step = "hv", marker = list(enabled = FALSE),
    tooltip = list(valueDecimals = 0)
  ) |>
  hc_add_series(
    data_mm, type = "line", hcaes(x = muestra, y = LIC), 
    useHTML = TRUE, color = "red", dashStyle = "Dash", 
    name = "<b><i>LIC</i></b>",
    step = "hv", marker = list(enabled = FALSE),
    tooltip = list(valueDecimals = 4)
  ) |>
  hc_add_series(
    data_mm, type = "line", color = "black",
    hcaes(
      x = muestra, y = m_i, color = color
    ),
    marker = list(radius = 3, color = color),
    useHTML = TRUE, 
    name = "<b><i>M<sub>i</sub><b><i>",
    dashStyle = "Dash",
    tooltip = list(valueDecimals = 4)
  ) |>
  hc_xAxis(
    title = list(
      useHTML = TRUE,
      text = "<b><i>Periodo (i)</i></b><br>"
    ),
    plotLines = list(
      list(
        color = "blue", value = 20.5, 
        width = 2, dashStyle = "Dash"
      )
    )
  ) |>
  hc_yAxis(
    title = list(
      useHTML = TRUE,
      text = "<b><i>Promedio móvil",
      startOnTick = FALSE
    )
  ) |>
  hc_tooltip(crosshairs = TRUE, shared = TRUE) |> 
  hc_plotOptions( series = list(animation = TRUE))
```
Figura 5.18: Diagrama de contol promedio móvil con el paquete highcharter
Código 
```{r}
#| label: fig-dcmm-plotly
#| fig-cap: "Diagrama promedio móvil con el paquete `plotly`"

vline <- function(x = 0, color = "blue") {
  list(
    type = "line",
    y0 = 0,
    y1 = 1,
    yref = "paper",
    x0 = x,
    x1 = x,
    line = list(color = color, dash = "dash")
  )
}

plot_ly(data_mm, x = ~muestra) |>
  add_lines(
    y = ~LSC, name = "<b><i>LSC</i></b>", 
    line = list(
      shape = "hvh", color = "red", dash = "dash"
    ),
    hovertemplate = paste(
      "<b><i>i = %{x}</i></b>",
      "<br><b><i>LSC<b><i> = %{y:.4f}<br><extra></extra>"
    )
  ) |>
  add_lines(
    y = 100, name = "<b><i>LC</i></b>", fill = "tonexty",
    fillcolor = "green",
    hovertemplate = paste(
      "<b><i>i = %{x}</i></b>",
      "<br><b><i>LC<b><i> = %{y:.0f}<br><extra></extra>"
    ),
    line = list(
      shape = "linear", color = "blue"
    )
  ) |>
  add_lines(
    y = ~LIC, name = "<b><i>LIC</i></b>", fill = "tonexty",
    fillcolor = "blue",
    hovertemplate = paste(
      "<b><i>i = %{x}</i></b>",
      "<br><b><i>LCI<b><i> = %{y:.4f}<br><extra></extra>"
    ),
    line = list(
      shape = "hvh", color = "red", dash = "dash"
    )
  ) |>
  add_trace(
    y = ~m_i,
    name = "<b><i>M<sub>i</sub><b><i>", type = "scatter",
    mode = "lines+markers",
    marker = list(color = color, symbol = symbol),
    hovertemplate = paste(
      "<b><i>i = %{x}</i></b>",
      "<br><b><i>M<sub>i</sub><b><i> = %{y:.4f}<br><extra></extra>"
    ),
    line = list(color = "black", dash = "dash")
  ) |>
  layout(
    legend = list(
      x = 0.5, y = -0.23, xref = "paper", yref = "paper",
      orientation = "h", xanchor = "center", showarrow = FALSE
    ),
    shapes = list(vline(20.5)),
    xaxis = list(
      title = "<b><i>Periodo</i></b><br>",
      zeroline = F
    ),
    yaxis = list(
      title = "<b><i>Promedio móvil (M<sub>i</sub>)</i></b><br>",
      zeroline = FALSE
    )
  ) |>
  plotly::config(locale = "es", mathjax = "cdn") |>
  fillOpacity(alpha = 0.3)
```
Figura 5.19: Diagrama promedio móvil con el paquete plotly

5.4 Diagrama de Control Multivariante

En los núcleos temáticos anteriores se planteó el problema del monitoreo y del control del proceso en la perspectiva básica de una sola variable; es decir, se ha supuesto que hay una sola variable de salida o característica de la calidad de interés. En la práctica, sin embargo, muchos, de no ser que la mayoría, de los escenarios de monitoreo y control incluyen varias variables relacionadas. Aun cuando una posible solución es la aplicación de cartas de control de una sola variable a cada variable individual, se verá que este enfoque es ineficiente y puede llevar a conclusiones equivocadas. En tal sentido se requieren diagramas de control multivariados que consideren a las variables en conjunto.

En esta sección se presentan diagramas \(T^2\) de Hotelling, tanto para observaciones individuales como para subgrupos, para monitorear de manera simultánea las medias de varias variables. Así como también, diagramas de control para monitorear la variabilidad en el caso de variables múltiples.

5.5 El Problema del Control de Calidad con Variables Múltiples

Hay muchas situaciones en las que es necesario el monitoreo o control simultáneo de dos o más características de la calidad relacionadas. Por ejemplo, suponer que un rodamiento tiene tanto un diámetro interior (\(x_1\)) como un diámetro exterior (\(x_2\)) que determinan en conjunto la utilidad de la pieza. Suponer que \(x_1\) y \(x_2\) tienen distribuciones normales independientes. Puesto que ambas características de la calidad son mediciones, podrían monitorearse aplicando la carta \(\bar{x}\) ordinaria a cada característica de la calidad, como se ilustra en la figura Figura 5.20. El proceso se considera bajo control si las medias muestrales \(\bar{x}_1\) y \(\bar{x}_2\) se localizan dentro de sus respectivos límites de control. Esto es equivalente a que el par de medias (\(\bar{x}_1\), \(\bar{x}_2\)) se localicen dentro de la región de control conjunta de la figura Figura 5.21.

Llevar a cabo el monitoreo independiente de estas dos características de la calidad puede ser muy engañoso. Por ejemplo, obsérvese de la figura Figura 5.21 que una de las observaciones parece ser un poco inusual con respecto a las demás. Este punto estaría dentro de los límites de control en las dos cartas \(\bar{x}\) de una variable para \(\bar{x}_1\) y \(\bar{x}_2\), no obstante que cuando ambas variables se examinan simultáneamente, el comportamiento inusual del punto es bastante evidente. Además, obsérvese que la probabilidad de que \(\bar{x}_1\) o bien \(\bar{x}_2\), excedan los límites de control tres sigma es 0,0027. Sin embargo, la probabilidad conjunta de que ambas variables excedan simultáneamente sus límites de control cuando ambas están bajo control es \((0,0027)(0,0027)=0,00000729\), que es considerablemente menor que \(0,0027\). Adicionalmente, la probabilidad de que tanto \(\bar{x}_1\) como \(\bar{x}_2\) se localicen simultáneamente dentro de los límites de control cuando el proceso en realidad está bajo control es \((0,9973)(0,9973)=0,99460729\). Por lo tanto, el uso de los diagramas \(\bar{x}\) independientes ha distorsionado el monitoreo simultáneo de \(\bar{x}_1\) y \(\bar{x}_2\), por cuanto el error tipo I y la probabilidad de que un punto se localice correctamente bajo control no son iguales a los niveles advertidos en los diagramas de control individuales.

Código 
```{r}
#| label: fig-dc_x1barra_x2barra
#| fig-cap: "Región de contol usando límites independiente para $\\bar{x}_1$ y $\\bar{x}_2$"
#| fig-subcap: 
#|   - "Diagrama $\\bar{x}_1$"
#|   - "Diagrama $\\bar{x}_2$"
#| layout-ncol: 2

vec_medias1 <- c(100, 110)
matriz_var_cov <- matrix(
  c(4, 0, 0, 4), nrow = 2 , byrow = TRUE
)
vec_medias2 <- c(101.5, 112)

library(mvtnorm)
set.seed(135711)
data_dnbv <- data.frame(
  mvtnorm::rmvnorm(
    n = 10000, mean = vec_medias1, sigma = matriz_var_cov
  )
) |>
  dplyr::slice_sample(n = 171) |>
  dplyr::bind_rows(
    data.frame(
      mvtnorm::rmvnorm(
        n = 9, mean = vec_medias2, sigma = matriz_var_cov
      )
    )
  ) |>
  dplyr::mutate(
    Muestra = rep(x = 1:20, each = 9), .before = X1
  ) |>
  # dplyr::relocate(Muestra, .before = X1) |>
  dplyr::group_by(Muestra) |>
  dplyr::summarize(
    n = dplyr::n(),
    X1_bar = mean(X1),
    X2_bar = mean(X2)
  )

######### Diagrama de control para la media de X1 ##########
LIC_1 <- qcc::limits.xbar(
  center = 100, std.dev = 2, sizes = 9, conf = 3
)[[1]]
LSC_1 <- qcc::limits.xbar(
  center = 100, std.dev = 2, sizes = 9, conf = 3
)[[2]]

dcxbar_1 <- ggplot(
  data = data_dnbv,
  aes(x = Muestra, y = X1_bar)
) +
  # geom_step(
  #   mapping = aes(x = Muestra, y = LIC_1),
  #   direction = "hv", color = "red", lty = "dashed"
  # ) +
  geom_step(
    mapping = aes(x = Muestra, y = 100),
    direction = "hv", color = "blue"
  ) +
  # geom_step(
  #   mapping = aes(x = Muestra, y = LSC_1),
  #   direction = "hv", color = "red", lty = "dashed"
  # ) +
  geom_point(
    mapping = aes(x = Muestra, y = X1_bar), color = "black"
  ) +
  geom_line(
    mapping = aes(x = Muestra, y = X1_bar), color = "black"
  ) +
  pammtools::geom_stepribbon(
    aes(ymin = LIC_1, ymax = 100), fill = "green",
    alpha = 0.2, direction = "hv"
  ) +
  pammtools::geom_stepribbon(
    aes(ymin = 100, ymax = LSC_1), fill = "green",
    alpha = 0.2, direction = "hv"
  ) +
  annotate(
    "text", x = 20 + 0.9, y = LIC_1,
    label = latex2exp::TeX("$LIC_{\\bar{x}_1}$", bold = TRUE),
    colour = "red"
  ) +
  annotate(
    "text", x = 20 + 0.75, y = 100,
    label = latex2exp::TeX("$LC_{\\bar{x}_1}$", bold = TRUE),
    colour = "blue"
  ) +
  annotate(
    "text", x = 20 + 0.9, y = LSC_1,
    label = latex2exp::TeX("$LSC_{\\bar{x}_1}$", bold = TRUE),
    colour = "red"
  ) +
  labs(
    x = latex2exp::TeX(
      "Muestra", bold = TRUE, italic = TRUE
    ),
    y = latex2exp::TeX(
      "$\\bar{x}_1$", bold = TRUE, italic = TRUE
    )
  ) +
  scale_x_continuous(
    breaks = seq(from = 1, to = 20, by = 2),
    limits = c(1, 22)
  ) +
  theme_bw() +
  theme(
    axis.title.y = element_text(
      angle = 0, vjust = 0.5, hjust = 0.5
    )
  )
dcxbar_1

# dcxbar_1 +
#   theme_minimal() +
#   theme(
#     axis.text.x = element_text(face = "bold"),
#     axis.text.y = element_text(face = "bold"),
#     axis.title.x = element_text(face = "bold"),
#     axis.title.y = element_text(face = "bold")
#   )
  
######### Diagrama de control para la media de X2 ##########
LIC_2 <- qcc::limits.xbar(
  center = 110, std.dev = 2, sizes = 9, conf = 3
)[[1]]
LSC_2 <- qcc::limits.xbar(
  center = 110, std.dev = 2, sizes = 9, conf = 3
)[[2]]

dcxbar_2 <- ggplot(
  data = data_dnbv,
  aes(x = Muestra, y = X2_bar)
) +
  geom_step(
    mapping = aes(x = Muestra, y = LIC_2),
    direction = "hv", color = "red", lty = "dashed"
  ) +
  geom_step(
    mapping = aes(x = Muestra, y = 110),
    direction = "hv", color = "blue"
  ) +
  geom_step(
    mapping = aes(x = Muestra, y = LSC_2),
    direction = "hv", color = "red", lty = "dashed"
  ) +
  geom_point(
    mapping = aes(x = Muestra, y = X2_bar), color = "black"
  ) +
  geom_line(
    mapping = aes(x = Muestra, y = X2_bar), color = "black"
  ) +
  pammtools::geom_stepribbon(
    aes(ymin = LIC_2, ymax = 110), fill = "green",
    alpha = 0.2, direction = "hv"
  ) +
  pammtools::geom_stepribbon(
    aes(ymin = 110, ymax = LSC_2), fill = "green",
    alpha = 0.2, direction = "hv"
  ) +
  annotate(
    "text", x = 20 + 0.9, y = LIC_2,
    label = latex2exp::TeX("$LIC_{\\bar{x}_2}$", bold = TRUE),
    colour = "red"
  ) +
  annotate(
    "text", x = 20 + 0.75, y = 110,
    label = latex2exp::TeX("$LC_{\\bar{x}_2}$", bold = TRUE),
    colour = "blue"
  ) +
  annotate(
    "text", x = 20 + 0.9, y = LSC_2,
    label = latex2exp::TeX("$LSC_{\\bar{x}_2}$", bold = TRUE),
    colour = "red"
  ) +
  labs(
    x = latex2exp::TeX(
      "Muestra", bold = TRUE, italic = TRUE
    ),
    y = latex2exp::TeX(
      "$\\bar{x}_2$", bold = TRUE, italic = TRUE
    )
  ) +
  scale_x_continuous(
    breaks = seq(from = 1, to = 20, by = 2),
    limits = c(1, 22)
  ) +
  theme_bw() +
  theme(
    axis.title.y = element_text(
      angle = 0, vjust = 0.5, hjust = 0.5
    )
  )
dcxbar_2

# dcxbar_2 +
#   theme_minimal() +
#   theme(
#     axis.text.x = element_text(face = "bold"),
#     axis.text.y = element_text(face = "bold"),
#     axis.title.x = element_text(face = "bold"),
#     axis.title.y = element_text(
#       face = "bold", angle = 0, vjust = 0.5, hjust = 0.5
#     )
#   )
```
(a) Diagrama \(\bar{x}_1\)
(b) Diagrama \(\bar{x}_2\)
Figura 5.20: Región de contol usando límites independiente para \(\bar{x}_1\) y \(\bar{x}_2\)
Código 
```{r}
#| label: fig-gd_xbarra1_xbarra2
#| fig-cap: "Diagramas de control $\\bar{x}$ independientes para las variables $X_1$ y $X_2$"

##################### Diagrama de dispersión ##################

# gd <- ggplot(
#   data = data_dnbv,
#   aes(x = X1_bar, y = X2_bar)
# ) +
#   # ggforce::geom_circle(
#   #   aes(
#   #     x0 = 100, y0 = 105,
#   #     r = qchisq(p = (1 - 0.0027)^2, df = 2)
#   #   ),
#   #   fill = "green", alpha = 0.01, color = "red"
#   # ) +
#   coord_fixed() +
#   geom_point(
#     color = ifelse(data_dnbv$Muestra <= 19, "black", "red"),
#     shape = ifelse(data_dnbv$Muestra <= 19, 16, 8),
#     size = 2.5
#   ) +
#   labs(
#     x = latex2exp::TeX(
#       "$\\bar{x}_1$",
#       bold = TRUE, italic = TRUE
#     ),
#     y = latex2exp::TeX(
#       "$\\bar{x}_2$",
#       bold = TRUE, italic = TRUE
#     )
#   ) +
#   theme_minimal() +
#   scale_y_continuous(position = "right") +
#   theme(
#     axis.title.y = element_text(
#       angle = 90, vjust = 0.5, hjust = 0.5
#     ),
#     axis.text.y = element_text(hjust = 0.5)
#   ) +
#   # geom_hline(
#   #   yintercept = c(LIC_2, LSC_2), color = "red",
#   #   linetype = "dashed", linewidth = 0.5
#   # ) +
#   # geom_vline(
#   #   xintercept = c(LIC_1, LSC_1), color = "red",
#   #   linetype = "dashed", linewidth = 0.5
#   # ) +
#   geom_ribbon(
#     aes(
#       x = seq(from = LIC_1, to = LSC_1, length.out = 20),
#         ymin = LIC_2, ymax = LSC_2
#     ),
#     fill = "green", alpha = 0.2
#   )
# 
# gd + theme_void()

knitr::include_graphics(path = "figuras/multivar.png")
```
Figura 5.21: Diagramas de control \(\bar{x}\) independientes para las variables \(X_1\) y \(X_2\)

Esta distorsión del procedimiento de monitoreo del proceso se incrementa conforme aumenta el número de características de la calidad. En general, si hay \(p\) características de la calidad estadísticamente independientes para un producto particular y si se mantiene un diagrama \(\bar{x}\) con \(P\left\{ \text{error tipo I} \right\} = \alpha\) en cada una, entonces la verdadera probabilidad del error tipo I para el procedimiento de control conjunto es

Probabilidad del error tipo I muultivariado

\[ \begin{equation} \alpha^{\prime} = 1 - (1 - \alpha)^p \end{equation}\:, \tag{5.30}\]

y, por complemento, la probabilidad de que las \(p\) medias se localicen simultáneamente dentro de sus límites de control cunado el proceso está bajo control es

\[ \begin{equation} P\left\{ \text{todas las} \: p \:\text{mdias se localicen bajo control} \right\} = \left( 1 - \alpha \right)^p \end{equation}\:, \tag{5.31}\]

Evidentemente, la distorsión del procedimiento de control conjunto puede ser grave, incluso para valores moderados de \(p\). Además, si las \(p\) características de la calidad no son independientes, lo cual sería generalmente el caso si estuvieran relacionadas con el mismo producto, las ecuaciones ecuación 5.30 y ecuación 5.31 no son válidas, y no se cuenta con una manera sencilla de medir la distorsión en el procedimiento de control conjunto.

A los problemas de monitoreo del proceso en que son varias las variables relacionadas de interés, en ocasiones se les llama problemas de control de calidad (o de monitoreo del proceso) con variables múltiples. Los trabajos originales en el control de calidad de variables múltiples fueron realizados por Hotelling (1947), quien aplicó sus procedimientos a datos de miras de bombardeo durante la Segunda Guerra mundial. Escritos posteriores que tratan los procedimientos de control para varias variables relacionadas incluyen a Hicks (1955); Jackson (1956; 1959, 1985); Crosier (1988); Hawkins (1991, 1993b); Lowry y otros (1992); Lowry y Montgomery (1995); Pignatiello y Runger (1990); Tracy, Young, and Mason (1992); Montgomery y Wadsworth (1972); y Alt (1985).

5.5.1 El Diagrama de Control \(T^2\) de Hotelling

El diagrama de control \(T^2\) de Hotelling es el procedimiento de monitoreo y control de procesos multivariables más ampliamente usado para para monitorear el vector de medias del proceso. Este diagrama es equivalente al diagrama \(\bar{x}\) de Shewhart para una variable. Se presentas dos versiones del diagrama \(T^2\) de Hotelling: una para datos subagrupados y otra para observaciones individuales.

5.5.1.1 El Diagrama de Control \(T^2\) de Hotelling para Datos Subagrupados

Supóngase que se tienen dos características de la calidad \(X_1\) y \(x_2\) que tienen distribución normal bivariadas. Sean \(\mu_1\) y \(\mu_2\) los valores medios de las características de la calidad, y \(\sigma_1\) y \(\sigma_2\) las desviaciones estándar de \(X_1\) y \(X_2\), respectivamente. La covarianza entre \(X_1\) y \(X_2\) se denota por \(\sigma_{12}\). Si se supone que \(\sigma_1\), \(\sigma_2\) y \(\sigma_{12}\) son conocidas. Si \(\bar{x}_1\) y \(\bar{x}_2\) son los promedios muestrales de las dos características de la calidad calculados a partir de una muestra de tamaño \(n\), entonces el estadístico

\[ \begin{align*} \chi_{0}^{2} = & \frac{n}{\sigma_{1}^{2}\sigma_{2}^{2} - \sigma_{12}^{2}}\Bigl[ \sigma_{2}^{2}\left( {\bar{x}}_{1} - \mu_{1} \right)^2 + \sigma_{1}^{2}\left( {\bar{x}}_{2} - \mu_{2} \right)^2 \\ &-2\sigma_{12}\left( {\bar{x}}_{1} - \mu_{1} \right)\left( {\bar{x}}_{2} - \mu_{2} \right) \Bigr] \end{align*} \tag{5.32}\]

tendrá una distribución ji-cuadrada con 2 grados de libertad. Esta ecuación puede usarse como base de un diagrama de control para las medias del proceso \(\mu_1\) y \(\mu_2\). Si las medias del proceso se mantienen en los valores \(\mu_1\) y \(\mu_2\), entonces los valores \(\chi_{0}^{2}\) deberán ser menores que el límite de control superior \(LSC = {\chi}_{_{1 - \alpha, 2}}^2\), donde \({\chi}_{_{1 - \alpha, 2}}^2\) es el percentil \(1 - \alpha\) de una variable aleatoria ji-cuadrada con dos grados de libertad. Si al menos una de las medias se corre a algún valor nuevo (fuera de control), entonces la probabilidad de que el estadístico \(\chi_{0}^{2}\) exceda el límite de control superior se incrementa.

Código 
```{r}
#| label: fig-elipse_vai
#| fig-cap: "Elipse de control para dos variables independientes"


# vec_medias1 <- c(100, 110)
# matriz_var_cov <- matrix(
#   c(4, 0, 0, 4), nrow = 2 , byrow = TRUE
# )
# vec_medias2 <- c(101.5, 112)
# 
# library(mvtnorm)
# set.seed(135711)
# data_dnbv_vai <- data.frame(
#   mvtnorm::rmvnorm(
#     n = 10000, mean = vec_medias1, sigma = matriz_var_cov
#   )
# ) |>
#   dplyr::slice_sample(n = 171) |>
#   dplyr::bind_rows(
#     data.frame(
#       mvtnorm::rmvnorm(
#         n = 9, mean = vec_medias2, sigma = matriz_var_cov
#       )
#     )
#   ) |>
#   dplyr::mutate(
#     Muestra = rep(x = 1:20, each = 9)
#   ) |>
#   dplyr::relocate(Muestra, .before = X1) |>
#   dplyr::group_by(Muestra) |>
#   dplyr::summarize(
#     n = dplyr::n(),
#     X1_bar = mean(X1),
#     X2_bar = mean(X2)
#   )
# 
# ######### Diagrama de control para la media de X1 ##########
# 
# LIC_1 <- qcc::limits.xbar(
#   center = 100, std.dev = 2, sizes = 9, conf = 3
# )[[1]]
# LSC_1 <- qcc::limits.xbar(
#   center = 100, std.dev = 2, sizes = 9, conf = 3
# )[[2]]
# 
# dcxbar_1 <- ggplot(
#   data = data_dnbv_vai,
#   aes(x = Muestra, y = X1_bar)
# ) +
#   # geom_step(
#   #   mapping = aes(x = Muestra, y = LIC_1),
#   #   direction = "hv", color = "red", lty = "dashed"
#   # ) +
#   geom_step(
#     mapping = aes(x = Muestra, y = 100),
#     direction = "hv", color = "blue"
#   ) +
#   # geom_step(
#   #   mapping = aes(x = Muestra, y = LSC_1),
#   #   direction = "hv", color = "red", lty = "dashed"
#   # ) +
#   geom_point(
#     mapping = aes(x = Muestra, y = X1_bar), color = "black"
#   ) +
#   geom_line(
#     mapping = aes(x = Muestra, y = X1_bar), color = "black"
#   ) +
#   pammtools::geom_stepribbon(
#     aes(ymin = LIC_1, ymax = 100), fill = "green",
#     alpha = 0.2, direction = "hv"
#   ) +
#   pammtools::geom_stepribbon(
#     aes(ymin = 100, ymax = LSC_1), fill = "green",
#     alpha = 0.2, direction = "hv"
#   ) +
#   annotate(
#     "text", x = 20 + 0.9, y = LIC_1,
#     label = latex2exp::TeX("$LIC_{\\bar{x}_1}$", bold = TRUE),
#     colour = "red"
#   ) +
#   annotate(
#     "text", x = 20 + 0.75, y = 100,
#     label = latex2exp::TeX("$LC_{\\bar{x}_1}$", bold = TRUE),
#     colour = "blue"
#   ) +
#   annotate(
#     "text", x = 20 + 0.9, y = LSC_1,
#     label = latex2exp::TeX("$LSC_{\\bar{x}_1}$", bold = TRUE),
#     colour = "red"
#   ) +
#   labs(
#     x = latex2exp::TeX(
#       "Muestra", bold = TRUE, italic = TRUE
#     ),
#     y = latex2exp::TeX(
#       "$\\bar{x}_1$", bold = TRUE, italic = TRUE
#     )
#   ) +
#   scale_x_continuous(
#     breaks = seq(from = 1, to = 20, by = 2),
#     limits = c(1, 21)
#   ) +
#   theme_bw() +
#   theme(
#     axis.title.y = element_text(
#       angle = 0, vjust = 0.5, hjust = 0.5
#     )
#   )
# dcxbar_1
# 
# dcxbar_1 +
#   theme_minimal() +
#   theme(
#     axis.text.x = element_text(face = "bold"),
#     axis.text.y = element_text(face = "bold"),
#     axis.title.x = element_text(face = "bold"),
#     axis.title.y = element_text(face = "bold")
#   )
# 
# ########## Diagrama de control para la media de X2 #########
# 
# LIC_2 <- qcc::limits.xbar(
#   center = 110, std.dev = 2, sizes = 9, conf = 3
# )[[1]]
# LSC_2 <- qcc::limits.xbar(
#   center = 110, std.dev = 2, sizes = 9, conf = 3
# )[[2]]
# 
# dcxbar_2 <- ggplot(
#   data = data_dnbv_vai,
#   aes(x = Muestra, y = X2_bar)
# ) +
#   # geom_step(
#   #   mapping = aes(x = Muestra, y = LIC_2),
#   #   direction = "hv", color = "red", lty = "dashed"
#   # ) +
#   geom_step(
#     mapping = aes(x = Muestra, y = 110),
#     direction = "hv", color = "blue"
#   ) +
#   # geom_step(
#   #   mapping = aes(x = Muestra, y = LSC_2),
#   #   direction = "hv", color = "red", lty = "dashed"
#   # ) +
#   geom_point(
#     mapping = aes(x = Muestra, y = X2_bar), color = "black"
#   ) +
#   geom_line(
#     mapping = aes(x = Muestra, y = X2_bar), color = "black"
#   ) +
#   pammtools::geom_stepribbon(
#     aes(ymin = LIC_2, ymax = 110), fill = "green",
#     alpha = 0.2, direction = "hv"
#   ) +
#   pammtools::geom_stepribbon(
#     aes(ymin = 110, ymax = LSC_2), fill = "green",
#     alpha = 0.2, direction = "hv"
#   ) +
#   annotate(
#     "text", x = 20 + 0.9, y = LIC_2,
#     label = latex2exp::TeX("$LIC_{\\bar{x}_2}$", bold = TRUE),
#     colour = "red"
#   ) +
#   annotate(
#     "text", x = 20 + 0.75, y = 110,
#     label = latex2exp::TeX("$LC_{\\bar{x}_2}$", bold = TRUE),
#     colour = "blue"
#   ) +
#   annotate(
#     "text", x = 20 + 0.9, y = LSC_2,
#     label = latex2exp::TeX("$LSC_{\\bar{x}_2}$", bold = TRUE),
#     colour = "red"
#   ) +
#   labs(
#     x = latex2exp::TeX(
#       "Muestra", bold = TRUE, italic = TRUE
#     ),
#     y = latex2exp::TeX(
#       "$\\bar{x}_2$", bold = TRUE, italic = TRUE
#     )
#   ) +
#   scale_x_continuous(
#     breaks = seq(from = 1, to = 20, by = 2),
#     limits = c(1, 21)
#   ) +
#   theme_bw() +
#   theme(
#     axis.title.y = element_text(
#       angle = 0, vjust = 0.5, hjust = 0.5
#     )
#   )
# dcxbar_2
# 
# dcxbar_2 +
#   theme_minimal() +
#   theme(
#     axis.text.x = element_text(face = "bold"),
#     axis.text.y = element_text(face = "bold"),
#     axis.title.x = element_text(face = "bold"),
#     axis.title.y = element_text(
#       face = "bold", angle = 0, vjust = 0.5, hjust = 0.5
#     )
#   )
# 
# ############# Elipse de control #################
# n <- 9
# mu_1 <- vec_medias1[1]
# mu_2 <- vec_medias1[2]
# angulo <- 0
# var_1 <- matriz_var_cov[1, 1]
# var_2 <- matriz_var_cov[2, 2]
# covar <- matriz_var_cov[1, 2]
# chi2 <- qchisq(p = (1 - 0.0027)^2, df = 2)
# a <- max(
#   sqrt(((var_1 * var_2 - covar^2) * chi2) / (n * var_2)),
#   sqrt(((var_1 * var_2 - covar^2) * chi2) / (n * var_1))
# )
# b <- min(
#   sqrt(((var_1 * var_2 - covar^2) * chi2) / (n * var_2)),
#   sqrt(((var_1 * var_2 - covar^2) * chi2) / (n * var_1))
# )
# 
# ec_var_ind <- ggplot(
#   data = data_dnbv_vai,
#   aes(x = X1_bar, y = X2_bar)
# ) +
#   # ggforce::geom_circle(
#   #   aes(
#   #     x0 = mu_1, y0 = mu_2,
#   #     r = sqrt(var_1 * chi2 / n)
#   #   ),
#   #   fill = "green", alpha = 0.01, color = "red"
#   # ) +
#   ggforce::geom_ellipse(
#     aes(x0 = mu_1, y0 = mu_2, a = a, b = b, angle = 0),
#     fill = "green", alpha = 0.2, color = "red"
#   ) +
#   coord_fixed() +
#   geom_point(
#     color = ifelse(
#       data_dnbv_vai$Muestra <= 19, "black", "red"
#     ),
#     shape = ifelse(data_dnbv_vai$Muestra <= 19, 16, 8),
#     size = 2.5
#   ) +
#   labs(
#     x = latex2exp::TeX(
#       "$\\bar{x}_1$",
#       bold = TRUE, italic = TRUE
#     ),
#     y = latex2exp::TeX(
#       "$\\bar{x}_2$",
#       bold = TRUE, italic = TRUE
#     )
#   ) +
#   theme_minimal() +
#   scale_y_continuous(position = "right") +
#   theme(
#     axis.title.y = element_text(
#       angle = 90, vjust = 0.5, hjust = 0.5
#     ),
#     axis.text.y = element_text(hjust = 0.5)
#   ) +
#   # geom_hline(
#   #   yintercept = c(LIC_2, LSC_2), color = "red",
#   #   linetype = "dashed", linewidth = 0.5
#   # ) +
#   # geom_vline(
#   #   xintercept = c(LIC_1, LSC_1), color = "red",
#   #   linetype = "dashed", linewidth = 0.5
#   # ) +
#   geom_ribbon(
#     aes(
#       x = seq(from = LIC_1, to = LSC_1, length.out = 20),
#         ymin = LIC_2, ymax = LSC_2
#     ),
#     fill = "green", alpha = 0.2
#   )
# 
# ec_var_ind +
#   theme_void()

knitr::include_graphics(path = "figuras/multivar_ind.png")
```
Figura 5.22: Elipse de control para dos variables independientes

El procedimiento de monitoreo del proceso puede representarse gráficamente. Considere el caso en que las dos variables aleatorias \(X_1\) y \(X_2\) son independientes; es decir, \(\sigma_{12} = 0\). Entonces, la ecuación ecuación 5.32 define una elipse con centro en \(\left( \mu_{0},\, \mu_{1} \right)\) y ejes principales paralelos a los ejes \(\left( \bar{x}_1,\, \bar{x}_2 \right)\), como se muestra en la Figura 5.22. Tomar \(\chi_{0}^{2}\) igual a \({\chi}_{_{1 - \alpha, 2}}^2\) en la ecuación ecuación 5.32 implica que un par de promedios muestrales \(\left( \bar{x}_1,\, \bar{x}_2 \right)\) que produzcan un valor de \(\chi_{0}^{2}\) que se localice dentro de la elipse indica que el proceso está bajo control, mientras que si el valor correspondiente de \(\chi_{0}^{2}\) se localiza fuera de la elipse, el proceso está fuera de control. A la Figura 5.22 acostumbra llamársele la elipse de control.

Cuando las dos características de la calidad están relacionadas (\(\sigma_{12} \neq 0\)) los ejes principales de la elipse no son paralelos a los ejes cartesianos \(\bar{x}_1\), \(\bar{x}_2\). La elipse de control que describe esta situación se muestra en la Figura 5.23. Esto es importante, porque ayuda a detectar problemas que podrían no ser evidentes al analizar los gráficos de control de cada variable por separado. En el ejemplo dado, el punto muestral número 20 se encuentra fuera de la elipse de control, lo que indica una posible causa asignable. Sin embargo, este punto está dentro de los límites de control de los gráficos individuales para \(\bar{x}_1\) y \(\bar{x}_2\). Esto significa que, aunque los gráficos individuales no muestren nada inusual, la elipse de control combinada revela un posible problema. Esto evidencia la importancia de considerar la relación entre las variables para identificar problemas que de otro modo podrían pasar desapercibidos.

Código 
```{r}
#| label: fig-elipse_vad
#| fig-cap: "Elipse de control para dos variables dependientes"

# vec_medias1 <- c(100, 110)
# matriz_var_cov <- matrix(
#   c(4, 3.6, 3.6, 9), nrow = 2 , byrow = TRUE
# )
# vec_medias2 <- c(101.7, 113.2)
# 
# library(mvtnorm)
# set.seed(135711)
# data_dnbv_vai <- data.frame(
#   mvtnorm::rmvnorm(
#     n = 10000, mean = vec_medias1, sigma = matriz_var_cov
#   )
# ) |>
#   dplyr::slice_sample(n = 171) |>
#   dplyr::bind_rows(
#     data.frame(
#       mvtnorm::rmvnorm(
#         n = 9, mean = vec_medias2, sigma = matriz_var_cov
#       )
#     )
#   ) |>
#   dplyr::mutate(
#     Muestra = rep(x = 1:20, each = 9, .before = X1)
#   ) |>
#   dplyr::group_by(Muestra) |>
#   dplyr::summarize(
#     n = dplyr::n(),
#     X1_bar = mean(X1),
#     X2_bar = mean(X2)
#   )
# 
# ######### Diagrama de control para la media de X1 ##########
# 
# LIC_1 <- qcc::limits.xbar(
#   center = 100, std.dev = 2, sizes = 9, conf = 3
# )[[1]]
# LSC_1 <- qcc::limits.xbar(
#   center = 100, std.dev = 2, sizes = 9, conf = 3
# )[[2]]
# 
# dcxbar_1 <- ggplot(
#   data = data_dnbv_vai,
#   aes(x = Muestra, y = X1_bar)
# ) +
#   # geom_step(
#   #   mapping = aes(x = Muestra, y = LIC_1),
#   #   direction = "hv", color = "red", lty = "dashed"
#   # ) +
#   geom_step(
#     mapping = aes(x = Muestra, y = 100),
#     direction = "hv", color = "blue"
#   ) +
#   # geom_step(
#   #   mapping = aes(x = Muestra, y = LSC_1),
#   #   direction = "hv", color = "red", lty = "dashed"
#   # ) +
#   geom_point(
#     mapping = aes(x = Muestra, y = X1_bar), color = "black"
#   ) +
#   geom_line(
#     mapping = aes(x = Muestra, y = X1_bar), color = "black"
#   ) +
#   pammtools::geom_stepribbon(
#     aes(ymin = LIC_1, ymax = 100), fill = "green",
#     alpha = 0.2, direction = "hv"
#   ) +
#   pammtools::geom_stepribbon(
#     aes(ymin = 100, ymax = LSC_1), fill = "green",
#     alpha = 0.2, direction = "hv"
#   ) +
#   annotate(
#     "text", x = 20 + 0.9, y = LIC_1,
#     label = latex2exp::TeX("$LIC_{\\bar{x}_1}$", bold = TRUE),
#     colour = "red"
#   ) +
#   annotate(
#     "text", x = 20 + 0.75, y = 100,
#     label = latex2exp::TeX("$LC_{\\bar{x}_1}$", bold = TRUE),
#     colour = "blue"
#   ) +
#   annotate(
#     "text", x = 20 + 0.9, y = LSC_1,
#     label = latex2exp::TeX("$LSC_{\\bar{x}_1}$", bold = TRUE),
#     colour = "red"
#   ) +
#   labs(
#     x = latex2exp::TeX(
#       "Muestra", bold = TRUE, italic = TRUE
#     ),
#     y = latex2exp::TeX(
#       "$\\bar{x}_1$", bold = TRUE, italic = TRUE
#     )
#   ) +
#   scale_x_continuous(
#     breaks = seq(from = 1, to = 20, by = 2),
#     limits = c(1, 21)
#   ) +
#   theme_bw() +
#   theme(
#     axis.title.y = element_text(
#       angle = 0, vjust = 0.5, hjust = 0.5
#     )
#   )
# dcxbar_1
# 
# dcxbar_1 +
#   theme_minimal() +
#   theme(
#     axis.text.x = element_text(face = "bold"),
#     axis.text.y = element_text(face = "bold"),
#     axis.title.x = element_text(face = "bold"),
#     axis.title.y = element_text(face = "bold")
#   )
# 
# ########## Diagrama de control para la media de X2 #########
# 
# LIC_2 <- qcc::limits.xbar(
#   center = 110, std.dev = 3, sizes = 9, conf = 3
# )[[1]]
# LSC_2 <- qcc::limits.xbar(
#   center = 110, std.dev = 3, sizes = 9, conf = 3
# )[[2]]
# 
# dcxbar_2 <- ggplot(
#   data = data_dnbv_vai,
#   aes(x = Muestra, y = X2_bar)
# ) +
#   # geom_step(
#   #   mapping = aes(x = Muestra, y = LIC_2),
#   #   direction = "hv", color = "red", lty = "dashed"
#   # ) +
#   geom_step(
#     mapping = aes(x = Muestra, y = 110),
#     direction = "hv", color = "blue"
#   ) +
#   # geom_step(
#   #   mapping = aes(x = Muestra, y = LSC_2),
#   #   direction = "hv", color = "red", lty = "dashed"
#   # ) +
#   geom_point(
#     mapping = aes(x = Muestra, y = X2_bar), color = "black"
#   ) +
#   geom_line(
#     mapping = aes(x = Muestra, y = X2_bar), color = "black"
#   ) +
#   pammtools::geom_stepribbon(
#     aes(ymin = LIC_2, ymax = 110), fill = "green",
#     alpha = 0.2, direction = "hv"
#   ) +
#   pammtools::geom_stepribbon(
#     aes(ymin = 110, ymax = LSC_2), fill = "green",
#     alpha = 0.2, direction = "hv"
#   ) +
#   annotate(
#     "text", x = 20 + 0.9, y = LIC_2,
#     label = latex2exp::TeX("$LIC_{\\bar{x}_2}$", bold = TRUE),
#     colour = "red"
#   ) +
#   annotate(
#     "text", x = 20 + 0.75, y = 110,
#     label = latex2exp::TeX("$LC_{\\bar{x}_2}$", bold = TRUE),
#     colour = "blue"
#   ) +
#   annotate(
#     "text", x = 20 + 0.9, y = LSC_2,
#     label = latex2exp::TeX("$LSC_{\\bar{x}_2}$", bold = TRUE),
#     colour = "red"
#   ) +
#   labs(
#     x = latex2exp::TeX(
#       "Muestra", bold = TRUE, italic = TRUE
#     ),
#     y = latex2exp::TeX(
#       "$\\bar{x}_2$", bold = TRUE, italic = TRUE
#     )
#   ) +
#   scale_x_continuous(
#     breaks = seq(from = 1, to = 20, by = 2),
#     limits = c(1, 21)
#   ) +
#   theme_bw() +
#   theme(
#     axis.title.y = element_text(
#       angle = 0, vjust = 0.5, hjust = 0.5
#     )
#   )
# dcxbar_2
# 
# dcxbar_2 +
#   theme_minimal() +
#   theme(
#     axis.text.x = element_text(face = "bold"),
#     axis.text.y = element_text(face = "bold"),
#     axis.title.x = element_text(face = "bold"),
#     axis.title.y = element_text(
#       face = "bold", angle = 0, vjust = 0.5, hjust = 0.5
#     )
#   )
# 
# ############# Elipse de control #################
# n <- 9
# mu_1 <- vec_medias1[1]
# mu_2 <- vec_medias1[2]
# var_1 <- matriz_var_cov[1, 1]
# var_2 <- matriz_var_cov[2, 2]
# covar <- matriz_var_cov[1, 2]
# angulo <- 0.5 * atan(2 * covar / (var_2 - var_1))
# chi2 <- qchisq(p = (1 - 0.0027)^2, df = 2)
# a <-  max(
#   sqrt(((var_1 * var_2 - covar^2) * chi2) / (n * var_2)),
#   sqrt(((var_1 * var_2 - covar^2) * chi2) / (n * var_1))
# )
# b <-  min(
#   sqrt(((var_1 * var_2 - covar^2) * chi2) / (n * var_2)),
#   sqrt(((var_1 * var_2 - covar^2) * chi2) / (n * var_1))
# )
# 
# ec_var_ind <- ggplot(
#   data = data_dnbv_vai,
#   aes(x = X1_bar, y = X2_bar)
# ) +
#   ggforce::geom_ellipse(
#     aes(
#       x0 = mu_1, y0 = mu_2, a = a , b = b,
#       angle = angulo
#     ),
#     fill = "green", alpha = 0.2, color = "red"
#   ) +
#   coord_fixed() +
#   geom_point(
#     color = ifelse(
#       data_dnbv_vai$Muestra <= 19, "black", "red"
#     ),
#     shape = ifelse(data_dnbv_vai$Muestra <= 19, 16, 8),
#     size = 2.5
#   ) +
#   labs(
#     x = latex2exp::TeX(
#       "$\\bar{x}_1$",
#       bold = TRUE, italic = TRUE
#     ),
#     y = latex2exp::TeX(
#       "$\\bar{x}_2$",
#       bold = TRUE, italic = TRUE
#     )
#   ) +
#   theme_minimal() +
#   scale_y_continuous(position = "right") +
#   theme(
#     axis.title.y = element_text(
#       angle = 90, vjust = 0.5, hjust = 0.5
#     ),
#     axis.text.y = element_text(hjust = 0.5)
#   ) +
#   # geom_hline(
#   #   yintercept = c(LIC_2, LSC_2), color = "red",
#   #   linetype = "dashed", linewidth = 0.5
#   # ) +
#   # geom_vline(
#   #   xintercept = c(LIC_1, LSC_1), color = "red",
#   #   linetype = "dashed", linewidth = 0.5
#   # ) +
#   geom_ribbon(
#     aes(
#       x = seq(from = LIC_1, to = LSC_1, length.out = 20),
#         ymin = LIC_2, ymax = LSC_2
#     ),
#     fill = "green", alpha = 0.2
#   )
# 
# ec_var_ind +
#   theme_void()

knitr::include_graphics(path = "figuras/elipse_control.png")
```
Figura 5.23: Elipse de control para dos variables dependientes

Hay dos desventajas asociadas con la elipse de control. La primera es que se pierde la secuencia temporal de los puntos graficados. Esto se puede superar numerando los puntos trazados o utilizando símbolos de trazado especiales para representar las observaciones más recientes. La segunda y más grave desventaja es que es difícil construir la elipse para más de dos características de calidad. Para evitar estas dificultades, es costumbre trazar los valores de \(\chi_0^2\) calculados a partir de la ecuación ecuación 5.32 para cada muestra en un diagrama de control con solo un límite de control superior en \(\chi_{1 - \alpha,\: 2}^2\) como se muestra en la Figura 5.24. Este diagrama de control se llama generalmente diagrama de control ji-cuadrado. Note que la secuencia temporal de los datos se conserva con este diagrama de control, de modo que se pueden investigar corridas u otros patrones no aleatorios. Además, tiene la ventaja adicional de que el “estado” del proceso se caracteriza por un solo número (el valor de del estadístico \(\chi_0^2\)). Esto es particularmente útil cuando se monitorean de manera simultanea dos o más características de la calidad, y más aún, si estas están relacionadas, lo cual es muy común en los procesos producción.

Código 
```{r}
#| label: fig-ji-cuadrada_vad
#| fig-cap: "Diagrama de control ji-cuadrada con $p = 2$ características de la caridad"

vec_medias1 <- c(100, 110)
matriz_var_cov <- matrix(
  c(4, 3.6, 3.6, 9),
  nrow = 2, byrow = TRUE
)
vec_medias2 <- c(101.7, 113.2)
n <- 9
library(mvtnorm)
set.seed(135711)

data_ji_cuadrada <- data.frame(
  mvtnorm::rmvnorm(
    n = 10000, mean = vec_medias1, sigma = matriz_var_cov
  )
) |>
  dplyr::slice_sample(n = 171) |>
  dplyr::bind_rows(
    data.frame(
      mvtnorm::rmvnorm(
        n = 9, mean = vec_medias2, sigma = matriz_var_cov
      )
    )
  ) |>
  dplyr::mutate(
    Muestra = rep(x = 1:20, each = 9), .before = X1
  ) |>
  dplyr::group_by(Muestra) |>
  dplyr::summarize(
    n = dplyr::n(),
    X1_bar = mean(X1),
    X2_bar = mean(X2)
  ) |>
  dplyr::mutate(
    Ji_cuadrada = (
      n / (
        matriz_var_cov[1, 1] * matriz_var_cov[2, 2] -
          matriz_var_cov[1, 2]^2
      )

    ) *
      (
        matriz_var_cov[2, 2] * (X1_bar - vec_medias1[1])^2 +
          matriz_var_cov[1, 1] * (X2_bar - vec_medias1[2])^2 -
          2 * matriz_var_cov[1, 2] * (X1_bar - vec_medias1[1])
            * (X2_bar - vec_medias1[2])
      )
  )

LSC <- qchisq(p = (1 - 0.0027)^2, df = 2)

ggplot(
  data = data_ji_cuadrada,
  aes(x = Muestra, y = Ji_cuadrada)
) +
  geom_step(
    mapping = aes(x = Muestra, y = LSC),
    direction = "hv", color = "red", lty = "dashed"
  ) +
  geom_point(
    mapping = aes(x = Muestra, y = Ji_cuadrada),
    color = ifelse(
      data_ji_cuadrada$Ji_cuadrada <= LSC, "black", "red"
    )
  ) +
  geom_line(
    mapping = aes(x = Muestra, y = Ji_cuadrada),
    color = "black"
  ) +
  pammtools::geom_stepribbon(
    aes(ymin = LSC, ymax = 0),
    fill = "green",
    alpha = 0.2, direction = "hv"
  ) +
  annotate(
    "text",
    x = 20 + 0.9, y = LSC,
    label = latex2exp::TeX("$LSC$", bold = TRUE),
    colour = "red"
  ) +
  labs(
    x = latex2exp::TeX(
      "Muestra",
      bold = TRUE, italic = TRUE
    ),
    y = latex2exp::TeX(
      "$\\chi_{0}^{2}",
      bold = TRUE, italic = TRUE
    )
  ) +
  scale_x_continuous(
    breaks = seq(from = 1, to = 20, by = 2),
    limits = c(1, 21)
  ) +
  theme_bw() +
  theme(
    axis.title.y = element_text(
      angle = 0, vjust = 0.5, hjust = 0.5
    )
  )
```
Figura 5.24: Diagrama de control ji-cuadrada con \(p = 2\) características de la caridad

Es posible ampliar estos resultados al caso en que se controlan al mismo tiempo \(p\) características de la calidad relacionadas. Suponga que se tiene un proceso de producción en el cual se monitorean de manera simultanea \(p\) características de la calidad y que la distribución de probabilidad conjunta de las \(p\) características de la calidad es la distribución normal \(p\) variada. Si se toman \(m\) muestras de tamaño \(n\) de cada característica de la calidad, de tal manera que \(X_{ijk}\) representa la observación \(j\) en la muestra \(i\) de la característica de la calidad \(k\). El procedimiento requiere calcular la media de cada una de las \(m\) muestras para cada una de las \(p\) características de la calidad a partir de una muestra de tamaño \(n\). En base a lo anterior, se define el conjunto de medias de la muestra \(i\) para cada una las \(p\) características de la calidad a través de un vector columna de orden \(p \times 1\), como sigue

\[ \overline{\mathbf{X}}_i = \begin{bmatrix} \overline{X}_{i1} \\ \overline{X}_{i2} \\ \vdots \\ \overline{X}_{ik} \\ \vdots \\ \overline{X}_{i(p-1)} \\ \overline{X}_{ip} \end{bmatrix} \]

Donde \(\overline{X}_{ik}\) representa la media de las \(n\) observaciones en la muestra \(i\) tomada de la característica de la calidad \(k\), con \(i = 1, 2, \dots, m\) y \(k = 1, 2, \dots, p\)

El estadístico graficado en el diagrama de control ji-cuadrada para cada muestra es

Estadístico ji-cuagrada

\[ \begin{equation} \chi_{i}^{2} = n\left(\overline{\mathbf{X}}_i - \pmb{\mu} \right)^{\pmb{'}}\pmb{\Sigma }^{-1}\left(\overline{\mathbf{X}}_i - \pmb{\mu} \right) \end{equation} \tag{5.33}\]

donde

\[ \pmb{\mu } = \begin{bmatrix} \mu_{1} \\ \mu_{2} \\ \vdots \\ \mu_{k} \\ \vdots \\ \mu_{p-1} \\ \mu_{p} \end{bmatrix} \] es el vector de medias de orden \(p\times 1\) que contiene las medias bajo control de las \(p\) características de la calidad y

\[ \pmb{\Sigma} = \begin{bmatrix} \sigma_1^2 & \sigma_{12} & \cdots & \sigma_{1k} & \cdots & \sigma_{1(p-1)} & \sigma_{1p} \\ \sigma_{21} & \sigma_2^2 & \cdots & \sigma_{2k} & \cdots & \sigma_{2(p-1)} & \sigma_{2p} \\ \vdots & \vdots & \ddots & \vdots & \ddots & \vdots & \vdots \\ \sigma_{k1} & \sigma_{k2} & \cdots & \sigma_{k}^{2} & \cdots & \sigma_{k(p-1)} & \sigma_{kp} \\ \vdots & \vdots & \ddots & \vdots & \ddots & \vdots & \vdots \\ \sigma_{(p-1)1} & \sigma_{(p-1)2} & \cdots & \sigma_{(p-1)k} & \cdots & \sigma_{(p-1)}^2 & \sigma_{(p-1)p} \\ \sigma_{p1} & \sigma_{p2} & \cdots & \sigma_{pk} & \cdots & \sigma_{p(p-1)} & \sigma_{p}^2 \end{bmatrix}. \] es la matriz de varianzas-covarianzas de orden \(p \times p\).

Los límites de control del diagrama ji-cuadrada es

Límite de control del diagrama ji-cuadrada

\[ \begin{align*} LSC & = \chi_{1 - \alpha, p}^{2} \\ LIC & = 0 \end{align*} \:. \tag{5.34}\]

Estimación de \(\pmb{\mu }\) y \(\pmb{\Sigma}\)

En la práctica, generalmente \(\pmb{\mu}\) y \(\pmb{\Sigma}\) son desconocidos por lo que deben ser estimados a partir del análisis de \(m\) muestras preliminares de tamaño \(n\), tomadas cuando se supone que el proceso está bajo control. De manera general, la media de la \(i\)-ésima muestra en la \(k\)-ésima característica de la calidad se define por medio del siguiente estadístico

\[ \begin{equation} \overline{X}_{ik} = \frac{\sum_{j = 1}^{n}X_{ijk}}{n}, \; \text{para} \: i=1,2,\dots, m \: \text{y}\: k=1,2, \dots, p \end{equation}. \tag{5.35}\] donde \(X_{ijk}\) es la \(j\)-ésima observación en la \(i\)-ésima muestra de la \(k\)-ésima característica de la calidad.

Mientras que, la varianza de la \(i\)-ésima muestra en la \(k\)-ésima característica de la calidad se define por el siguiente estadístico:

\[ \begin{equation} S_{ik}^{2} = \frac{\sum_{j = 1}^{n}\left ( X_{ijk} - \overline{X}_{ik} \right )^2}{n - 1}, \; \text{para} \: i = 1, 2, \dots, m \: \text{y}\: k = 1, 2, \dots, p \end{equation} \tag{5.36}\]

Por otro lado, la covarianza entre la característica de la calidad \(r\) y la característica de la calidad \(s\) en la muestra \(i\) se define a través del siguiente estadístico

\[ \begin{equation} S_{irs} = \frac{\sum_{j = 1}^{n}\left (X_{ijr} - \overline{X}_{ir} \right )\left (X_{ijs} - \overline{X}_{is} \right )}{n - 1}, \; \text{para} \: i = 1, 2, \dots, m \: \text{y}\: r \neq s \end{equation}. \tag{5.37}\]

Al promediar las \(m\) medias muestrales para cada característica de la calidad obtenidas con la ecuación 5.35 se obtiene el estimador del vector de media bajo control. Como sigue,

\[ \begin{equation} \overline{\overline{X}}_{k} = \frac{\sum_{i = 1}^{m}\overline{X}_{ik}}{m}, \; \text{para} \: j = 1,2,\dots,p \end{equation}. \tag{5.38}\]

Mientras que, promediando las \(m\) varianzas y covarianzas para cada característica de la calidad obtenidas con las ecuaciones 5.36, y 5.37, por medio de las siguientes ecuaciones, resultan los elementos de la matriz de varianzas-covarianzas estimada:

\[ \begin{equation} \overline{S}_{i}^{2} = \frac{\sum_{i = 1}^{m}S_{ik}^{2}}{m}, \: \text{para}\:k = 1, 2, \dots, p \end{equation} \tag{5.39}\]

\[ \begin{equation} \overline{S}_{rs} = \frac{\sum_{i = 1}^{m}S_{irs}}{m}, \: \text{para}\:r \neq s \end{equation}. \tag{5.40}\]

A partir de los resultados anteriores se construye el estimador del vector de medias \(\pmb{\mu }\) y de la matriz de varianzas-covarianzas \(\pmb{\Sigma}\) cuando el proceso está bajo control, como se muestra a continuación:

Vector de medias \(\pmb{\mu}\) estimado

\[ \begin{equation} \pmb{\widehat{\mu}} = \overline{\overline{\mathbf{X}}} = \begin{bmatrix} \overline{\overline{X}}_1 \\ \overline{\overline{X}}_2 \\ \vdots \\ \overline{\overline{X}}_{k} \\ \vdots \\ \overline{\overline{X}}_{p-1} \\ \overline{\overline{X}}_{p} \end{bmatrix} \end{equation} \tag{5.41}\]

Matriz de varianzas-covarianzas \(\pmb{\Sigma}\) estimada

\[ \begin{equation} \pmb{\widehat{\Sigma}} = \mathbf{S} = \begin{bmatrix} \overline{S}_1^2 & \overline{S}_{12} & \cdots & \overline{S}_{1k} & \cdots & \overline{S}_{1(p-1)} & \overline{S}_{1p} \\ \overline{S}_{21} & \overline{S}_2^2 & \cdots & \overline{S}_{2k} & \cdots & \overline{S}_{2(p-1)} & \overline{S}_{2p} \\ \vdots & \vdots & \ddots & \vdots & \ddots & \vdots & \vdots \\ \overline{S}_{k1} & \overline{S}_{k2} & \cdots & \overline{S}_{k}^{2} & \cdots & \overline{S}_{k(p-1)} & \overline{S}_{kp} \\ \vdots & \vdots & \ddots & \vdots & \ddots & \vdots & \vdots \\ \overline{S}_{(p-1)1} & \overline{S}_{(p-1)2} & \cdots & \overline{S}_{(p-1)k} & \cdots & \overline{S}_{(p-1)}^2 & \overline{S}_{(p-1)p} \\ \overline{S}_{p1} & \overline{S}_{p2} & \cdots & \overline{S}_{pk} & \cdots & \overline{S}_{p(p-1)} & \overline{S}_{p}^2 \end{bmatrix} \end{equation} \tag{5.42}\]

Suponga ahora que se usa \(\overline{\overline{\mathbf{X}}}\) de la ecuación de la ecuación 5.41 como estimador del vector de medias del proceso bajo control y \(\mathbf{S}\) de la ecuación 5.42 como estimador de la matriz de varianzas-covarianzas \(\pmb{\Sigma}\). Si se sustituye \(\pmb{\mu}\) con \(\overline{\overline{\mathbf{X}}}\) y \(\pmb{\Sigma}\) con en la ecuación 5.33, el estadístico graficado queda reescrito de la siguiente manera

Estadístico \(T^2\) de Hotelling

\[ \begin{equation} T_{i}^{2} = n\left(\overline{\mathbf{X}}_i - \overline{\overline{\mathbf{X}}} \right)^{\pmb{'}}\mathbf{S}^{-1}\left(\overline{\mathbf{X}}_i - \overline{\overline{\mathbf{X}}} \right) \end{equation} \tag{5.43}\]

En esta forma, al procedimiento suele llamársele diagrama de control \(T^2\) de Hotelling. Este diagrama de control es direccionalmente invariante; es decir, su capacidad para detectar un cambio en el vector de medias solo depende de la magnitud del cambio, y no de su dirección.

Alt (1985) ha señalado que en las aplicaciones de control de la calidad de variables múltiples debe tenerse cuidado al seleccionar los límites de control para el estadístico \(T^2\) de Hotelling (ecuación 5.43) con base en la forma en que se esté usando el diagrama. Hay dos fases distintas en la utilización de los diagramas de control. La fase I es el uso del diagrama para establecer el control; es decir, para probar si el proceso se encontraba bajo control cuando se sacaron los \(m\) subgrupos preliminares y se calcularon los estadísticos \(\overline{\overline{\mathbf{X}}}\) y \(\mathbf{S}\). El objetivo de la fase I es obtener un conjunto de observaciones bajo control a fin de poder establecer los límites de control para la fase II, la cual consiste en el monitoreo de la producción futura. En ocasiones se llama análisis retrospectivo.

Los límites de control del diagrama de control \(T^2\) de Hotelling en la fase I están dados por

Límites de control del diagrama \(T^2\) de Hotelling en la fase I

\[ \begin{align*} LSC & = \frac{p\left(m -1 \right)\left(n - 1\right)}{mn- m - p +1}F_{1 - \alpha, \:p, \:mn - m - p +1} \\ LIC & = 0 \end{align*}\:. \tag{5.44}\]

En la fase II, cuando la carta se usa para monitorear la producción futura, los límites de control son los siguientes:

Límites de control del diagrama \(T^2\) de Hotelling en la fase II

\[ \begin{align*} LSC & = \frac{p\left(m + 1\right)\left(n - 1\right)}{mn- m - p + 1}F_{1 - \alpha, \:p, \:mn - m - p +1} \\ LIC & = 0 \end{align*} \tag{5.45}\]

Nota que el \(LSC\) en la ecuación 5.45 es simplemente el \(LSC\) en la ecuación 5.44 multiplicado por \((m + 1) / (m - 1)\), por lo que el \(LSC\) en la fase II será ligeramente mayor que el obtenido en la fase I. Esta diferencia se aproxima a cero para valores grandes de \(m\).

Cuando \(\pmb{\mu}\) y \(\pmb{\Sigma}\) se estiman a partir de un gran número de muestras preliminares, es habitual utilizar \(LSC = {\chi}_{1 - \alpha, \: p}^{2}\) como el límite de control superior tanto en la fase I como en la fase II. El análisis retrospectivo de las muestras preliminares para verificar el control estadístico y establecer límites de control también ocurre cuando se establecen diagramas de control univariados. Para el diagrama \(\bar{x}\), es un hecho muy conocido que si se usan \(m ≥ 20\) o \(25\) muestras preliminares, suele ser innecesario hacer la distinción entre los límites de control de la fase I y la fase II, ya que dichos límites prácticamente coincidirán. Sin embargo, deberá tenerse cuidado con los diagramas de control multivariados. En un artículo de revisión, Jensen et al. (2006) señalan que se requieren tamaños de muestra aún mayores para asegurar que el rendimiento de la longitud media de corrida (LMC) de la fase II sea realmente cercano a los valores anticipados. Ellos recomiendan usar tantas muestras de la fase I como sea posible para estimar los límites de la fase II. Sin embargo, deberá tenerse cuidado con los diagramas de control para variables múltiples.

Lowry y Montgomery (1995) muestran que en muchas situaciones se requieren un gran número de muestras preliminares antes de que los límites de control ji-cuadrada sean una buena aproximación de los límites de control de la fase II. Estos autores presentan tablas que indican el valor mínimo recomendado de \(m\) para tamaños de muestra de \(n = 3, \:5\) y \(10\) y para \(p = 2, \:3, \:4, \:5, \:10\) y \(20\) características de la calidad. Los valores recomendados de \(m\) son siempre superiores a \(20\) muestras preliminares, y a menudo más de \(50\) muestras. Jensen et al. (2006) observan que estos tamaños de muestra recomendados son probablemente demasiado pequeños. Por lo que sugieren el número de muestras tomados en la fase I sea de al menos \(200\) (\(m \geq 200\)) para asegurar que los límites de la fase II sean correctamente estimados.

Ejemplo 5.6 (Aplicación del Diagrama de Control \(T^2\) de Hotelling para Datos Subagrupados)
En una planta química dedicada a la producción de polímeros sintéticos , el proceso comienza con la alimentación de materias primas (como el etileno, polietileno o polipropileno), junto con catalizadores, en un reactor. Estas materias primas son sometidas a una reacción de polimerización, donde las moléculas pequeñas (monómeros) se unen para formar cadenas largas de polímeros. Esta reacción es altamente sensible a las condiciones operativas, especialmente a la temperatura y la presión, que deben mantenerse dentro de rangos específicos para garantizar la calidad del producto final.

La temperatura es una variable crítica porque afecta directamente la velocidad y la eficiencia de la reacción de polimerización. Si la temperatura es demasiado alta, la reacción puede volverse incontrolable, generando polímeros con propiedades no deseadas, como baja resistencia mecánica o inconsistencia en la estructura molecular. Por otro lado, si la temperatura es demasiado baja, la reacción puede ser demasiado lenta, reduciendo la eficiencia del proceso y aumentando los costos de producción. Además, un aumento excesivo de la temperatura puede dañar el equipo o incluso provocar riesgos de seguridad, como fugas o explosiones. Mientra que, la presión está directamente relacionada con la cantidad de gases en el reactor y su expansión térmica. A medida que la temperatura aumenta, los gases se expanden, lo que incrementa la presión en el reactor. Si la presión es demasiado alta, puede superar los límites de diseño del equipo, provocando fallas mecánicas o riesgos de explosión. Por el contrario, si la presión es demasiado baja, la reacción puede no ocurrir de manera eficiente, lo que resulta en un producto de baja calidad o en la interrupción del proceso. Por lo tanto, mantener la presión dentro de un rango óptimo es esencial para la seguridad y la eficiencia del proceso.

La temperatura y la presión están intrínsecamente relacionadas en este proceso. Un aumento en la temperatura provoca un aumento en la presión debido a la expansión de los gases, y viceversa. Esta correlación lineal hace que sea crucial monitorear ambas variables simultáneamente. Si solo se monitoreara una de ellas, podrían pasarse por alto problemas en el proceso. Por ejemplo, un aumento anormal en la temperatura podría no detectarse a tiempo si no se considera su efecto sobre la presión, lo que podría llevar a condiciones peligrosas o a la producción de un producto defectuoso.

Para garantizar que el proceso esté bajo control, se monitorean continuamente la temperatura y la presión. Para ello, se toman muestras periódicas (cada hora) de 5 mediciones consecutivas de temperatura y presión. A continuación se muestran los resultados de 20 muestras de cada variable.

Código 
```{r}
#| label: tbl-data-reaccion
#| tbl-cap: "Datos del proceso de producción de polímeros sintéticos"
#| table-caption: left

vec_medias <- c(102, 204)
matriz_var_cov <- matrix(
  c(2, 3, 3, 5),
  nrow = 2, byrow = TRUE
)

library(mvtnorm)
set.seed(135711)

data_reaccion <- data.frame(
  mvtnorm::rmvnorm(
    n = 10000, mean = vec_medias, sigma = matriz_var_cov
  )
) |>
  dplyr::slice_sample(n = 100) |>
  tibble::as_tibble() |> 
  dplyr::mutate(
    Muestra = as.factor(rep(x = 1:20, each = 5)), 
    .before = X1
  )

data_reaccion <- dplyr::bind_cols(
  with(data = data_reaccion, qcc::qcc.groups(X1, Muestra)),
  with(data = data_reaccion, qcc::qcc.groups(X2, Muestra))
) |>
  dplyr::rename(
    Xi11 = "...1", Xi21 = "...2", Xi31 = "...3", 
    Xi41 = "...4", Xi51 = "...5", Xi12 = "...6", 
    Xi22 = "...7", Xi32 = "...8", Xi42 = "...9", 
    Xi52 = "...10"
  ) |>
  dplyr::mutate(
    Muestra = 1:20, .before = Xi11
  )

# La función MASS::mvrnorm() también simula distribuciones normales multivariadas

# gt::gt(data_reaccion) |>
#   gt::cols_label(
#     Muestra = html("<span style='font-size: 14px;'>\\( \\text{Muestra}\\:(i) \\)</span>"),
#     Xi11 = html("<span style='font-size: 14px;'>\\( X_{i11} \\)</span>"),
#     Xi21 = html("<span style='font-size: 14px;'>\\( X_{i21} \\)</span>"),
#     Xi31 = html("<span style='font-size: 14px;'>\\( X_{i31} \\)</span>"),
#     Xi41 = html("<span style='font-size: 14px;'>\\( X_{i41} \\)</span>"),
#     Xi51 = html("<span style='font-size: 14px;'>\\( X_{i51} \\)</span>"),
#     Xi12 = html("<span style='font-size: 14px;'>\\( X_{i12} \\)</span>"),
#     Xi22 = html("<span style='font-size: 14px;'>\\( X_{i22} \\)</span>"),
#     Xi32 = html("<span style='font-size: 14px;'>\\( X_{i32} \\)</span>"),
#     Xi42 = html("<span style='font-size: 14px;'>\\( X_{i42} \\)</span>"),
#     Xi52 = html("<span style='font-size: 14px;'>\\( X_{i52} \\)</span>"),
#     Xi1_barra = html("<span style='font-size: 14px;'>\\( \\overline{X}_{i1} \\)</span>"),
#     Xi2_barra = html("<span style='font-size: 14px;'>\\( \\overline{X}_{i2} \\)</span>"),
#     var_i1 = html("<span style='font-size: 14px;'>\\( S_{i1}^2 \\)</span>"),
#     var_i2 = html("<span style='font-size: 14px;'>\\( S_{i2}^2 \\)</span>"),
#     cov_12 = html("<span style='font-size: 14px;'>\\( S_{i12} \\)</span>")
#   ) |>
#   tab_spanner(
#     label =  html(
#       "<span style='font-size: 14px;'>\\( \\text{Temperatura (ºC)} \\)</span>"
#     ),
#     columns = c("Xi11", "Xi21", "Xi31", "Xi41", "Xi51")
#   ) |>
#   tab_spanner(
#     label =  html(
#       "<span style='font-size: 14px;'>\\( \\text{Presión (kPa)} \\)</span>"
#     ),
#     columns = c("Xi12", "Xi22", "Xi32", "Xi42", "Xi52")
#   ) |>
#   tab_spanner(
#     label =  html(
#       "<span style='font-size: 14px;'>\\( \\text{Medias muestrales} \\)</span>"
#     ),
#     columns = c("Xi1_barra", "Xi2_barra")
#   ) |>
#   tab_spanner(
#     label =  html(
#       "<span style='font-size: 14px;'>\\( \\text{varianzas y covarianzas muestrales} \\)</span>"
#     ),
#     columns = c("var_i1", "var_i2", "cov_12")
#   ) |>
#   fmt_number(decimals = 2, sep_mark = ".", dec_mark = ",")

options(
  reactable.language = reactable::reactableLang(
    searchPlaceholder = "Búsqueda...",
    pageNext = "Siguiente",
    pagePrevious = "Anterior",
    noData = "No existe",
    pageInfo = "{rowStart} de {rowEnd} de {rows} filas"
  )
)

reactable::reactable(
  data_reaccion,
  columns = list(
    Muestra = colDef(
      name = "\\( i \\)",
      sticky = "left",
      align = "right",
      style = list(
        borderRight = "1px solid #eee", fontWeight = "bold"
      ),
      headerStyle = list(borderRight = "1px solid #eee"),
      format = colFormat(
        separators = TRUE, digits = 0, locales = "es-ES"
      )
    ),
    Xi11 = colDef(
      name = "\\( X_{i11} \\)",
      format = colFormat(
        separators = TRUE, digits = 2, locales = "es-ES"
      )
    ),
    Xi21 = colDef(
      name = "\\( X_{i21} \\)",
      format = colFormat(
        separators = TRUE, digits = 2, locales = "es-ES"
      )
    ),
    Xi31 = colDef(
      name = "\\( X_{i31} \\)",
      format = colFormat(
        separators = TRUE, digits = 2, locales = "es-ES"
      )
    ),
    Xi41 = colDef(
      name = "\\( X_{i41} \\)",
      format = colFormat(
        separators = TRUE, digits = 2, locales = "es-ES"
      )
    ),
    Xi51 = colDef(
      name = "\\( X_{i51} \\)",
      format = colFormat(
        separators = TRUE, digits = 2, locales = "es-ES"
      )
    ),
    Xi12 = colDef(
      name = "\\( X_{i12} \\)",
      format = colFormat(
        separators = TRUE, digits = 2, locales = "es-ES"
      )
    ),
    Xi22 = colDef(
      name = "\\( X_{i22} \\)",
      format = colFormat(
        separators = TRUE, digits = 2, locales = "es-ES"
      )
    ),
    Xi32 = colDef(
      name = "\\( X_{i32} \\)",
      format = colFormat(
        separators = TRUE, digits = 2, locales = "es-ES"
      )
    ),
    Xi42 = colDef(
      name = "\\( X_{i42} \\)",
      format = colFormat(
        separators = TRUE, digits = 2, locales = "es-ES"
      )
    ),
    Xi52 = colDef(
      name = "\\( X_{i52} \\)",
      format = colFormat(
        separators = TRUE, digits = 2, locales = "es-ES"
      )
    )
  ),
  columnGroups = list(
    colGroup(
      name = "\\( \\text{Muestra} \\)", 
      columns = c("Muestra")
    ),
    colGroup(
      name = "\\( \\text{Temperatura}\\left( {^\\circ}\\text{C} \\right) \\)", 
      columns = c("Xi11", "Xi21", "Xi31", "Xi41", "Xi51")
    ),
    colGroup(
      name = "\\( \\text{Presión (kPa)} \\)", 
      columns = c("Xi12", "Xi22", "Xi32", "Xi42", "Xi52")
    )
  ),
  defaultColDef = colDef(
    minWidth = 100,
  ),
  fullWidth = TRUE,
  # width = "100%",
  # style = list(width = "100%"),
  resizable = TRUE, 
  wrap = FALSE, 
  bordered = TRUE,
  defaultPageSize = 10,
  striped = TRUE,
  highlight = TRUE,
  filterable = TRUE
)
```
Tabla 5.14: Datos del proceso de producción de polímeros sintéticos

Dado que en este proceso de producción se están monitoreando de manera simultanea dos características relacionadas con la calidad del producto (\(p = 2\)). Entonces, el estadístico \(T_{i}^{2}\) dado en la ecuación 5.43 se puede reescribir de la siguiente manera

\[ \begin{align*} T_{i}^{2} & = 5\left[ \begin{matrix} \bar{x}_{i1} - \bar{\bar{x}}_{1}\\ \bar{x}_{i2} - \bar{\bar{x}}_{2} \end{matrix} \right]^{\pmb{`}} \begin{bmatrix} \bar{S}_{1}^{2} & \bar{S}_{12}\\ \bar{S}_{12} & \bar{S}_{2}^{2} \end{bmatrix}^{-1} \begin{bmatrix} \bar{x}_{i1} - \bar{\bar{x}}_{1}\\ \bar{x}_{i2} - \bar{\bar{x}}_{2} \end{bmatrix}\\ & = 5\left[ \begin{matrix} \bar{x}_{i1} - 101,87\\ \bar{x}_{i2} - 203,75 \end{matrix} \right]^{\pmb{`}} \begin{bmatrix} 2,26 & 3,45\\ 3,45 & 5,71 \end{bmatrix}^{-1} \begin{bmatrix} \bar{x}_{i1} - 101,87\\ \bar{x}_{i2} - 203,75 \end{bmatrix}\\ & = 5\left[ \begin{matrix} \bar{x}_{i1} - 101,87 & \bar{x}_{i2} - 203,75 \end{matrix} \right] \begin{bmatrix} 5,73 & -3,46\\ -3,46 & 2,27 \end{bmatrix} \begin{bmatrix} \bar{x}_{i1} - 101,87\\ \bar{x}_{i2} - 203,75 \end{bmatrix}\\ \end{align*} \]

o de manera equivalente,

\[ \begin{align*} T_{i}^{2} & = \frac{n}{\bar{S}_{1}^{2}\bar{S}_{2}^{2} - \bar{S}_{12}^{2}}\left[\bar{S}_{2}^{2}\left( {\bar{x}}_{i1} - \bar{\bar{x}}_{1} \right)^2 + \bar{S}_{1}^{2}\left( {\bar{x}}_{i2} - \bar{\bar{x}}_{2} \right)^2 \\ -2\,\bar{S}_{12}\left( {\bar{x}}_{i1} - \bar{\bar{x}}_{1} \right)\left( {\bar{x}}_{i2} - \bar{\bar{x}}_{2} \right) \right]\\ & = \frac{5*5,71}{2,26*5,71 - 3,45^{2}}\left( {\bar{x}}_{i1} - 101,87 \right)^2 + \frac{5*2,26}{2,26*5,71 - 3,45^{2}}\left( {\bar{x}}_{i2} - 203,75 \right)^2\\ &\quad \: - \frac{2*5*3,45}{2,26*5,71 - 3,45^{2}}\left( {\bar{x}}_{i1} - 101,87 \right)\left( {\bar{x}}_{i2} - 203,75 \right)\\ & = 28,65\left( {\bar{x}}_{i1} - 101,87 \right)^2 + 11,33\left( {\bar{x}}_{i2} - 203,75 \right)^2 - 34,61\left( {\bar{x}}_{i1} - 101,87 \right)\left( {\bar{x}}_{i2} - 203,75 \right)\\ \end{align*} \]

A continuación, se muestra para la primera muestra, el cálculo de la media, la varianza, la covarianza, el estadístico \(T^{2}\) de Hotelling y el determínate de la matriz de varianzas-covarianzas. Los valores de estos estadísticos para las 20 muestran se pueden observar en la Tabla 5.15.

  • \[ \begin{align*} \overline{x}_{11} &= \frac{x_{111} + x_{121}+ x_{131} + x_{141} +x_{151}}{5} \\ & = \frac{105,01 + 103,56 + 100,77 + 104,55 + 100,15}{5}\\ & = \frac{514,03}{5} = 102,81 \end{align*} \]

  • \[ \begin{align*} \overline{x}_{12} &= \frac{x_{122} + x_{122}+ x_{132} + x_{142} +x_{152}}{5} \\ & = \frac{207,16 + 207,15 + 202,45 + 207,29 + 201,48}{5}\\ & = \frac{1.025,53}{5} = 205,11 \end{align*} \]

  • \[ \begin{align*} S_{1}^{2} & = \frac{\left ( x_{111} - \overline{x}_{11} \right )^2 + \left ( x_{121} - \overline{x}_{11} \right )^2 + \left ( x_{131} - \overline{x}_{11} \right )^2 + \left ( x_{141} - \overline{x}_{11} \right )^2 + \left ( x_{151} - \overline{x}_{11} \right )^2}{5 - 1} \\ & = \frac{\left (105,01 - 102,81 \right )^2 + \left (103,56 - 102,81 \right )^2 + \cdots + \left (100,15 - 102,81 \right )^2}{5 - 1} \\ & = \frac{19,65}{4} = 4,91 \end{align*} \]

  • \[ \begin{align*} S_{2}^{2} & = \frac{\left ( x_{112} - \overline{x}_{12} \right )^2 + \left ( x_{122} - \overline{x}_{12} \right )^2 + \left ( x_{132} - \overline{x}_{12} \right )^2 + \left ( x_{142} - \overline{x}_{12} \right )^2 + \left ( x_{152} - \overline{x}_{12} \right )^2}{5 - 1} \\ & = \frac{\left (207,16 - 205,11 \right )^2 + \left (207,15 - 205,11 \right )^2 + \cdots + \left (201,48 - 205,11 \right )^2}{5 - 1} \\ & = \frac{33,38}{4} = 8,35 \end{align*} \]

  • \[ \begin{align*} S_{112} & = \frac{\left (x_{111} - \overline{x}_{11} \right )\left ( x_{112} - \overline{x}_{12} \right) + \cdots + \left (x_{151} - \overline{x}_{11} \right)\left (x_{152} - \overline{x}_{12} \right)}{5 - 1}\\ & = \frac{\left (105,01 - 100,15 \right )\left ( 207,16 - 205,11 \right) + \cdots + \left (100,15 - 100,15 \right)\left (201,48 - 205,11 \right)}{5 - 1}\\ & = \frac{24,9}{4} = 6,23 \end{align*} \]

  • \[ \begin{align*} T_{1}^{2} & = 5\begin{bmatrix} \bar{x}_{11} - \bar{\bar{x}}_{1} \\ \bar{x}_{12} - \bar{\bar{x}}_{2} \end{bmatrix}^{\pmb{`}}\begin{bmatrix} \bar{S}_{1}^{2} & \bar{S}_{12}\\ \bar{S}_{12} & \bar{S}_{2}^{2} \end{bmatrix}^{-1}\begin{bmatrix} \bar{x}_{11} - \bar{\bar{x}}_{1} \\ \bar{x}_{12} - \bar{\bar{x}}_{2} \end{bmatrix}\\ & = 5\left[ \begin{matrix} 102,81 - 101,87\\ 205,11 - 203,75\\ \end{matrix} \right]^{\pmb{`}}\left[ \begin{matrix} 2,26 & 3,45\\ 3,45 & 5,71 \end{matrix} \right]^{-1}\left[ \begin{matrix} 102,81 - 101,87\\ 205,11 - 203,75 \end{matrix} \right]\\ & = \left[ \begin{matrix} 4,68 & 6,78\\ \end{matrix} \right]\left[ \begin{matrix} 5,73 & -3,46\\ -3,46 & 2,27 \end{matrix} \right]\left[ \begin{matrix} 0,94\\ 1,36 \end{matrix} \right] = \left[ \begin{matrix} 3,34 & -0,83 \end{matrix} \right]\left[ \begin{matrix} 0,94 \\ 1,36 \end{matrix} \right]\\ & = 2 \end{align*} \]

  • \[ \begin{equation} \left| S_{1} \right| = \left| \begin{matrix} S_{11}^{2} & S_{112}\\ S_{112} & S_{12}^{2} \end{matrix} \right| = \left| \begin{matrix} 4,91 & 6,23\\ 6,23 & 8,35 \end{matrix} \right| = 2,23 \end{equation} \]

El siguiente código muestra el vector de medias y la matriz de varianzas-covarianzas.

Código 
```{r}
p <- 2
n <- 5
m <- 20

q <- qcc::mqcc(
  data = list(
    `Presión` = as.matrix(
      data_reaccion[, 2:6], nrow = m, ncol = p 
    ),
    `Temperatura` = as.matrix(
      data_reaccion[, 7:11], nrow = m, ncol = p 
    )
  ),
  type = "T2",
  plot = FALSE
)

## Vector de medias muestral
Vector_medias <- data_reaccion |>
  dplyr::mutate(
    Xi1_barra = rowMeans(dplyr::across(Xi11:Xi51)),
    Xi2_barra = rowMeans(dplyr::across(Xi12:Xi52))
  ) |>
  dplyr::summarise(
    `Presión` = mean(Xi1_barra),
    `Temperatura` = mean(Xi2_barra)
  ) |> 
  matrix(nrow = p, ncol = 1, byrow = TRUE,
    dimnames = list(
      c("Temperatura", "Presión")
    )
  )

## Matriz de varianzas-covarianzas muestral
matriz_var_cov <- data_reaccion |>
  dplyr::rowwise() |>
  dplyr::transmute(
    var_i1 = var(dplyr::c_across(Xi11:Xi51)),
    cov_i1 = cov(
      x = dplyr::c_across(Xi11:Xi51),
      y = dplyr::c_across(Xi12:Xi52)
    ),
    var_i2 = var(dplyr::c_across(Xi12:Xi52)),
    cov_i2 = cov(
      x = dplyr::c_across(Xi11:Xi51),
      y = dplyr::c_across(Xi12:Xi52)
    )
  ) |> 
  dplyr::ungroup() |> 
  dplyr::summarise(
    mean(var_i1), mean(cov_i1), mean(cov_i2), mean(var_i2)
  ) |>
  as.double() |> 
  matrix(
    nrow = 2, ncol = 2, byrow = TRUE,
    dimnames = list(
      c("Temperatura", "Presión"),
      c("Temperatura", "Presión")
    )
  )

q$center   #Vector de medias usando la función qcc::mqcc
q$cov      #Matriz de varianzas-covarianzas qcc::mqcc
```
#>     Presión Temperatura 
#>    101.8708    203.7496 
#>              Presión Temperatura
#> Presión     2.257963    3.449796
#> Temperatura 3.449796    5.712185

En base a las muestras tomadas del proceso (Tabla 5.14), las medias muestrales (calculadas con la ecuación 5.35), las varianzas muestrales (calculadas con la ecuación 5.36), las covarianzas muestrales (calculadas con la ecuación 5.37), los estadísticos \(T_{i}^2\) (calculados con la ecuación 5.43) y los determinantes de las matrices de varianzas-covarianzas se muestra en la siguiente tabla.

Código 
```{r}
#| label: tbl-data-reaccion-estadisticos
#| tbl-cap: "Estadísticos para los datos del proceso de producción de polímeros sintéticos"
#| table-caption: left

data_reaccion_estad <- data_reaccion |>
  dplyr::rowwise() |>
  dplyr::transmute(
    Muestra = Muestra,
    Xi1_barra = mean(dplyr::c_across(Xi11:Xi51)),
    Xi2_barra = mean(dplyr::c_across(Xi12:Xi52)),
    var_i1 = var(dplyr::c_across(Xi11:Xi51)),
    var_i2 = var(dplyr::c_across(Xi12:Xi52)),
    cov_i12 = cov(
      x = dplyr::c_across(Xi11:Xi51),
      y = dplyr::c_across(Xi12:Xi52)
    ),
    det_mvc = det(
      cov(
        data.frame(
          x = c(c_across(Xi11:Xi51)), 
          y = c(c_across(Xi12:Xi52))
        )
      )
    )
  ) |> 
  dplyr::ungroup() |> 
  dplyr::mutate(
    T2_i = q$statistics, .before = det_mvc
  )

options(
  reactable.language = reactable::reactableLang(
    searchPlaceholder = "Búsqueda...",
    pageNext = "Siguiente",
    pagePrevious = "Anterior",
    noData = "No existe",
    pageInfo = "{rowStart} de {rowEnd} de {rows} filas"
  )
)

reactable::reactable(
  data_reaccion_estad,
  columns = list(
    Muestra = colDef(
      name = "\\( i \\)",
      sticky = "left",
      align = "right",
      style = list(
        borderRight = "1px solid #eee", fontWeight = "bold"
      ),
      headerStyle = list(borderRight = "1px solid #eee"),
      format = colFormat(
        separators = TRUE, digits = 0, locales = "es-ES"
      )
    ),
    Xi1_barra = colDef(
      name = "\\( \\overline{X}_{i1} \\)",
      format = colFormat(
        separators = TRUE, digits = 2, locales = "es-ES"
      )
    ),
    Xi2_barra = colDef(
      name = "\\( \\overline{X}_{i2} \\)",
      format = colFormat(
        separators = TRUE, digits = 2, locales = "es-ES"
      )
    ),
    var_i1 = colDef(
      name = "\\( S_{i1}^2 \\)",
      format = colFormat(
        separators = TRUE, digits = 2, locales = "es-ES"
      )
    ),
    var_i2 = colDef(
      name = "\\( S_{i2}^2 \\)",
      format = colFormat(
        separators = TRUE, digits = 2, locales = "es-ES"
      )
    ),
    cov_i12 = colDef(
      name = "\\( S_{i12} \\)",
      format = colFormat(
        separators = TRUE, digits = 2, locales = "es-ES"
      )
    ),
    det_mvc = colDef(
      name = "\\( \\left| S_{i} \\right| \\)",
      format = colFormat(
        separators = TRUE, digits = 2, locales = "es-ES"
      )
    ),
    T2_i = colDef(
      name = "\\( T_{i}^{2} \\)",
      format = colFormat(
        separators = TRUE, digits = 2, locales = "es-ES"
      )
    )
  ),
  columnGroups = list(
    colGroup(
      name = "\\( \\text{Muestra} \\)", 
      columns = c("Muestra")
    ),
    colGroup(
      name = "\\( \\text{Media muestral} \\)", 
      columns = c("Xi1_barra", "Xi2_barra")
    ),
    colGroup(
      name = "\\( \\text{Varianza y covarianza muestral} \\)", 
      columns = c("var_i1", "var_i2", "cov_i12")
    ),
    colGroup(
      name = "\\( \\text{Estadístco del diagrama de control} \\)", 
      columns = c("T2_i", "det_mvc")
    )
  ),
  defaultColDef = colDef(
    minWidth = 100,
  ),
  fullWidth = TRUE,
  resizable = TRUE, 
  wrap = TRUE, 
  bordered = TRUE,
  defaultPageSize = 10,
  striped = TRUE,
  highlight = TRUE,
  filterable = TRUE
)
```
Tabla 5.15: Estadísticos para los datos del proceso de producción de polímeros sintéticos

El límite de control superior del diagrama de control \(T^2\) de Hotelling en la fase I se muestra a continuación, para \(\alpha^{\prime} = 1 - (1 - \alpha)^p = 1 - (1 - 0,\!0027)^2 = 0,0054\)

\[ \begin{align*} LSC & = \frac{p\left(m - 1\right)\left( n -1\right)}{mn -m -p +1}F_{1-\alpha, \:p, \: mn - m - p + 1} = \frac{2\left(20 - 1\right)\left(5 -1\right)}{20*5 - 20 - 2 + 1}F_{1 - 0,0054, \:2, \: 20*5 - 20- 2 + 1}\\ & = \frac{2\left(19\right)\left(4\right)}{79}F_{0,9946, \:2, \: 79} = 1,9241 * 5,5744 = 10,7434 \end{align*} \] Los límites de control calculados con se muestran a continuación.

Código 
```{r}
#| label: lsc-hotelling

LIC <- 0
LSC <- p * (m - 1) * (n - 1) / (m * n - m - p + 1) * qf((1 - 0.0027)^2, p, m * n - m - p + 1)

# Límites de control con el paquete qcc
paste("LSC =", formato(LSC), "LIC =", LIC)
q$limits # q es un objeto qcc::mqcc()
```
#> [1] "LSC = 10,7434 LIC = 0"
#>  LCL      UCL
#>    0 10.74335

A continuación se muestran los diagramas de control \(T^2\) de Hotelling con los paquetes ggplot2, qcc, qcr, highcharter y plotly. Como se puede observar en estos gráficos, el proceso de producción de polímeros sintéticos está bajo control, ya que todos los estadístico \(T^2\) de Hotelling graficados se encuentran dentro de los límites de control y no muestran un comportamiento no aleatorio o sistemático.

Código 
```{r}
#| label: fig-reaccion-ggplot2
#| fig-cap: "Diagrama $T^2$ de Hotelling para los datos de @tbl-data-reaccion-estadisticos con el paquete `ggplot2`"


data_reaccion_estad <- data_reaccion_estad |> 
  dplyr::mutate(
    violacion = ifelse(
      T2_i <= LSC, "No", "Si"
    ),
    color = ifelse(
      T2_i <= LSC, "blue", "red"
    ),
    LIC = LIC,
    LSC = LSC
  )


ggplot(
  data = data_reaccion_estad, 
  mapping = aes(x = Muestra, y = T2_i),
  show.legend = FALSE
) +
  geom_step(
    mapping = aes(x = Muestra, y = LIC), 
    direction = "hv", color = "red", lty = "dashed"
  ) +
  geom_step(
    mapping = aes(x = Muestra, y = LSC),
    direction = "hv", color = "red", lty = "dashed"
  ) +
  geom_point(
    mapping = aes(
      x = Muestra, y = T2_i, color = factor(violacion), 
      shape = factor(violacion)
    ),
    size = 2
  ) +
  scale_color_manual(
    values = c("No" = "blue", "Si" = "red")
  ) +
  scale_shape_manual(values = c(16, 8)) + 
  geom_line(
    mapping = aes(x = Muestra, y = T2_i, group = 1), 
    color = "black"
  ) +
  pammtools::geom_stepribbon(
    mapping = aes(ymin = LIC, ymax = LSC, group = 1), 
    fill = "green", alpha = 0.2
  ) +
  pammtools::geom_stepribbon(
    mapping = aes(ymin = LSC, ymax = Inf, group = 1), 
    fill = "red", alpha = 0.2
  ) +
  pammtools::geom_stepribbon(
    mapping = aes(ymin = LIC, ymax = -Inf, group = 1), 
    fill = "red", alpha = 0.2
  ) +
  annotate(
    "text",
    x = nrow(data_reaccion_estad) + 1, y = LIC, 
    label = "LIC", color = "red"
  ) +
  annotate(
    "text",
    x = nrow(data_reaccion_estad) + 1, y = LSC, 
    label = "LSC", color = "red"
  ) +
  labs(
    x = latex2exp::TeX(
      "$Muestra \\, \\left(i\\right)$", 
      bold = TRUE, italic = TRUE
    ),
    y = latex2exp::TeX(
      "$T^{2}_{i}$",
      bold = TRUE, italic = TRUE
    ),
    shape = "Fuera de control", color = "Fuera de control"
  ) +
  scale_x_continuous(
    breaks = seq(
      from = 1, to = nrow(data_reaccion_estad), by = 2
    ),
    limits = c(1, nrow(data_reaccion_estad) + 1)
  ) +
  theme_bw() +
  theme(
    legend.position = "bottom",
    legend.text = element_text(face = "bold"),
    legend.title.position = "top"
  ) 
```
Figura 5.25: Diagrama \(T^2\) de Hotelling para los datos de Tabla 5.15 con el paquete ggplot2
Código 
```{r}
#| label: fig-reaccion-qcc
#| fig-cap: "Diagrama $T^2$ de Hotelling para los datos de @tbl-data-reaccion-estadisticos con el paquete `qcc`"

# Restablecer parámetros gráficos a los valores por defecto
op <- par(no.readonly = TRUE)  # Guardar los parámetros actuales
on.exit(par(op))                # Asegurarse de que se restablezcan al salir
# Ajustar márgenes
par(mar = c(5, 6, 4, 2) + 0.1)
# Crear el gráfico
plot(
  x = q, add.stats = TRUE, label.limits = c("LIC", "LSC"),
  title = "", xlab = "Periodo", 
  ylab = expression(T[i]^2),
  axes.las = 0, 
  digits = getOption("digits"), restore.par = TRUE
)
```
Figura 5.26: Diagrama \(T^2\) de Hotelling para los datos de Tabla 5.15 con el paquete qcc
Código 
```{r}
#| label: fig-reaccion-qcr
#| fig-cap: "Diagrama $T^2$ de Hotelling para los datos de @tbl-data-reaccion-estadisticos con el paquete `qcr`"

# Función para transformar la tibble
convert_to_array <- function(tibble_data, m, p, n) {
  # Convertir la tibble a matriz para facilitar el procesamiento
  data_matrix <- as.matrix(tibble_data)

  # Verificar que las dimensiones sean compatibles
  if (ncol(data_matrix) != p * n) {
    stop(
      "Las dimensiones de las columnas no coinciden con p * n"
    )
  }

  # Crear el arreglo con la dimensión deseada
  array_data <- array(data = NA, dim = c(m, p, n))

  # Llenar el arreglo con los datos de la tibble
  for (i in 1:n) {
    array_data[, , i] <- data_matrix[
      1:m, i - 1 + seq(from = 1, to = n * p, by = n)
    ]
  }

  return(array_data)
}

data.mqcd <- mqcd(
  convert_to_array(data_reaccion[-1], m = 20, p = 2, n = 5)
)

res.mqcs <- mqcs.t2(
  x = data.mqcd, alpha = 1 - (1 - 0.0027)^2,
  phase = 1 ,plot = FALSE
)

# Crear el gráfico
plot(
  res.mqcs, title = "", xlab = "Periodo", 
  axes.las = 0, 
  digits = getOption("digits"), restore.par = TRUE
)

# Estadísticos T^2 de Hotelling
# res.mqcs$statistics
# Violación de las reglas
# límites de control
# res.mqcs$limits
# res.mqcs$violations
# Resumen estadístico
# summary(res.mqcs)
```
Figura 5.27: Diagrama \(T^2\) de Hotelling para los datos de Tabla 5.15 con el paquete qcr
Código 
```{r}
#| label: fig-reaccion-highcharter
#| fig-cap: "Diagrama $T^2$ de Hotelling para los datos de @tbl-data-reaccion-estadisticos con el paquete `highcharter`"

color <- ifelse(
  data_reaccion_estad$Muestra == unique(q$violations),
  "#FF0000", "#0000FF"
)
symbol <- ifelse(
  data_reaccion_estad$Muestra == unique(q$violations),
  "asterisk-open", "circle"
)

highchart() |>
  hc_add_series(
    data = data_reaccion_estad,
    hcaes(x = Muestra, low = LIC, high = LSC),
    step = "hv",
    marker = list(enabled = FALSE, visible = FALSE),
    type = "arearange", 
    name = "Zona de control", 
    color = "green", fillOpacity = 0.4,
    showInLegend = TRUE,
    enableMouseTracking = FALSE #
  ) |>
  hc_add_series(
    data_reaccion_estad,
    type = "line",
    hcaes(x = Muestra, y = LSC),
    useHTML = TRUE, color = "red", dashStyle = "Dash",
    name = "LSC",
    step = "hv", marker = list(enabled = FALSE),
    tooltip = list(
      valueDecimals = 3
    )
  ) |>
  hc_add_series(
    data_reaccion_estad,
    type = "line",
    hcaes(x = Muestra, y = LIC),
    useHTML = TRUE, color = "red", dashStyle = "Dash",
    name = "LIC", step = "hv", 
    marker = list(enabled = FALSE),
    tooltip = list(
      useHTML = TRUE,
      valueDecimals = 0
    )
  ) |>
  hc_add_series(
    data_reaccion_estad,
    type = "line", color = "blue",
    hcaes(
      x = Muestra, y = T2_i, color = violacion
    ),
    marker = list(
      enabled = TRUE, 
      symbol = "circle", fillColor = "white", radius = 3,
      lineWidth = 2, lineColor = "blue"
    ),
    useHTML = TRUE, name = "T²",
    dashStyle = "Dash",
    tooltip = list(
      valueDecimals = 3
    )
  ) |>
  hc_xAxis(
    title = list(
      text = "<span class='mathjax-process'>\\(\\text{Muestra } (i)\\)</span>", 
      useHTML = TRUE
    )
  ) |>
  hc_yAxis(
    title = list(
      text = "<span class='mathjax-process'>\\(T_{i}^{2}\\)</span>", 
      useHTML = TRUE,
      align = "middle", # alineació: "low", "high"
      x = -10, #desplazar el título del eje x a la izquierda
      rotation = 0 # Rotar el título
    )
  ) |>
  hc_plotOptions(series = list(animation = TRUE)) |>
  hc_tooltip(
    crosshairs = TRUE,
    borderWidth = 2,
    sort = TRUE,
    table = TRUE
  ) |>
  hc_legend(
    labelFormat = "{name}",
    useHTML = TRUE
  )
```
Figura 5.28: Diagrama \(T^2\) de Hotelling para los datos de Tabla 5.15 con el paquete highcharter
Código 
```{r}
#| label: fig-reacion-plotly
#| fig-cap: "Diagrama $T^2$ de Hotelling para los datos de @tbl-data-reaccion-estadisticos con el paquete `plotly`"

plot_ly(data_reaccion_estad, x = ~Muestra) |>
  add_lines(
    y = LSC, name = "<b><i>LSC</i></b>", 
    line = list(
      shape = "hvh", color = "red", dash = "dash"
    ),
    hovertemplate = paste(
      "<b><i>Muestra = %{x}</i></b>",
      "<br><b><i>LSC<b><i> = %{y:.4f}<br><extra></extra>"
    )
  ) |>
  add_lines(
    y = LIC, name = "<b><i>LC</i></b>", fill = "tonexty", 
    fillcolor = 'green',
    hovertemplate = paste(
      "<b><i>Muestra = %{x}</i></b>",
      "<br><b><i>LIC<b><i> = %{y:.0f}<br><extra></extra>"
    ),
    line = list(
      shape = "linear", color = "red", dash = "dash"
    )
  ) |>
  add_trace(
    y = ~T2_i,
    name = "<b><i>T<sup>2</sup><sub>i</sub><b><i>", 
    type = "scatter",
    mode = "lines+markers",
    marker = list(color = ~color, symbol = ~symbol),
    hovertemplate = paste(
      "<b><i>Muestra = %{x}</i></b>",
      "<br><b><i>T<sup>2</sup><sub>%{x}</sub><b><i> = %{y:.3f}<br><extra></extra>"
    ),
    line = list(color = "black", dash = "dash")
  ) |>
  layout(
    legend = list(
      x = 0.5, y = -0.23, xref = 'paper', yref = 'paper',
      orientation = "h", xanchor = 'center', showarrow = F
    ),
    xaxis = list(
      title = "<b><i>Muestra (i) </i></b><br>",
      zeroline = F
    ),
    yaxis = list(
      title = "<b><i>T<sup>2</sup></i></b>",
      zeroline = FALSE
    )
  ) |>
  plotly::config(locale = "es", mathjax = "cdn") |>
  fillOpacity(alpha = 0.3)
```
Figura 5.29: Diagrama \(T^2\) de Hotelling para los datos de Tabla 5.15 con el paquete plotly

El límite de control superior apropiado ha ser usado en la fase II, de acuerdo con la ecuación 5.45, sería:

\[ \begin{align*} LSC & = \frac{p\left(m + 1\right)\left( n -1\right)}{mn -m -p +1}F_{1-\alpha, \:p, \: mn - m - p + 1} = \frac{2\left(20 + 1\right)\left(5 -1\right)}{20*5 - 20 - 2 + 1}F_{1 - 0,0054, \:2, \: 20*5 - 20- 2 + 1}\\ & = \frac{2\left(21\right)\left(4\right)}{79}F_{0,9946, \:2, \: 79} = 2,1266 * 5,5744 = 11,8742 \end{align*} \]