Estadística Industrial

require(tidyverse)
#install.packages('gridExtra')
require(gridExtra) #Libreria para hacer paneles de gráficas con ggplot2

1 Introducción

La Estadística Industrial nace en el contexto del área del control de calidad para medir, monitorear y anticiparse a posibles desviaciones en los procesos.

En este campo vamos a distinguir distintos conceptos:

• Defectos: es cualquier no conformidad de un producto respecto a las especificaciones de diseño de un producto, proceso o servicio

• Defectivo: es cualquier producto, proceso o servicio con al menos un defecto

• Límites de especificación: es la variación permitida por el cliente o diseño funcional de un producto, proceso o servicio

• Límites de control: es el rango de valores bajo los cuales se espera que un proceso fluctúe con cierto grado o nivel de probabilidad

En este módulo vamos a aprender:

1. Gráficos de control y su interpretación
2. Análisis de capacidad del proceso
3. Funciones definidas por el usuario
4. Ciclos for, if, while

2 Gráficos de Control para Variables Continuas

La base de un gráfico de control es el análisis de series de tiempo. Una serie temporal o cronológica es una sucesión de datos medidos en determinados momentos y ordenados cronológicamente. Los datos pueden estar espaciados a intervalos iguales (como la temperatura en un observatorio meteorológico en días sucesivos al mediodía) o desiguales (como el peso de una persona en sucesivas mediciones en el consultorio médico, la farmacia, etc.).

Uno de los usos más habituales de las series de datos temporales es su análisis para predicción y pronóstico (así se hace por ejemplo con los datos climáticos, las acciones de bolsa, o las series de datos demográficos). Resulta difícil imaginar una rama de las ciencias en la que no aparezcan datos que puedan ser considerados como series temporales. Las series temporales se estudian en estadística, procesamiento de señales, econometría y muchas otras áreas.

2.1 Componentes de las series de tiempo

El análisis más clásico de las series temporales se basa en que los valores que toma la variable de observación es la consecuencia de cinco componentes, cuya actuación conjunta da como resultado los valores medidos, estos componentes son:

• Tendencia secular o regular, que indica la marcha general y persistente del fenómeno observado, es una componente de la serie que refleja la evolución a largo plazo. Por ejemplo, la tendencia creciente del índice de reciclado de basuras en los países desarrollados, o el uso creciente de Internet en la sociedad, independientemente de que en un mes concreto en un país, por determinadas causas, haya una baja en la utilización de Internet.

• Variación estacional o variación cíclica regular, el movimiento periódico de corto plazo. Se trata de una componente causal debida a la influencia de ciertos fenómenos que se repiten de manera periódica en un año (las estaciones), una semana (los fines de semana) o un día (las horas puntas) o cualquier otro periodo. Recoge las oscilaciones que se producen en esos períodos de repetición.

• Variación cíclica, el componente de la serie que recoge las oscilaciones periódicas de amplitud superior a un año. movimientos normalmente irregulares alrededor de la tendencia, en las que a diferencia de las variaciones estacionales, tiene un período y amplitud variables, pudiendo clasificarse como cíclicos, cuasicíclicos o recurrentes.

• Variación aleatoria o ruido, accidental, de carácter errático, también denominada residuo, no muestran ninguna regularidad (salvo las regularidades estadísticas), debidos a fenómenos de carácter ocasional como pueden ser tormentas, terremotos, inundaciones, huelgas, guerras, avances tecnológicos, etcétera.

• Variación transitoria, accidental, de carácter errático, debido a fenómenos aislados que pueden llegar a modificar el comportamiento de la serie (tendencia, estacionalidad variaciones cíclicas y aleatoria).

Una serie de tiempo \(y_t\) se define como \(y_t = T_t+S_t+\epsilon_t\) en donde \(T_t\) es la tendencia, \(S_t\) es el componente cíclico o estacional regular y \(\epsilon_t\) es la variación aleatoria.

2.2 Proceso Estacionario

Una serie de tiempo se dice que es un proceso estacionario si la media \(\mu\) y la varianza \(\sigma^2\) se mantienen constantes a lo largo del tiempo.

Si \(y_t\approx N(\mu,\sigma^2)\) y el proceso es estacionario, entonces, para cada momento \(t \in T\) la distribución de probabilidades de \(y_t\) es \(N(\mu,\sigma^2)\).

Más aún, se dice que un proceso es estacionario fuerte si todos los momentos se mantienen constantes en para cada momento de tiempo \(t \in T\).

2.2.1 Ruido Blanco

El ruido blanco o sonido blanco es una señal aleatoria (proceso estocástico) que se caracteriza por el hecho de que sus valores de señal en dos tiempos diferentes no guardan correlación estadística. Como consecuencia de ello, su densidad espectral de potencia (PSD, siglas en inglés de power spectral density) es una constante.

#Ejemplo de ruido blanco
set.seed(1234)
x = rnorm(1000, mean = 0, sd = 1)
plot(x, type = 'l')

Una forma de medir si una señal es ruido blanco es medir la correlación temporal de la serie. Por ejemplo \(\text{cor}(x_t, x_{t-1})\):

df = data.frame(x = x, x_lag1 = lag(x))
df = na.omit(df) #Eliminar NA´s
ggplot(df, aes(x = x, y = x_lag1)) +
  geom_point()

#Imprimir la correlacion
print(cor(df$x, df$x_lag1))

## [1] 0.02615475

En este caso, el índice de correlación es muy cercano a \(0\) por lo que podríamos inferir que esta serie de ruido blanco es estacionaria. Sin embargo, también tendríamos que obtener \(\text{cor}(x_t, x_{t-2})\)

