Control estadístico de procesos
Previsión
Jordi Cuadros, Lucinio González
Noviembre de 2018

Fundamentos de la previsión

Previsión

En la práctica de la gestión (aunque no solo en la industria) a menudo es necesario poder anticipar cuál será el resultado de una determinada situación para poder actuar en previsión.

  • Prever las ventas o las compras, es decir, necesidades de producción del próximo periodo.
  • Prever que pasará ante un cambio de las condiciones de producción.
  • Anticipar la duración de una pieza o evaluar la posibilidades de fallo de un proceso en función de datos observables.

Aproximaciones metodológicas

  • Usar métodos que no requieren el uso de funciones matemáticas porque se basan en el conocimiento de expertos (técnicas cualitativas).
  • Utilización de procedimientos estadísticos para predecir valores futuros de una serie temporal en base a tendencias históricas (técnicas cuantitativas).
    • Uso de funciones modelo: Ajustar una función a los datos y usar esta función para hacer la previsión.
    • Técnicas que no usan funciones modelo.
  • Combinar técnicas cualitativas y cuantitativas

Recomendaciones generales

  • Comprobar que la ventana de modelización no incluye cambios estructurales.
  • Preferir modelos simples.
  • No inferir relaciones causa-efecto de los modelos.
  • El pronóstico siempre lleva asociado un cierto error, esta incertidumbre es mayor cuanto más alejado del momento actual se haga el pronóstico.

  • Para reducir la incertidumbre del pronóstico se recomienda combinar varios métodos de pronóstico.
  • Utilizar el conocimiento experto aunque no sea cuantitativo.
  • Recalcular el modelo cada vez que se disponga de nuevos datos.

Ajuste de modelos

Proceso de ajuste de un modelo

  • Determinar el modelo adecuado
    • A partir del conocimiento de las variables de los datos
    • A partir de la forma de los datos
  • Ajustar el modelo ≡ Determinar los coeficientes del modelo
    • Habitualmente usando la técnico de mínimos cuadrados ordinarios (OLS)
  • Comprobación del modelo y del ajuste
    • ¿Representa adecuadamente los datos?
    • ¿Los residuales son independientes y siguen una distribución normal?
    • ¿Representa adecuadamente datos no usados para construir el modelo?

Ajustes de modelos en R por OLS

Una vez decidido un modelo, el ajuste por el método de mínimos cuadrados ordinarios se realiza usando la notación de fórmulas y la función lm.

Notación de fórmulas

Una fórmula en R indica una relación lineal entre variables o combinaciones de variables.

  • Recta: y ~ x
  • Recta (sin ordenada en origen): y ~ x - 1 o y ~ 0 + x
  • Expresiones lineales multivariantes: y ~ x1 + x2
  • Expresiones multivariantes con interacciones: y ~ (X1 + x2) ^ 2
  • Polinomio: y ~ x + I(x ^ 2) o y ~ poly(x, 3)
  • Expresiones no lineales: y ~ exp(x) o log(y) ~ log(x)

Uso de la función lm

x <- seq(1, 2, length.out = 11)
y <- round(3 * x + 5 + rnorm(11) / 10, 3)
datos <- data.frame(x = x, y = y)
x y
1.0 7.944
1.1 8.277
1.2 8.756
1.3 8.907
1.4 9.213
1.5 9.672
1.6 9.846
1.7 9.973
1.8 10.331
1.9 10.655
2.0 11.122

ajuste <- lm(y ~ x, datos)
summary(ajuste)
## 
## Call:
## lm(formula = y ~ x, data = datos)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.14286 -0.06881 -0.01278  0.06913  0.15418 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   5.0325     0.1495   33.67 8.87e-11 ***
## x             2.9902     0.0975   30.67 2.04e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1023 on 9 degrees of freedom
## Multiple R-squared:  0.9905, Adjusted R-squared:  0.9895 
## F-statistic: 940.5 on 1 and 9 DF,  p-value: 2.041e-10
datos$yc <- ajuste$fitted
x y yc
1.0 7.944 8.022727
1.1 8.277 8.321746
1.2 8.756 8.620764
1.3 8.907 8.919782
1.4 9.213 9.218800
1.5 9.672 9.517818
1.6 9.846 9.816836
1.7 9.973 10.115854
1.8 10.331 10.414873
1.9 10.655 10.713891
2.0 11.122 11.012909
ggplot(datos, aes(x = x, y = y)) +
  geom_point() +
  geom_line(aes(y = yc)) + 
  theme_classic()