df = data.frame(x = x, x_lag2 = lag(x, n = 2))
df = na.omit(df) #Eliminar NA´s
ggplot(df, aes(x = x, y = x_lag2)) +
  geom_point()

#Imprimir la correlacion
print(cor(df$x, df$x_lag2))

## [1] -0.03030294

y también vemos que el coeficiente de correlaciones es cercano a \(0\).

2.2.2 Función de Autocorrelación

La autocorrelación o dependencia secuencial es una herramienta estadística utilizada frecuentemente en el procesado de señales.

La función de autocorrelación se define como la correlación cruzada de la señal consigo misma. La función de autocorrelación resulta de gran utilidad para encontrar patrones repetitivos dentro de una señal, como la periodicidad de una señal enmascarada bajo el ruido o para identificar la frecuencia fundamental de una señal que no contiene dicha componente, pero aparecen numerosas frecuencias armónicas de esta.

Los procesos de raíz unitaria, autorregresivos, de tendencia estacionaria y los Modelos de Medias Móviles; son ejemplos de procesos con autocorrelación.

En estadística, la autocorrelación de una serie temporal discreta de un proceso X_t no es más que simplemente la correlación de dicho proceso con una versión desplazada en el tiempo de la propia serie temporal.

Si X_t representa un proceso estacionario de segundo orden con un valor principal de μ se define entonces:

\(R(k) = \frac{E[X_i-\mu]E[X_{i-k}-\mu]}{\sigma^2}\)

donde \(E\) es el valor esperado y \(k\) el desplazamiento temporal considerado (normalmente denominado desfase). Esta función varía dentro del rango [−1, 1], donde 1 indica una correlación perfecta (la señal se superpone perfectamente tras un desplazamiento temporal de \(k\)) y −1 indica una anticorrelación perfecta.

Es una práctica común en muchas disciplinas el abandonar la normalización por σ² y utilizar los términos autocorrelación y autocovarianza de manera intercambiable.

En R podemo calcularla función de autocorrelación con la función acf(x):

acf(df$x)

Esta gráfica lo que mide es la autocorrelación de \(x_t, x_{t-k}\) con \(k=30\) . El proceso es estacionario si para cada periodo \(k>0\) el tamaño de la línea no sobrepasa los límites de confianza (líneas azules).

En este caso es evidente que el proceso \(x_t\) es no autocorrelacionado y hay indicios de que se trata de un proceso estacionario.

2.2.3 Caminata Aleatoria

La caminata aleatoria o paseo aleatorio o camino aleatorio, abreviado en inglés como RW (Random Walks), es una formalización matemática de la trayectoria que resulta de hacer sucesivos pasos aleatorios. Por ejemplo, la ruta trazada por una molécula mientras viaja por un líquido o un gas, el camino que sigue un animal en su búsqueda de comida, el precio de una acción fluctuante y la situación financiera de un jugador pueden tratarse como una caminata aleatoria. El término caminata aleatoria fue introducido por Karl Pearson en 1905.​ Los resultados del estudio de las caminatas aleatorias han sido aplicados a muchos campos como la computación, la física, la química, la ecología, la biología, la psicología o la economía.

Una caminata aleatoria se define como:

\(X(t+\tau) = X(t)+\phi(\tau)\)

Vamos a simular una caminata aleatoria:

set.seed(1234)

n <- 1000
x <- cumsum(rnorm(n))
plot(x, type = 'l')

Ahora vamos a revisar el supuesto de autocorrelación:

acf(x)

Ahora podemos ver que el proceso está autocorrelacionado. Esto quiere decir que quizás el proceso no es estacionario y que su distribución de probabilidades podría no ser constante en el tiempo.

En realidad, para saber si el proceso es estacionario, tendríamos que hacer varias cosas adicionales pero por ahora, supondremos que este proceso es no estacionario aun cuando tal vez si lo sea.

2.2.4 Importancia de la estacionariedad en control de procesos

La importancia de la estacionariedad radica en que si se cumple que un proceso es estacionario, entonces, la media y la varianza son constantes en el tiempo, lo que quiere decir que la distribución de probabilidades es constante en el tiempo (la misma en cada momento) y por lo tanto, con cierto grado de error o incertidumbre podemos predecir la salida del un proceso.

Generalmente, para medir la incertidumbre utilizamos la definición de un intervalo de confianza:

\(\bar{x}-z_{\alpha/2}\frac{s}{\sqrt(n)} \leq \hat{\mu}\leq \bar{x}+z_{\alpha/2}\frac{s}{\sqrt(n)}\)

En donde \(\alpha\) es la significancia que usualmente es \(0.05\) y \(z\) el el valor en la curva normal estándar al cuál corresponde un área bajo la curva de \(\alpha/2\)

Tomando el ejemplo anterior de ruido blanco, el intervalo de confianza es:

set.seed(1234)
x = rnorm(1000, mean = 0, sd = 1)

a = 0.05
z = -qnorm(a/2, mean = 0, sd = 1)
l_inf = mean(x) - z * sd(x) #asume que el tamaño de la muestra es de 1
l_sup = mean(x) + z * sd(x) #asume que el tamaño de la muestra es de 1
print(c(l_inf, l_sup))

## [1] -1.981343  1.928149

#Graficat los resultados
plot(x, type = 'l')
abline(h = l_inf, col = 'red')
abline(h = l_sup, col = 'red')

Esto quiere decir que. 95 de cada 100 observaciones estarán entre \([-1.98 , 1.92]\). Así, si una observación cae por fuera de este intervalo, podríamos pensar en la existencia de una anomalía en el proceso.

2.3 Gráfico IMR o de observaciones individuales

Este tipo de gráfico se utiliza cuando sólo es posible obtener 1 dato de forma secuencial fuera de un grupo racional. Un grupo o subgrupo racional es un conjunto de datos de tamaño \(n\) que han sido recolectados bajo las mismas circunstancias. Por ejemplo:

1. Mismo turno
2. Mismo operador
3. Misma máquina
4. Misma hora del día

La idea general de hacer muestreos con grupos racionales es minimizar el ruido ocasionado por estas fuentes de variación común en los procesos. Cuando no es posible asegurar un grupo racional o las mediciones toman demasiado tiempo, utilizamos el gráfico IMR. La letra I es de individual y las letras MR de moving range.

En este caso, los límites o intervalos de confianza se computan como:

\(LCL = \bar{x_t}-3\frac{\bar{MR}}{1.128}\)

\(UCL = \bar{x_t}+3\frac{\bar{MR}}{1.128}\)

\(MR = |x_t-x_{t-1}|\)

Por lo general, para observar la estacionariedad del proceso, graficamos la media del proceso y la desviación estándar. En el contexto de las gráficas de control, el rango móvil o \(MR\) es una aproximación de la variaición del proceso. Los límites de control del rango móvil son:

\(LCL_s=0\)

\(UCL_s=3.267\bar{MR}\)

Ahora vamos a generar las gráficas de control:

#Algoritmo para una gráfica IMR
#Vamos a utilizar datos simulados:
set.seed(1234)
x = rnorm(100, mean = 20, sd = 5)
xbar = mean(x)
d2 = 1.128

#Paso 1 computar el rango movil MR
MR = c(NA, abs(diff(x)))
MRbar = mean(MR, na.rm = T)


#Paso 2 poner todo en un solo data frame
IMR = data.frame(id = seq(1, length(x)),x = x, MR = MR, 
                 xbar = xbar, LCLx = NA,UCLx = NA, 
                 MRbar = MRbar, LCLs = NA, UCLs = NA)

#Computar los Limites de Control
IMR = IMR %>%
  mutate(LCLx = xbar - 3 * (MRbar / d2),
         UCLx = xbar + 3 * MRbar,
         LCLs = 0,
         UCLs = 3.267 * (MRbar / d2))

#Paso 3 Graficar

g_I = ggplot(IMR, aes(x = id, y = x)) +
  geom_line(col = 'blue4') +
  geom_line(aes(y = LCLx), col = 'red') +
  geom_line(aes(y = UCLx), col = 'red') +
  geom_line(aes(y = xbar), col = 'grey') +
  theme_bw() +
  labs(title = 'Grafico de observaciones Individuales')

g_MR = ggplot(IMR, aes(x = id, y = MR)) +
  geom_line(col = 'blue4') +
  geom_line(aes(y = LCLs), col = 'red') +
  geom_line(aes(y = UCLs), col = 'red') +
  geom_line(aes(y = MRbar), col = 'grey') +
  theme_bw() +
  labs(title = 'Grafico de rangos móviles')

grid.arrange(g_I, g_MR)

## Warning: Removed 1 row(s) containing missing values (geom_path).

Ahora vamos a crear una función que haga todas estas tareas:

IMR = function(data, chart = TRUE){
  x =data
  xbar = mean(x)
  d2 = 1.128
  
  #Paso 1 computar el rango movil MR
  MR = c(NA, abs(diff(x)))
  MRbar = mean(MR, na.rm = T)
  
  
  #Paso 2 poner todo en un solo data frame
  IMR = data.frame(id = seq(1, length(x)),x = x, MR = MR, 
                   xbar = xbar, LCLx = NA,UCLx = NA, 
                   MRbar = MRbar, LCLs = NA, UCLs = NA)
  
  #Computar los Limites de Control
  IMR = IMR %>%
    mutate(LCLx = xbar - 3 * (MRbar / d2),
           UCLx = xbar + 3 * MRbar,
           LCLs = 0,
           UCLs = 3.267 * (MRbar / d2))
  
  #Computar las señales
  #Señal 1: puntos fuera de los límites de control
  IMR$S1_x = ifelse(IMR$x > IMR$UCLx | IMR$x < IMR$LCLx, 'S', 'NS')
  IMR$S1_s = ifelse(IMR$x > IMR$UCLs | IMR$x < IMR$LCLs, 'S', 'NS')
  
  
    #Paso 3 Graficar solo si chart == True
  if (chart == TRUE){
    
    g_I = ggplot(IMR, aes(x = id, y = x)) +
      geom_line(col = 'blue4') +
      geom_point(aes(col = S1_x)) +
      geom_line(aes(y = LCLx), col = 'red') +
      geom_line(aes(y = UCLx), col = 'red') +
      geom_line(aes(y = xbar), col = 'grey') +
      theme_bw() +
      labs(title = 'Grafico de observaciones Individuales')
    
    g_MR = ggplot(IMR, aes(x = id, y = MR)) +
      geom_line(col = 'blue4') +
      geom_point(aes(col = S1_s)) +
      geom_line(aes(y = LCLs), col = 'red') +
      geom_line(aes(y = UCLs), col = 'red') +
      geom_line(aes(y = MRbar), col = 'grey') +
      theme_bw() +
      labs(title = 'Grafico de rangos móviles')
    
    p3 = grid.arrange(g_I, g_MR)
  }
  
  return(IMR)
      
}

#correr la función
control_chrt = IMR(x, chart = TRUE)

## Warning: Removed 1 row(s) containing missing values (geom_path).

## Warning: Removed 1 rows containing missing values (geom_point).

head(control_chrt)

2.4 Gráfico XR o de rangos y medias de tamaño \(n\)

Cuando podemos hacer un muestreo con grupos racionales es recomendable tomar muestras pequeñas de tamaño \(n\) en forma consecutiva. Entonces, para cada grupo racional computamos el promedio y su rango móvil y lo que graficamos es la media de cada muestra y la aproximación a la desviación estándar de cada grupo racional con el rango móvil