Análisis estadístico de los modelos de ajuste

Un ajuste es adecuado si

  • Se ajusta a los valores experimentales y su funcionalidad
    • Gráfico con el modelo ajustado
    • Análisi de los residuales
      • Residuales sin funcionalidad y centrado en cero
    • Prueba F
  • No hay puntos excesivamente influyentes
    • Vigilar con los puntos anómalos y los puntos extremos
  • Los errores están normalmente distribuidos
    • Prueba de normalidad de los residuales
    • Gráfico de normalidad de los residuales
  • Sus coeficientes son significativos
    • Estadísticamente
    • Teóricamente / científicamente

summary(ajuste)
## 
## Call:
## lm(formula = y ~ x, data = datos)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.14286 -0.06881 -0.01278  0.06913  0.15418 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   5.0325     0.1495   33.67 8.87e-11 ***
## x             2.9902     0.0975   30.67 2.04e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1023 on 9 degrees of freedom
## Multiple R-squared:  0.9905, Adjusted R-squared:  0.9895 
## F-statistic: 940.5 on 1 and 9 DF,  p-value: 2.041e-10
plot(ajuste, which=1:2)

shapiro.test(ajuste$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  ajuste$residuals
## W = 0.93572, p-value = 0.4716

Predicción

Si el modelo es correcto (y solo en este caso), este puede usarse para estimar valores de la variable dependiente.

predict(ajuste,data.frame(x = c(1.12,1.23)))
##        1        2 
## 8.381549 8.710469

A esta predicción puede añadirse intervalos de confianza (por defecto, a un nivel de confianza del 95 %). Estos se calculan de forma distinta cuando se quiere predecir el valor medio (mejor estimación) o un nuevo dato.

Si x es 1,12, la estimación del promedio es

predict(ajuste,data.frame(x = c(1.12)), interval="confidence")
##        fit      lwr      upr
## 1 8.381549 8.272507 8.490592

y la estimación para un nuevo dato es

predict(ajuste,data.frame(x = c(1.12)), interval="prediction")
##        fit      lwr      upr
## 1 8.381549 8.125803 8.637295

Ejemplo

A partir de los datos que figuran en el archivo “cu.txt”, determina el coeficiente de expansión térmica del cobre (k) a las temperaturas (T) de 600 K y de 200 K.

 

http://www.itl.nist.gov/div898/handbook/pmd/section6/pmd641.htm

Ajuste de series temporales

Serie temporal

Se denomina serie temporal a una secuencia de puntos indexados temporalmente. Es decir, obtenidos a intervalos de tiempo constantes.

En estos casos, el índice temporal (o posición) se incorpora como variable de modelización. El tiempo no se usa como variable independiente.

Solo se presentará el tratamiento del ajuste de serie temporales completas.

Fuentes de variación en una serie temporal

  • Tendencia. Evolución de los datos a largo plazo
  • Estacionalidad. Variaciones estacionales a corto plazo (periodo de menos de un año)
    • Día de la semana
    • Día del mes
    • Día del año
    • Semana del año
    • Semana del mes
    • Mes del año
  • Componentes aleatorios no previsibles
  • Autocorrelación. Dependencia de los periodos inmediatamente anteriores.

Previsión en series temporales (pronóstico, forecasting)

  • Ajuste de modelos (con variables fictícias)
  • Métodos de previsión sin funciones de ajuste (para datos desestacionalizados)

Pronóstico mediante ajuste de modelos

Ejemplo

Sean las ventas de un determinado producto entre enero de 2013 y mediados de julio de 2016 las que figuran en el archivo “ventas.csv”, ¿cuál seria una previsión razonable de las ventas de julio y agosto de 2016?

ano mes ventas
2013 1 3038
2013 2 3420
2013 3 3138
2013 4 4137
2013 5 4302
2013 6 3374

Visualización de los datos

datos$fecha <- paste(datos$ano,
                     formatC(datos$mes, width=2,
                             format="d", flag="0"),
                     sep="-")
ggplot(datos,aes(x=fecha,y=ventas,group=1)) + geom_line() +
  theme_classic()+
  theme(axis.text.x = element_text(angle=90,vjust=0.5))

Modelos habituales (sin estacionalidad)

  • Constante ( y ~ 1 )
  • Recta ( y ~ x )
  • Recta sin ordenada ( y ~ 0 + x )

También pueden usarse funciones no lineales (polinomios, logaritmos, exponenciales…).

Ajuste de una constante

datosAj <- datos[1:42,]
datosAj$i <- 1:42
datosAj <- rbind(datosAj,
                 list(ano=c(2016,2016),
                      mes=7:8,
                      ventas=c(NA,NA),
                      fecha=c("2016-07","2016-08"),
                      i=43:44))

ajuste <- lm(ventas~1, datosAj[1:42,])

datosAj <- cbind(datosAj,
                 predict(ajuste, datosAj, interval="prediction"))
ggplot(datosAj,aes(x=fecha,y=ventas,group=1)) + geom_line() +
  geom_line(aes(y=fit), col="green")+
  geom_ribbon(aes(ymin=lwr, ymax=upr), fill="green", alpha=0.05)+
  theme_classic()+
  theme(axis.text.x = element_text(angle=90,vjust=0.5))

print(datosAj[43:44,c(4,6:8)])
##      fecha      fit      lwr      upr
## 43 2016-07 3331.548 1740.936 4922.159
## 44 2016-08 3331.548 1740.936 4922.159

Ajuste de una recta

datosAj <- datos[1:42,]
datosAj$i <- 1:42
datosAj <- rbind(datosAj,
                 list(ano=c(2016,2016),
                      mes=7:8,
                      ventas=c(NA,NA),
                      fecha=c("2016-07","2016-08"),
                      i=43:44))

ajuste <- lm(ventas~i, datosAj[1:42,])

datosAj <- cbind(datosAj,
                 predict(ajuste, datosAj, interval="prediction"))
ggplot(datosAj,aes(x=fecha,y=ventas,group=1)) + geom_line() +
  geom_line(aes(y=fit), col="green")+
  geom_ribbon(aes(ymin=lwr, ymax=upr), fill="green", alpha=0.05)+
  theme_classic()+
  theme(axis.text.x = element_text(angle=90,vjust=0.5))

print(datosAj[43:44,c(4,6:8)])
##      fecha      fit      lwr      upr
## 43 2016-07 2976.685 1364.645 4588.726
## 44 2016-08 2960.180 1342.917 4577.443

Incorporar la estacionalidad

La estacionalidad se incorpora a los modelos a través de variables ficticias

  • Se decide el periodo y el paso de la periodicidad
    • Día en mes
    • Día en semana
    • Mes en año
    • Trimestre en año
  • Se usa una variable ficticia (lógica, 0 o 1) para cada elemento del periodo
  • Se introducen al modelo todas las variables ficticias excepto una (que corresponderá a la referencia del modelo).

Modelizar la estacionalidad en R

  • Definir una variable para la estacionalidad (valores cíclicos).
  • Se establece esta variable como factor y se introduce como variable del modelo.

Ajuste de una recta más una periodicidad anual (por meses)

datosAj <- datos[1:42,]
datosAj$i <- 1:42
datosAj <- rbind(datosAj,
                 list(ano=c(2016,2016),
                      mes=7:8,
                      ventas=c(NA,NA),
                      fecha=c("2016-07","2016-08"),
                      i=43:44))
datosAj$mes <- factor(datosAj$mes)

ajuste <- lm(ventas~i+mes, datosAj[1:42,])

datosAj <- cbind(datosAj,
                 predict(ajuste, datosAj, interval="prediction"))
ggplot(datosAj,aes(x=fecha,y=ventas,group=1)) + geom_line() +
  geom_line(aes(y=fit), col="green")+
  geom_ribbon(aes(ymin=lwr, ymax=upr), fill="green", alpha=0.05)+
  theme_classic()+
  theme(axis.text.x = element_text(angle=90,vjust=0.5))

print(datosAj[43:44,c(4,6:8)])
##      fecha       fit       lwr      upr
## 43 2016-07 3996.2143 3373.7414 4618.687
## 44 2016-08  932.2143  309.7414 1554.687

Comprobación del ajuste

summary(ajuste)
## 
## Call:
## lm(formula = ventas ~ i + mes, data = datosAj[1:42, ])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -463.64 -101.52  -22.13  128.66  538.91 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3071.358    141.707  21.674  < 2e-16 ***
## i             -18.769      3.274  -5.732 3.33e-06 ***
## mes2          404.769    180.088   2.248 0.032375 *  
## mes3          530.038    180.177   2.942 0.006356 ** 
## mes4          737.557    180.326   4.090 0.000313 ***
## mes5          785.575    180.534   4.351 0.000153 ***
## mes6          736.344    180.801   4.073 0.000328 ***
## mes7         1731.917    194.485   8.905 8.56e-10 ***
## mes8        -1313.314    194.512  -6.752 2.07e-07 ***
## mes9         1593.788    194.595   8.190 4.96e-09 ***
## mes10        1171.557    194.733   6.016 1.52e-06 ***
## mes11         750.325    194.925   3.849 0.000601 ***
## mes12        1098.761    195.173   5.630 4.42e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 254.6 on 29 degrees of freedom
## Multiple R-squared:  0.9243, Adjusted R-squared:  0.893 
## F-statistic: 29.51 on 12 and 29 DF,  p-value: 5.199e-13
plot(ajuste, which = 1:2)

shapiro.test(ajuste$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  ajuste$residuals
## W = 0.98166, p-value = 0.7257

Pronósticos a corto plazo (para datos desestacionalizados)

Datos desestacionalizados

  • No tienen tendencia ni funcionalidad.
  • No presentan variabilidad estacional.

Presentan autocorrelación, es decir dependencia de los valores anteriores.

Métodos simples para predecir a corto plazo

  • Método Naïve (o Naïve I)
  • Método Naïve II
  • Método de la media móvil simple

Método Naïve

Consiste en pronosticar el siguiente valor de la serie temporal como igual al último valor disponible en la misma

\[ \\ \hat y_{n+1} = y_n \]

Método Naïve II

Consiste en pronosticar el siguiente valor de la serie temporal como igual al último valor disponible en la misma más un término dependiente de la variación entre los dos valores anteriores.

\[ \\ \hat y_{n+1} = y_n + k (y_n - y_{n-1}) \]

Si k vale 1, consiste en el ajuste de una recta a los dos datos anteriores.

Media móvil simple

Consiste en pronosticar el siguiente valor de la serie temporal como media de los m valores anteriores.

\[ \\ \hat{y}_{n+1} = \sum\limits_{i=n-m+1}^{n} \frac{y_i}{m} \]

Comprobación del ajuste

El ajuste se puede comprobar aplicando el método a los valores anteriores de la serie y calculando los residuales correspondientes.

El análisis de residuales sigue siendo válido en estos casos.

Estimación del intervalo de predicción

Si los residuales estan normalmente distribuidos, el intervalo de predicción puede calcularse a partir de la desviación estándar de los residuales.

\[ \\ y = \hat y \pm z_{1-\alpha/2} \ s_{res} \]

Ejemplo

Disponemos de los precios históricos (promedios mensuales) de uno de los productos de partida de nuestro proceso de fabricación, ¿podemos estimar el precio que nos costará este recurso en el siguiente mes?

Los datos figuran en el archivo “coconutOil.txt”.

Criterios de error para la comparación de modelos

Para establecer la calidad de un modelo y tener criterios que permitan comparar modelos, se usan distintos indicadores

  • Coeficiente de determinación ajustado, R2
  • Distintas medidas del error, calculadas a partir de los residuales (MAPE/MAPD , MSE/MSD…)

 

https://support.minitab.com/en-us/minitab/17/topic-library/modeling-statistics/time-series/time-series-models/what-are-mape-mad-and-msd/

Referencias adicionales