En este caso, la línea central es:

\(\bar{\bar{x}} = \frac{1}{m}\sum_{j=1}^m\bar{x_1}+\bar{x_2}+...+\bar{x_j}\)

\(\bar{x_j} = \frac{1}{n}\sum_{i=1}^nx_i\)

\(\bar{R}=\frac{1}{m}\sum_{j=1}^mR_j\)

Para la gráfica de medias:

\(LCL_x=\bar{\bar{x}} - A_2 \bar{R}\)

\(UCL_x=\bar{\bar{x}} + A_2 \bar{R}\)

Para la gráfica de Rangos

\(LCL_R=D_3\bar{R}\)

\(UCL_R=D_4\bar{R}\)

Las constantes \(A_2,D_3,D_4\) provienen de las tablas que pueden consultarse aquí:

controlChartsValues = data.frame(
  SampleSize = c(2L,3L,4L,5L,6L,7L,8L,9L,10L,11L,
                 12L,13L,14L,15L,16L,17L,18L,19L,20L,21L,22L,23L,24L,
                 25L),
           A = c(2.12132034,1.73205081,1.5,1.34164079,
                 1.22474487,1.13389342,1.06066017,1,0.9486833,0.90453403,
                 0.8660254,0.83205029,0.80178373,0.77459667,0.75,0.72760688,
                 0.70710678,0.6882472,0.67082039,0.65465367,0.63960215,
                 0.62554324,0.61237244,0.6),
          A2 = c(1.87996982,1.02332595,0.72859745,
                 0.57681907,0.48324654,0.41928346,0.37252746,0.33669694,0.30826327,
                 0.28508386,0.26577752,0.24941705,0.23535081,0.22310904,
                 0.21234548,0.20279577,0.1942569,0.18656944,0.17960626,
                 0.17326489,0.1674623,0.1621284,0.15720601,0.15264729),
          A3 = c(2.65869597,1.95440327,1.62809881,
                 1.42729262,1.28713217,1.18191461,1.09909554,1.0316617,0.97534935,
                 0.9273943,0.88590511,0.84954237,0.81733766,0.78853802,
                 0.76259812,0.73905483,0.71757622,0.69787084,0.67970413,
                 0.66288672,0.64725923,0.63269267,0.61906452,0.60628107),
          c4 = c(0.79788,0.88623,0.92132,0.93999,0.95153,
                 0.95937,0.96503,0.96931,0.97266,0.97535,0.97756,0.97941,
                 0.98097,0.98232,0.98348,0.98451,0.98541,0.98621,0.98693,
                 0.98758,0.98817,0.9887,0.98919,0.98964),
          B3 = c(0,0,0,0,0.03033238,0.11769902,
                 0.18508399,0.23912371,0.28371547,0.32127905,0.35352149,0.38162379,
                 0.40622232,0.42826302,0.44782765,0.46573899,0.48184817,
                 0.49656133,0.51014902,0.52272128,0.53440456,0.54513703,
                 0.555274,0.5647777),
          B4 = c(3.26656757,2.56814425,2.26602646,
                 2.08895411,1.96966762,1.88230098,1.81491601,1.76087629,1.71628453,
                 1.67872095,1.64647851,1.61837621,1.59377768,1.57173698,
                 1.55217235,1.53426101,1.51815183,1.50343867,1.48985098,
                 1.47727872,1.46559544,1.45486297,1.444726,1.4352223),
          B5 = c(0,0,0,0,0.02886217,0.11291691,
                 0.17861161,0.23178501,0.27595869,0.31335952,0.34558847,0.37376615,
                 0.3984919,0.42069133,0.44042954,0.4585247,0.474818,
                 0.48971375,0.50348137,0.51622908,0.52808255,0.53897699,
                 0.54927149,0.5589266),
          B6 = c(2.60632893,2.27596648,2.0877355,
                 1.96359598,1.87419783,1.80582309,1.75144839,1.70683499,1.66936131,
                 1.63734048,1.60953153,1.58505385,1.5634481,1.54394867,
                 1.52653046,1.5104953,1.496002,1.48270625,1.47037863,
                 1.45893092,1.44825745,1.43842301,1.42910851,1.4203534),
          d2 = c(1.12838,1.69257,2.05875,2.32593,2.53441,
                 2.70436,2.8472,2.97003,3.07751,3.17287,3.25846,3.33598,
                 3.40676,3.47183,3.53198,3.58788,3.64006,3.68896,3.73495,
                 3.77834,3.81938,3.85832,3.89535,3.93063),
          d3 = c(0.8525,0.88837,0.87981,0.86408,0.84804,
                 0.8332,0.81983,0.80783,0.79705,0.78731,0.77847,0.77041,
                 0.76302,0.7562,0.7499,0.74404,0.73858,0.73347,0.72868,
                 0.72416,0.7199,0.71588,0.71206,0.70843),
          D1 = c(0,0,0,0,0,0.20476,0.38771,0.54654,
                 0.68636,0.81094,0.92305,1.02475,1.1177,1.20323,1.28228,
                 1.35576,1.42432,1.48855,1.54891,1.60586,1.65968,1.71068,
                 1.75917,1.80534),
          D2 = c(3.68588,4.35768,4.69818,4.91817,5.07853,
                 5.20396,5.30669,5.39352,5.46866,5.5348,5.59387,5.64721,
                 5.69582,5.74043,5.78168,5.82,5.8558,5.88937,5.92099,
                 5.95082,5.97908,6.00596,6.03153,6.05592),
          D3 = c(0,0,0,0,0,0.07571477,0.13617238,
                 0.18401834,0.22302446,0.25558564,0.28327799,0.3071811,
                 0.32808299,0.34656939,0.36304849,0.37787217,0.39129025,0.40351481,
                 0.41470702,0.42501734,0.43454173,0.44337432,0.45160769,
                 0.45930042),
          D4 = c(3.26652369,2.57459367,2.28205464,
                 2.11449614,2.00383127,1.92428523,1.86382762,1.81598166,1.77697554,
                 1.74441436,1.71672201,1.6928189,1.67191701,1.65343061,
                 1.63695151,1.62212783,1.60870975,1.59648519,1.58529298,
                 1.57498266,1.56545827,1.55662568,1.54839231,1.54069958)
)

head(controlChartsValues)

Vamos a hacer una función para agrupar los datos en grupos de tamaño \(n\) :

crea_grupos = function(datos, size){
  #Datos = dataframe con el id del dato y x=observaciones
  #size = tamaño de los subgrupos
  
  n = size
  #Sacamos cuantos grupos de tamaño n pueden existir en los m datos
  k = m / n
  #Creamos la lista vacía
  grupo = c()
  #Iteramos para crear los grupos
  for(i in 1 : k){
    temp = rep(i, n)
    grupo = c(grupo, temp)
  
  }

  datos$grupo = grupo
  
  return(datos)
}

Ahora vamos a probar la función:

#Generamos los datos de forma aleatoria
set.seed(1234)
n = 5
m = 1000
x = rnorm(m, mean = 20, sd = 5)
XR_df = data.frame(id = seq(1, m), x = x)

#Computamos los subgrupos
XR_df = crea_grupos(XR_df, size = 5)
head(XR_df)

El siguiente paso es crear una función que pueda usarse para computar los límites de control y mosrar los gráficos:

control_chrtXR = function(datos, size, chart = TRUE){
  
  XR_df = datos
  n = size
  XR_df = crea_grupos(XR_df, size = n)

  #Agrupamos y sacamos el rango y el promedio por grupo
  XR_df = XR_df %>%
    group_by(grupo) %>%
    summarise(xbar = mean(x, na.rm = T),
              xmin = min(x),
              xmax = max(x)) %>%
    mutate(R = xmax - xmin)
  
  #Computamos el Rango promedio
  XR_df = XR_df %>%
    mutate(Rbar = mean(R, na.rm = T))
  
  #Sacamos el valor de las tablas para computar los ímites de control
  valores = controlChartsValues %>%
    filter(SampleSize == n)
  
  A2 = valores$A2
  D3 = valores$D3
  D4 = valores$D4
  
  #computamos los límites de control
  XR_df = XR_df %>%
    mutate(LC = mean(xbar, na.rm = T)) %>%
    mutate(LCLx = LC - A2 * Rbar,
           UCLx = LC + A2 * Rbar,
           LCLr = D3 * Rbar,
           UCLr = D4 * Rbar)
  
  #Graficamos la carta de control
  
  if(chart == TRUE){
    g_X = ggplot(XR_df, aes(x = grupo, y = xbar)) +
      geom_line(col = 'blue4') +
      geom_point() +
      geom_line(aes(y = LCLx), col = 'red') +
      geom_line(aes(y = UCLx), col = 'red') +
      geom_line(aes(y = LC), col = 'grey') +
      theme_bw() +
      labs(title = 'Grafico X')
    
    g_MR = ggplot(XR_df, aes(x = grupo, y = R)) +
      geom_line(col = 'blue4') +
      geom_point() +
      geom_line(aes(y = LCLr), col = 'red') +
      geom_line(aes(y = UCLr), col = 'red') +
      geom_line(aes(y = Rbar), col = 'grey') +
      theme_bw() +
      labs(title = 'Grafico de rangos móviles')
    
    grid.arrange(g_X, g_MR)
    
  }
  
  
  
  return(XR_df)
  
}

set.seed(1234)
n = 5
m = 1000
x = rnorm(m, mean = 20, sd = 5)
XR_df = data.frame(id = seq(1, m), x = x)

control_chrtXR(XR_df, size = 10, chart = TRUE)

2.5 Gráfico XS o de desviación estándar y medias

La gráfica XR es particularmente útil cuando queremos implementar directamente en el proceso y típicamente en papel y lápiz ya que hacer promedios y calcular rangos es una tarea relativamente fácil. Cuando el tamaño del grupo es de \(n \leq9\) la aproximación de la desviación estándar con rangos es un buen proxy.

Cuando el tamaño del grupo \(n>9\) es recomendable computar la desviación estándar para cada grupo.

La carta de control está definida como:

\(\bar{\bar{x}} = \frac{1}{m}\sum_{j=1}^m\bar{x_1}+\bar{x_2}+...+\bar{x_j}\)

\(\bar{x_j} = \frac{1}{n}\sum_{i=1}^nx_i\)

\(\bar{s}=\frac{1}{m}\sum_{j=1}^ms_j\)

Los límites de control para la gráfica de medias son:

\(LCL_x=\bar{\bar{x}} - \frac{3 \bar{s}}{c_4 \sqrt{n}}\)

\(LCL_x=\bar{\bar{x}} + \frac{3 \bar{s}}{c_4 \sqrt{n}}\)

Los límites de control para la gráfica de desviaciones estándar son:

\(LCL_s=B_3\bar{s}\)

\(UCL_s=B_4\bar{s}\)

La implementación de una función en R es:

control_chrtXS = function(datos, size, chart = TRUE){
  
  XR_df = datos
  n = size
  XR_df = crea_grupos(XR_df, size = n)

  #Agrupamos y sacamos el rango y el promedio por grupo
  XR_df = XR_df %>%
    group_by(grupo) %>%
    summarise(xbar = mean(x, na.rm = T),
              s = sd(x, na.rm = T))

  
  #Computamos la desviación promedio
  XR_df = XR_df %>%
    mutate(sbar = mean(s, na.rm = T))
  
  #Sacamos el valor de las tablas para computar los ímites de control
  valores = controlChartsValues %>%
    filter(SampleSize == n)
  
 c4 = valores$c4
 B3 = valores$B3
 B4 = valores$B4
  
  #computamos los límites de control
  XR_df = XR_df %>%
    mutate(LC = mean(xbar, na.rm = T)) %>%
    mutate(LCLx = LC - (3 * sbar / (c4*sqrt(n))),
           UCLx = LC + (3 * sbar / (c4*sqrt(n))),
           LCLs = B3 * sbar,
           UCLs = B4 * sbar)
  
  #Graficamos la carta de control
  
  if(chart == TRUE){
    g_X = ggplot(XR_df, aes(x = grupo, y = xbar)) +
      geom_line(col = 'blue4') +
      geom_point() +
      geom_line(aes(y = LCLx), col = 'red') +
      geom_line(aes(y = UCLx), col = 'red') +
      geom_line(aes(y = LC), col = 'grey') +
      theme_bw() +
      labs(title = 'Grafico X')
    
    g_MS = ggplot(XR_df, aes(x = grupo, y = s)) +
      geom_line(col = 'blue4') +
      geom_point() +
      geom_line(aes(y = LCLs), col = 'red') +
      geom_line(aes(y = UCLs), col = 'red') +
      geom_line(aes(y = sbar), col = 'grey') +
      theme_bw() +
      labs(title = 'Grafico de desviaciones móviles')
    
    grid.arrange(g_X, g_MS)
    
  }
  
  
  
  return(XR_df)
  
}

Vamos a probar la función

set.seed(1234)
n = 5
m = 1000
x = rnorm(m, mean = 20, sd = 5)
XR_df = data.frame(id = seq(1, m), x = x)

control_chrtXS(XR_df, size = 10, chart = TRUE)

2.6 El Caso de un proceso No Estacionario

Los gráficos de control que hemos visto hasta ahora asumen que \(x_t\) es un proceso estacionario. Existén diversas razones por las cuales un proceso no es estacionario:

1. La media \(\mu\) depende del tiempo y no es constante
2. La varianza \(\sigma^2\) no es constante y depende del tiempo

2.6.1 Gráficos de control cuando \(\mu(t)\)

Cuando el la media del proceso se va moviendo en el tiempo, podemos utilizar una técnica conocida como ARIMA o Autorregresive Integrated Moving Average.

2.6.2 Proceso Autorregresivo AR(1)

El modelo AR(1), es un modelo autorregresivo que está construido únicamente sobre un retardo.

En otras palabras, la autoregresión de primer orden, AR(1), retrocede la autorregresión en un período de tiempo. La forma genérica de representar un AR(1) sería la siguiente:

\(x_t = c + \theta y_{t-1} + \epsilon_t\)

Este modelo asume que el parámetro \(|\theta| < 1\)

En R podemos estimar un AR(1) con la libreria forecast

#Simular un proceso AR(1)
yt <- arima.sim(list(order=c(1,0,0), ar = 0.5), n=500) + 10
plot(yt, type = 'l')

Vamos a estimar el modelo AR(1):

require(forecast)

## Loading required package: forecast

## Warning: package 'forecast' was built under R version 4.0.5

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

ar1 = Arima(yt, order = c(1,0,0))
summary(ar1)

## Series: yt 
## ARIMA(1,0,0) with non-zero mean 
## 
## Coefficients:
##          ar1     mean
##       0.5033  10.0624
## s.e.  0.0387   0.0840
## 
## sigma^2 estimated as 0.8772:  log likelihood=-675.85
## AIC=1357.7   AICc=1357.75   BIC=1370.34
## 
## Training set error measures:
##                        ME      RMSE       MAE        MPE     MAPE      MASE
## Training set -0.001619711 0.9346973 0.7412751 -0.9176503 7.544718 0.8733068
##                    ACF1
## Training set 0.01544549

2.6.3 Proceso Autorregresivo AR(p)

En general, podemos escribir un modelo autorregresivo de orden p o AR(p) como:

\(x_t = c + \sum_{i = 1}^p\theta_p y_{t-p} + \epsilon_t\)

Los valores \(\theta_p\) deben caer dentro de las raíces que definin un círculo unitario.

#Simular un proceso AR(4)
yt <- arima.sim(list(order=c(4,0,0), ar = c(0.7, 0.2, -0.1, -0.3)), n=500) + 10
plot(yt, type = 'l')

vamos a estimar un modelo AR(4):

ar4 = Arima(yt, order = c(4,0,0))
summary(ar4)

## Series: yt 
## ARIMA(4,0,0) with non-zero mean 
## 
## Coefficients:
##          ar1     ar2      ar3      ar4     mean
##       0.6968  0.2002  -0.0273  -0.3323  10.0314
## s.e.  0.0423  0.0528   0.0528   0.0425   0.0998
## 
## sigma^2 estimated as 1.074:  log likelihood=-725.62
## AIC=1463.25   AICc=1463.42   BIC=1488.54
## 
## Training set error measures:
##                       ME     RMSE       MAE      MPE     MAPE      MASE
## Training set 0.004054513 1.031202 0.8276152 -1.11532 8.609222 0.8416459
##                      ACF1
## Training set -0.006712821

2.6.4 Proceso de Media Móvil MA(1)

En análisis de serie del tiempo, el modelo de medias móviles (MA) es una aproximación común para series de tiempo univariadas. El modelo de medias móviles especifica que la variable de salida depende linealmente del valor actual y varios de los anteriores de un término estocástico.

Junto con el modelo autorregresivo (AR), es un caso especial y componente clave de los modelos de serie de tiempo más generales ARMA y ARIMA, los cuales tienen una estructura estocástica más compleja.

Al modelo de medias móviles no se le ha de confundir con la media móvil, un concepto distinto, a pesar de sus semejanzas.

Contrariamente a los modelos autorregresivos, el MA es siempre estacionario.

El modelo de media movil MA(1) se define como:

\(y_t = \mu + \epsilon_t + \gamma_1 \epsilon_{t-1}\)

#Simular un proceso MA(1)
yt <- arima.sim(list(order=c(0,0,1), ma = 0.9), n=500) + 10
plot(yt, type = 'l')

Vamos a estimar el modelo MA(1):

ma1 = Arima(yt, order = c(0,0,1))
summary(ma1)

## Series: yt 
## ARIMA(0,0,1) with non-zero mean 
## 
## Coefficients:
##          ma1     mean
##       0.8856  10.0444
## s.e.  0.0235   0.0855
## 
## sigma^2 estimated as 1.035:  log likelihood=-717.88
## AIC=1441.76   AICc=1441.81   BIC=1454.41
## 
## Training set error measures:
##                       ME     RMSE       MAE       MPE     MAPE     MASE
## Training set 0.001195331 1.015407 0.8066083 -1.133913 8.335987 0.783975
##                    ACF1
## Training set 0.04114502

De la misma manera podemos generalizar el modelo MA(q):

\(y_t = \mu + \epsilon_t + \sum_{j=1}^q\gamma_q \epsilon_{t-q}\)

2.6.5 Modelos ARIMA(p,d,q)

Podemos generalizar un modelo AR(p) y MA(q) para estimar una serie de tiempo:

\(x_t = \mu + \sum_{i=1}^p \theta_p y_{t-p} + \sum_{j=1}^q \gamma_q \epsilon_{t-q} + \epsilon _t\)

Podemos usar una función de la libreria forecast para estimar el orden del modelo ARMA(p,q). Esta función es auto.arima(xt).

#Simular un proceso ARMA(4, 2)
ar_par = c(0.7, 0.2, -0.1, -0.3)
ma_par = c(0.3, -0.6)
yt <- arima.sim(list(order=c(4,1,2), ar = ar_par, ma = ma_par), n=500) + 10
plot(yt, type = 'l')

Estimación del modelo ARMA(p,q)

arima_pq = auto.arima(yt)
summary(arima_pq)

## Series: yt 
## ARIMA(5,1,3) with drift 
## 
## Coefficients:
##          ar1      ar2      ar3      ar4     ar5      ma1      ma2     ma3
##       1.2221  -0.1864  -0.0351  -0.4051  0.1985  -0.2441  -0.8178  0.2011
## s.e.  0.3561   0.3260   0.0816   0.0768  0.1144   0.3590   0.0393  0.2840
##        drift
##       0.0713
## s.e.  0.0302
## 
## sigma^2 estimated as 0.9744:  log likelihood=-700.24
## AIC=1420.47   AICc=1420.92   BIC=1462.62
## 
## Training set error measures:
##                      ME      RMSE       MAE       MPE     MAPE      MASE
## Training set 0.00173744 0.9772325 0.7939182 0.1761182 4.931384 0.5587361
##                     ACF1
## Training set 0.001659735

2.6.6 Gráfico de Control con errores ARIMA

Ya que hemos revisado los conceptos básicos, podemos aprovechar un modelo ARIMA para hacer un gráfico de control de un proceso no estacionario. La teoría dice que los errores de un proceso ARIMA deben cumplir con las siguientes propiedades:

1. Los errores \(\epsilon_t \rightarrow \text{NID}(0, \sigma^2)\)
2. Los errores son no autocorrelacionados

Podemos aprovecharnos de estas 2 propiedades para trabajar con cualquier carta de control XR, XS o IMR. En este caso, vamos a utilizar la gráfica IMR. Ahora tomará sentido haber programado una función para hacer el análisis ya que no tenemos que programar toda la gráfica desde 0.

#Simular un proceso ARMA(3, 1) con drift
ar_par = c(0.7, 0.2, -0.1)
ma_par = c(0.87)
yt <- arima.sim(list(order=c(3,1,1), ar = ar_par, ma = ma_par), n=500) + 10
plot(yt, type = 'l')

Si aplicamos las gráficas de control que hemos revisado, obtenemos las siguientes conclusiones:

#Contruir el data frame. Recordemos que nuestrs funciones piden como input un
#data frame con la columna de id

yt = as.numeric(yt)
IMR(yt)

## Warning: Removed 1 row(s) containing missing values (geom_path).

## Warning: Removed 1 rows containing missing values (geom_point).

De acuerdo a esta gráfica, el proceso está totalmente fuera de control. Cuando aplicamos la carta de control a los errores del modelo ARIMA obtenemos:

arima_pq = auto.arima(yt)
errores = as.numeric(arima_pq$residuals)
IMR(errores)

## Warning: Removed 1 row(s) containing missing values (geom_path).

## Warning: Removed 1 rows containing missing values (geom_point).

ahora lo que observamos es un proceso estacionario centrado en 0 que no tiene puntos fuera de los límites de control y que podemos decir que está en control estadístico.

Si bien, esta técnica para determinar si un proceso está en control, vale la pena preguntarnos:

1. El proceso debería ser estacionario en el tiempo?
2. Si el proceso depende del tiempo siempre debemos preguntarnos porqué
3. Existe alguna variable que sistemáticamente esté afectando el proceso ocasionando cambios en la media? Si es asi, deberíamos eliminar esa causa?

Es importante utilzar este enfoque cuando sabemos que la naturaleza del proceso es no estacionaria. Por ejemplo, sistemas de producción con procesos feedback loop, temporalidades o estacionalidades que se presentan por los cambios de clima o temperatura, mismo por cambios en los patrones de consumo. Sistemas que dependen de organismos o micro organismos tienden a ser no estacionarios (maceración, fermentación, etc). Recuerden que la teoría mata el efuerzo computacional.

Una forma efectiva de lidiar con este problema, es diseñar un plan de grupo racionales que ayude a obtener observaciones con el menor ruido posible.

3 Gráficos de Control para Variables Discretas

Existen 2 tipos de variables discretas que interesan desde el punto de vista del control de calidad:

1. Proporciones
2. Conteos

3.1 Gráficos para Proporciones p

3.1.1 La Distribución Binomial

En estadística, la distribución binomial una distribución de probabilidad discreta que cuenta el número de éxitos en una secuencia \(n\) nensayos de Bernoulli independientes entre sí con una probabilidad fija \(p\) de ocurrencia de éxito entre los ensayos.

Un experimento de Bernoulli se caracteriza por ser dicotómico, esto es, solo dos resultados son posibles, a uno de estos se le denomina éxito y tiene una probabilidad de ocurrencia \(p\) y al otro se le denomina fracaso y tiene una probabilidad \(q = 1-p\).

En general, una variable aleatoria discreta \(X\) tiene una distribución binomial con parámetros \(n\) y \(0\leq p \leq 1\) y escribimos \(x\backsim\text{Binom}(n, p)\) con función de probabilidades:

\(f(x) = {n\choose x}p^x(1-p)^{n-x}\)

La media y la varianza son:

\(\mu = np\)

\(\sigma^2 = np(1-p)\)

3.1.2 Fracción de defectivos p

Si \(D\) es el número de unidades o piezas defectuosas (un defectuoso es un item que al menos tiene un defecto) y \(n\) es el tamaño de la muestra, entonces la fracción de producto no conforme es:

\(\hat{p} = \frac{D}{n}\)

Entonces la media y la varianza son:

\(\mu_p=p\)

\(\sigma_p^2=\frac{p(1-p)}{n}\)

Los límites de control superior e inferior son:

\(LC =\bar{p}\)

\(LCL = \bar{p}-3\sqrt{\frac{\bar{p}(1-\bar{p})}{n}}\)

\(UCL = \bar{p}+3\sqrt{\frac{\bar{p}(1-\bar{p})}{n}}\)

Vamos a simular un conjunto de datos suponiendo que la probabilidad de un defectivo es 10% y que los muestreos en cada ejercicio de inspección podrían variar entre 20 y 50 muestras.

N = 200 #Numer de observaciones
p = 0.20 #Proporcion de defectivos
aleatorio = round(runif(N, min = 20, max = 50),0)
x = data.frame(id = 1 : N, sample_size = aleatorio)
x = x %>%
  mutate(defectives = rbinom(sample_size, sample_size, p))
x

La gráfica de control quedaría de la siguiente forma:

#Proporcion promedio, LC
x$LC = sum(x$defectives / sum(x$sample_size))

#La desviación estádar y la proporcion de defectos es:
x = x %>%
  mutate(p = defectives / sample_size) %>%
  mutate(sigma = sqrt(p * (1 - p) / sample_size))

#Los límites de control Inferior y Superior son:
x = x %>%
  mutate(LCL = LC - 3 * sigma,
         UCL = LC + 3 * sigma)

#Graficamos

g_p = ggplot(x, aes(x = id, y = p)) +
      geom_line(col = 'blue4') +
      geom_point() +
      geom_step(aes(y = LCL), col = 'red') +
      geom_step(aes(y = UCL), col = 'red') +
      geom_line(aes(y = LC), col = 'grey') +
      theme_bw() +
      labs(title = 'Grafico p')
g_p

Un hecho interesante de ls gráfica de control para proporciones es que los límites de control pueden variar con el tamaño de la muestra. La razón es que la varianza, \(p(1-p)/n\) depende del tamaño de la muestra: A mayor tamaño de muestra, menor es la varianza

plot(x$sample_size, x$sigma)

Al igual que antes, podemos hacer una función:

control_chrtp = function(data, chart = TRUE){
  #data es un dataframe con la columna id, el tamaño de la muestra, sample_size y
  #el numero de piezas defectuosas o defectives
  x = data
  
  #Proporcion promedio, LC
  x$LC = sum(x$defectives / sum(x$sample_size))
  
  #La desviación estádar y la proporcion de defectos es:
  x = x %>%
    mutate(p = defectives / sample_size) %>%
    mutate(sigma = sqrt(p * (1 - p) / sample_size))
  
  #Los límites de control Inferior y Superior son:
  x = x %>%
    mutate(LCL = LC - 3 * sigma,
           UCL = LC + 3 * sigma)
  
  #Correginr LCL < 0
  x$LCL = ifelse(x$LCL < 0, 0, x$LCL)
  
  #Graficamos
  
  if(chart == TRUE){
    g_p = ggplot(x, aes(x = id, y = p)) +
          geom_line(col = 'blue4') +
          geom_point() +
          geom_step(aes(y = LCL), col = 'red') +
          geom_step(aes(y = UCL), col = 'red') +
          geom_line(aes(y = LC), col = 'grey') +
          theme_bw() +
          labs(title = 'Grafico p')
    print(g_p)
  }
  
  return(x)
}

Ahora probamos la función:

#Simular observacion
N = 70 #Numer de observaciones
p = 0.20 #Proporcion de defectivos
aleatorio = round(runif(N, min = 20, max = 50),0)
x = data.frame(id = 1 : N, sample_size = aleatorio)
x = x %>%
  mutate(defectives = rbinom(sample_size, sample_size, p))

control_chrtp(x, chart = TRUE)