Cómo funciona la predicción de suavizado exponencial

La herramienta Predicción de suavizado exponencial utiliza el método de suavizado exponencial de Holt-Winters para descomponer las series temporales de cada ubicación de un cubo de espacio-tiempo en componentes estacionales y de tendencia para pronosticar eficazmente los periodos de tiempo futuros de cada ubicación. La salida principal es un mapa del periodo de tiempo predicho final, además de mensajes informativos y gráficos emergentes. También puede crear un nuevo cubo de espacio-tiempo que contenga los datos del cubo original junto con valores predichos anexados. También tiene la posibilidad de detectar valores atípicos en cada serie temporal para identificar posiciones y tiempos que se desvíen significativamente de los patrones y tendencias del resto de las series temporales.

El modelo de suavizado exponencial se utiliza para predecir periodos de tiempo futuros.

Aplicaciones potenciales

El suavizado exponencial es uno de los métodos de predicción de series temporales más antiguos y estudiados. Su máxima eficacia se obtiene cuando los valores de las series temporales siguen una tendencia gradual y presentan un comportamiento estacional en el cual los valores siguen un patrón cíclico repetido a lo largo de un número determinado de periodos de tiempo.

Por ejemplo, puede utilizar esta herramienta en las siguientes aplicaciones:

• Un responsable de salud municipal puede usar esta herramienta para predecir la temperatura por hora del centro urbano durante una ola de calor para prepararse para las enfermedades relacionadas con el calor.
• Una cadena de tiendas puede utilizar esta herramienta para predecir la demanda de elementos individuales para cada día de la semana siguiente.

Previsión y validación

La herramienta crea dos modelos al realizar la predicción de cada serie temporal. El primero es el modelo de previsión, que se utiliza para realizar la previsión de los valores de futuros periodos de tiempo. El segundo es el modelo de validación, que se utiliza para validar los valores previstos.

Modelo de previsión

El modelo de predicción se construye realizando un suavizado exponencial con los valores de la serie temporal en cada ubicación del cubo de espacio-tiempo. A continuación, este modelo se utiliza para predecir periodos de tiempo futuros. El ajuste del modelo de suavizado exponencial a cada serie temporal se mide mediante el error cuadrático medio (RMSE) de predicción, que es igual a la raíz cuadrada de la diferencia cuadrada promedio entre el modelo de suavizado exponencial y los valores de la serie temporal.

, donde T es el número de periodos de tiempo, c_t es el valor ajustado del suavizado exponencial y r_t es el valor sin procesar de la serie temporal en el tiempo t.

La siguiente imagen muestra los valores brutos de una serie temporal y un modelo de suavizado exponencial ajustado a la serie temporal junto con predicciones para dos futuros periodos de tiempo. El RMSE de predicción mide cuánto varían los valores ajustados del modelo con respecto a los valores brutos de la serie temporal.

El RMSE de predicción solo mide el grado con el que el modelo de suavizado exponencial se ajusta a los valores brutos de la serie temporal. No mide lo bien que el modelo de predicción realiza realmente la predicción de valores futuros. Es habitual que los modelos se ajusten fielmente a una serie temporal, aunque no proporcionen predicciones precisas cuando se extrapolan. Este problema se aborda con el modelo de validación.

Modelo de validación

El modelo de validación se utiliza para determinar lo bien que el modelo de predicción puede realizar la predicción de valores futuros de cada serie temporal. Se construye excluyendo algunos de los últimos periodos de tiempo de cada serie temporal y ajustando el modelo de suavizado exponencial a los datos no excluidos. A continuación, el modelo se utiliza para realizar la predicción de los valores de los datos que se retuvieron, y los valores predichos se comparan con los valores brutos que se ocultaron. De forma predeterminada, se retiene el 10 por ciento de los periodos de tiempo para la validación, pero este número se puede modificar con el parámetro Número de periodos de tiempo que se excluirán para la validación. El número de periodos de tiempo excluidos no puede superar el 25 por ciento del número de periodos de tiempo, y no se realiza ninguna validación si se especifica 0. La precisión de las predicciones se mide calculando una estadística de RMSE de validación, que es igual a la raíz cuadrada de la diferencia cuadrada promedio entre los valores predichos y brutos de los periodos de tiempo excluidos.

, donde T es el número de periodos de tiempo, m es el número de periodos de tiempo retenidos para validación, c_t es el valor predicho de los primeros T-m periodos de tiempo y r_t es el valor bruto de la serie temporal retenida para validación en el tiempo t.

La siguiente imagen muestra un modelo de suavizado exponencial ajustado a la primera mitad de una serie temporal y usado para predecir la segunda mitad de la serie temporal. El RMSE de validación mide cuánto varían los valores predichos de los valores brutos en los periodos de tiempo retenidos.

El modelo de validación es importante porque puede comparar directamente valores predichos con valores brutos para medir lo bien que puede predecir el modelo de suavizado exponencial. Aunque no se utiliza realmente para realizar predicción, se utiliza para justificar el modelo de predicción.

Nota:

La validación en la predicción de serie temporal es similar aunque no idéntica a una técnica habitual denominada validación cruzada. La diferencia es que la validación de predicción siempre excluye los periodos de tiempo finales para validación, y la validación cruzada excluye un subconjunto aleatorio de los datos o excluye cada valor secuencialmente.

Interpretación

Existen varias consideraciones al interpretar los valores de RMSE de predicción y RMSE de validación.

• Los valores de RMSE no son directamente comparables entre sí porque miden cosas diferentes. El RMSE de predicción mide el ajuste del modelo a los valores brutos de la serie temporal, y el RMSE de validación mide lo bien que el modelo puede realizar la predicción de valores futuros. Dado que el RMSE de predicción utiliza más datos y no extrapola, generalmente es menor que el RMSE de Validación.
• Ambos valores de RMSE están en las unidades de los datos. Por ejemplo, si sus datos son mediciones de temperatura en grados Celsius, un RMSE de validación de 50 es muy alto, porque significa que los valores predichos difirieron de los valores reales en una media de aproximadamente 50 grados. No obstante, si sus datos son ingresos diarios en dólares estadounidenses de una tienda minorista grande, el mismo RMSE de validación de 50 es muy pequeño porque significa que los ingresos diarios predichos solo varían con respecto a los valores reales en 50 $ al día en promedio.

Construcción del modelo de suavizado exponencial

Existen distintos tipos de suavizado exponencial, pero todos se basan en separar las series temporales en varios componentes. Los valores de cada componente se estiman ponderando exponencialmente los componentes de los periodos de tiempo anteriores, de manera que la influencia de cada periodo de tiempo disminuya exponencialmente a medida que se avanza en el tiempo. Cada componente se define recursivamente a través de un planteamiento de modelo de estado-espacio, y cada componente depende de todos los demás componentes. La estimación de todos los parámetros se realiza mediante una estimación de máxima probabilidad.

En esta herramienta, todos los componentes son aditivos, de modo que el modelo de predicción es la suma de los componentes individuales. Si se utiliza un componente estacional, la herramienta utiliza el método estacional amortiguado de Holt-Winters. Si no se utiliza ningún componente estacional, la herramienta utiliza el método de tendencia amortiguada. Encontrará los detalles de estos componentes y las ecuaciones que definen los modelos de estado-espacio en el libro de texto que se indica en la sección Referencias adicionales.

Componente de tendencia

El primer componente del modelo de suavizado exponencial es el componente de tendencia. Este componente se utiliza para modelar cambios graduales y sistemáticos en los valores de las series temporales. La estimación se realiza ponderando exponencialmente los valores de cada periodo de tiempo por la diferencia entre su valor y el valor del periodo de tiempo anterior. El componente de tendencia se utiliza directamente al realizar predicciones a partir de la última tendencia detectada por el modelo. Sin embargo, para evitar que las predicciones sigan ilimitadamente la tendencia final, la tendencia se amortigua de manera que se aplana gradualmente a medida que se avanza en el tiempo. La amortiguación aplana la tendencia multiplicando la pendiente del valor de la tendencia de cada periodo de tiempo por un valor exponencialmente decreciente. El modelo realiza una estimación del nivel de amortiguación, por lo que la tendencia puede aplanarse más rápida o lentamente (o no aplanarse en absoluto en el caso más extremo) en unos modelos que en otros al realizar predicciones para momentos futuros más lejanos.

Componente estacional

El segundo componente del suavizado exponencial es el componente estacional, que se utiliza para modelar patrones en los datos que se repiten a lo largo de un número determinado de periodos de tiempo. La forma y la magnitud del patrón dentro de cada estación puede cambiar con el tiempo, pero la duración de cada estación debe ser la misma para toda la serie temporal. Por ejemplo, la temperatura muestra un comportamiento estacional que corresponde a días y noches con temperaturas más bajas durante la noche y más altas durante el día. Aunque el amanecer se produce a horas distintas del día a lo largo del año (y, por tanto, la forma y la magnitud del patrón de temperatura cambian dentro de un mismo día), la duración de una estación es siempre un único día.

Al igual que con el componente de tendencia, el componente estacional de un determinado periodo de tiempo se determina ponderando exponencialmente los valores estacionales de los periodos de tiempo anteriores. Sin embargo, en lugar de utilizar los periodos de tiempo inmediatamente anteriores, solo pondera los periodos de tiempo anteriores correspondientes al mismo punto de un ciclo estacional. Por ejemplo, si la duración de una estación es de cuatro periodos de tiempo, el componente estacional pondera exponencialmente los valores correspondientes a 4 periodos de tiempo anteriores, 8 periodos de tiempo anteriores, 12 periodos de tiempo anteriores, y así sucesivamente.

Si conoce el número de periodos de tiempo que corresponden a una estación en sus datos, puede especificarlos en el parámetro Duración de estación, y este valor se utilizará con todas las ubicaciones del cubo de espacio-tiempo. Si no conoce la duración de una estación o si la duración estacional es diferente para diferentes ubicaciones, el valor del parámetro puede dejarse vacío, y se estimará una duración estacional óptima para cada ubicación utilizando una función de densidad espectral. Para disponer de detalles sobre esta función, consulte la sección Recursos adicionales.

En el caso de una ubicación individual, si la duración óptima de la estación determinada por el análisis espectral es mayor que uno y menor que un tercio del número de periodos de tiempo de la ubicación, la duración de la estación se establece en este valor óptimo. De lo contrario, la ubicación no utiliza ningún componente estacional. La duración de estación utilizada en la ubicación se guarda en el campo Duración de estación de las entidades de salida. Si no se utiliza ningún componente estacional, el valor de este campo es 1. Este flujo de trabajo se resume en la siguiente imagen:

La duración de la estación se puede especificar o estimar para cada ubicación.

Componente de nivel

El componente de nivel del suavizado exponencial representa el valor de línea base de la serie temporal teniendo en cuenta la estacionalidad y la tendencia. Al ajustar el modelo de predicción al cubo de espacio-tiempo de entrada, el nivel de un periodo de tiempo se calcula ponderando exponencialmente los niveles anteriores y teniendo en cuenta la estacionalidad y la tendencia. Al hacer predicciones a futuro, el componente de nivel se mantiene igual al componente de nivel del último periodo de tiempo medido y, en este caso, las predicciones reales dependen de los componentes de tendencia y estacional.

Componente residual e intervalos de confianza

El último componente es el componente residual (o de error). Este componente es la diferencia entre el valor real y el valor estimado por todos los demás componentes. Representa la incertidumbre y el error restantes en los datos una vez modelados los componentes de tendencia, estación y nivel. Este componente es importante porque constituye la base de los intervalos de confianza.

Para cada periodo de tiempo predicho, la herramienta calcula los límites superiores e inferiores de un intervalo de confianza del 90 por ciento para el valor predicho. El valor predicho en cada periodo de tiempo representa la mejor estimación individual del valor futuro, pero se puede utilizar el intervalo de confianza para visualizar la incertidumbre y el rango probable del valor futuro verdadero. Los límites superior e inferior se guardan como campos y se muestran en gráficos emergentes de las Entidades de salida.

En la estimación de los intervalos de confianza, se presupone que los residuales del modelo están distribuidos normalmente, de forma independiente e idéntica. Partiendo de esta suposición, se pueden derivar fórmulas para los intervalos de confianza. Encontrará las fórmulas y sus derivaciones en el libro de texto indicado en la sección Referencias adicionales.

Visualizar los componentes

Es posible visualizar los componentes de su modelo de suavizado exponencial creando un cubo de espacio-tiempo de salida. Use este cubo en la herramienta Visualizar cubo de espacio-tiempo en 3D con la opción Resultados de la predicción en el parámetro Tema de visualización. Se crea un gráfico para las entidades de salida, y es posible activar y desactivar los distintos componentes del modelo de suavizado exponencial en el panel Propiedades de gráfico. Al sumar estos componentes, construyen el modelo de predicción y las predicciones para los periodos de tiempo futuros. La imagen que aparece a continuación muestra los componentes individuales del modelo de suavizado exponencial que se muestra en la primera imagen de este tema:

El modelo de suavizado exponencial está separado en varios componentes.

Identificación de valores atípicos de series temporales

Los valores atípicos de los datos de series temporales son valores que varían significativamente de los patrones y tendencias de los demás valores en series temporales. Por ejemplo, las grandes cifras de compras en línea por las vacaciones o las altas cifras de accidentes de tráfico durante fuertes tormentas probablemente se detecten como valores atípicos en sus series temporales. Errores simples de introducción de datos, como la omisión de la parte decimal de un número, representan otra fuente común de valores atípicos. La identificación de valores atípicos en la previsión de series temporales es importante porque los valores atípicos influyen en el modelo de previsión que se utiliza para prever futuros valores, e incluso un pequeño número de valores atípicos en las series temporales de una ubicación pueden reducir significativamente la precisión y fiabilidad de las previsiones. Las ubicaciones con valores atípicos, particularmente valores atípicos en el comienzo y final de las series temporales, pueden generar previsiones erróneas, y la identificación de estas ubicaciones ayuda a determinar el grado de confianza que debe tener en los valores previstos en cada ubicación.

Los valores atípicos no se determinan simplemente mediante sus valores sin procesar, sino por cuánto varían sus valores con respecto a los valores ajustados del modelo de previsión. Eso significa que la posibilidad de que se determine que un valor es atípico será algo contextual y dependerá de su posición y tiempo. El modelo de previsión define el valor que se espera en función de toda la serie temporal, y los valores atípicos son los valores que se desvían significativamente de esta referencia. Por ejemplo, considere una serie temporal de temperatura media anual. Ya que las temperaturas promedio han aumentado durante las últimas décadas, el modelo de previsión ajustado de temperatura también aumentará con el tiempo para reflejar este aumento. Eso significa que un valor de temperatura que se consideraría típico y no un valor atípico en 1950 es probable que se considere un valor atípico si se produjera la misma temperatura en 2020. En otras palabras, una temperatura típica de 1950 se consideraría muy baja para los estándares de 2020.

Puede elegir detectar valores atípicos de series temporales en cada ubicación usando el parámetro Identificar valores atípicos. Si se especifica, la prueba Extreme Studentized Deviate (ESD) generalizada se realiza para cada ubicación con el fin de probar valores atípicos de series temporales. El nivel de confianza de la prueba puede especificarse con el parámetro Nivel de confianza, y se utiliza el 90 por ciento de confianza como el valor predeterminado. La prueba ESD generalizada prueba iterativamente un único valor atípico, dos valores atípicos, tres valores atípicos, y así sucesivamente, en cada ubicación hasta el valor del parámetro Número máximo de valores atípicos (de forma predeterminada, 5 por ciento del número de periodos de tiempo, redondeado a la baja), y se devuelve el número de valores atípicos estadísticamente significativo más grande. El número de valores atípicos en cada ubicación puede verse en una tabla de atributos de las entidades de salida, y los valores atípicos individuales pueden verse en los gráficos emergentes de la serie temporal que se analizan en la siguiente sección.

Conozca más información sobre valores atípicos en el análisis de series temporales, la prueba ESD generalizada y cómo interpretar los resultados

Salidas de la herramienta

La salida principal de esta herramienta es una clase de entidad 2D que muestra cada ubicación en el Cubo de espacio-tiempo de entrada simbolizado por el periodo de tiempo predicho final con las predicciones para todos los demás periodos de tiempo almacenados como campos. Aunque se realiza la predicción de cada ubicación de forma independiente y las relaciones espaciales no se tienen en cuenta, el mapa puede mostrar patrones espaciales para áreas con series temporales similares.

Gráficos emergentes

Al hacer clic en cualquier entidad del mapa usando la herramienta de navegación Explorar, se muestra un gráfico en el panel Emergente que muestra los valores del cubo de espacio-tiempo junto con el modelo de suavizado exponencial ajustado y los valores predichos junto con intervalos de confianza del 90 por ciento para cada predicción. Los valores del cubo de espacio-tiempo se muestran en azul y se conectan con una línea azul. Los valores ajustados se muestran en naranja y se conectan mediante una línea naranja discontinua. Los valores predichos se muestran en naranja y se conectan mediante una línea naranja continua que representa la predicción del modelo. Los límites de confianza rojo claro se dibujan alrededor de cada valor predicho. Puede situar el cursor en cualquier punto del gráfico para ver la fecha y valor del punto. Adicionalmente, si elige detectar valores atípicos en series temporales, los valores atípicos se muestran como puntos morados grandes.

El gráfico emergente muestra los valores originales, valores ajustados, valores predichos e intervalos de confianza.

Nota:

Los gráficos emergentes no se crean cuando las entidades de salida se guardan como un shapefile (.shp).

Mensajes de geoprocesamiento

La herramienta proporciona varios mensajes con información sobre la ejecución de la herramienta. Los mensajes tienen tres secciones principales.

La sección Detalles de cubo de espacio-tiempo de entrada muestra propiedades del cubo de espacio-tiempo de entrada junto con información sobre el número de periodos de tiempo, número de ubicaciones y número de bins espacio-tiempo. Las propiedades mostradas en esta primera sección dependen de cómo se creó originalmente el cubo, por lo que la información varía en función del cubo.

La sección Detalles del análisis muestra propiedades de los resultados de la predicción, incluidos el número de periodos de tiempo predichos, el número de periodos de tiempo excluidos para la validación, el porcentaje de ubicaciones que presentan estacionalidad e información sobre los periodos de tiempo predichos. Si no se proporciona ningún valor para el parámetro Duración de estación, se muestran estadísticas resumidas de la duración estimada de la estación, incluidas la mínima, máxima, media, mediana y desviación estándar.

La sección Resumen de precisión a través de ubicaciones muestra estadísticas resumen para el RMSE de previsión y RMSE de validación entre todas las ubicaciones. Para cada valor se muestra el mínimo, máximo, valor medio, mediana y desviación estándar.

La sección Resumen de los valores atípicos de series temporales aparece si elige detectar valores atípicos de series temporales usando el parámetro Opción de valores atípicos. Esta sección muestra información que incluye el número y porcentaje de ubicaciones que contienen valores atípicos, el periodo de tiempo que contiene la mayoría de valores atípicos y estadísticas resumen del número de valores atípicos por ubicación y por periodo de tiempo.

Nota:

Los mensajes de geoprocesamiento aparecen en la parte inferior del panel Geoprocesamiento durante la ejecución de la herramienta. Puede acceder a los mensajes desplazándose sobre la barra de progreso, haciendo clic en el botón emergente o expandiendo la sección de mensajes en el panel Geoprocesamiento. También puede acceder a los mensajes de una herramienta ejecutada anteriormente con el historial de geoprocesamiento.

Campos de las entidades de salida

Además de Id. de objeto, campos de geometría y el campo que contiene los gráficos emergentes, las Entidades de salida tendrán los siguientes campos:

• Id. de ubicación (LOCATION): el Id. de ubicación de la ubicación correspondiente del cubo de espacio-tiempo.
• Predicción para (variables de análisis) en (periodo de tiempo) (FCAST_1, FCAST_2 y así sucesivamente): el valor predicho de cada periodo de tiempo futuro. El alias de campo muestra el nombre de la Variable de análisis y la fecha de la predicción. Se crea un campo de este tipo para cada periodo de tiempo previsto.
• Intervalo alto para (variable de análisis) en (periodo de tiempo) (HIGH_1, HIGH_2 y así sucesivamente): el límite superior de un intervalo de confianza del 90 por ciento para el valor predicho de cada periodo de tiempo futuro. El alias de campo muestra el nombre de la Variable de análisis y la fecha de la predicción. Se crea un campo de este tipo para cada periodo de tiempo previsto.
• Intervalo bajo para (variable de análisis) en (periodo de tiempo) (LOW_1, LOW_2 y así sucesivamente): el límite inferior de un intervalo de confianza del 90 por ciento para el valor predicho de cada periodo de tiempo futuro. El alias de campo muestra el nombre de la Variable de análisis y la fecha de la predicción. Se crea un campo de este tipo para cada periodo de tiempo previsto.
• Error cuadrático medio de predicción (F_RMSE): el RMSE de predicción.
• Error cuadrático medio de validación (V_RMSE): el RMSE de validación. Si no se excluyeron periodos de tiempo para validación, este campo no se crea.
• Duración de estación (SEASON): el número de periodos de tiempo correspondientes a una estación en la ubicación. Un valor de 1 en este campo significa que no hay estacionalidad.
• Método de predicción (METHOD): un campo de texto que muestra el modelo empleado en la ubicación. Para esta herramienta, el valor es siempre suavizado exponencial. Este campo le permite identificar los modelos que se utilizan en la herramienta Evaluar predicciones por ubicación.
• Número de valores atípicos de ajuste de modelo (N_OUTLIERS): el número de valores atípicos detectados en la serie temporal de la ubicación. Este campo solo se crea si elige detectar valores atípicos con el parámetro Opción de valores atípicos.

Cubo de espacio-tiempo de salida

Si se especifica un Cubo de espacio-tiempo de salida, el cubo de salida contiene todos los valores originales del cubo de espacio-tiempo de entrada con los valores predichos incorporados. Este nuevo cubo de espacio-tiempo se puede visualizar con las herramientas Visualizar cubo de espacio-tiempo en 2D o Visualizar cubo de espacio-tiempo en 3D, y se puede utilizar como entrada para las herramientas de la caja de herramientas Minería de patrones en espacio-tiempo, tales como Análisis de puntos calientes emergentes y Clustering de serie temporal.

Con la herramienta Evaluar previsiones por ubicación se pueden comparar y fusionar varios cubos de espacio-tiempo previstos. De este modo, puede crear varios cubos de predicción con distintos parámetros y herramientas de predicción, y la herramienta identifica la mejor predicción para cada ubicación mediante el RMSE de Predicción o de Validación.

Prácticas recomendadas y limitaciones

Al decidir si esta herramienta resulta apropiada para sus datos y los parámetros que debe elegir, deben tenerse en cuenta varias cosas.

• En comparación con otras herramientas de predicción del conjunto de herramientas Predicción de serie temporal, esta herramienta se recomienda para datos que presenten tendencias moderadas y un comportamiento estacional marcado. En el modelo exponencial se presupone que es posible separar el comportamiento estacional de la tendencia, por lo que resulta más eficaz con datos cuya tendencia cambie gradualmente y siga patrones estacionales constantes a lo largo del tiempo. El componente estacional del modelo es opcional, por lo que esta herramienta se puede utilizar para datos que no presenten estacionalidad, pero funciona mejor en presencia de un comportamiento estacional marcado.
• La decisión sobre cuántos periodos de tiempo deben excluirse para la validación es importante. Cuantos más periodos de tiempo se excluyan, menores periodos de tiempo se mantienen para estimar el modelo de validación. No obstante, si se excluyen muy pocos periodos de tiempo, el RMSE de validación se estima usando una pequeña cantidad de datos y puede llevar a engaño. Se recomienda excluir el mayor número de periodos de tiempo que sea posible mientras aún se mantienen suficientes periodos de tiempo para estimar el modelo de validación. También se recomienda retener al menos el mismo número de periodos de tiempo para validación que el número de periodos de tiempo de los que pretende realizar la predicción, si su cubo de espacio-tiempo tiene suficientes periodos de tiempo para permitir esto.

Recursos adicionales

Para obtener más información sobre la predicción con suavizado exponencial con un planteamiento de estado-espacio, consulte el siguiente libro de texto:

• Hyndman R, Koehler A, Ord K, and Snyder R (2008). "Forecasting with Exponential Smoothing. The State Space Approach." https://doi.org/10.1007/978-3-540-71918-2

Para más información sobre la función de densidad espectral empleada para la estimación de la duración de una estación, consulte la función findfrequency en las siguientes referencias:

• Hyndman R, Athanasopoulos G, Bergmeir C, Caceres G, Chhay L, O'Hara-Wild M, Petropoulos F, Razbash S, Wang E and Yasmeen F (2019). "Forecasting functions for time series and linear models". Paquete R versión 8.7, https://pkg.robjhyndman.com/forecast.
• Hyndman RJ and Khandakar Y (2008). "Automatic time series forecasting: the forecast package for R". Journal of Statistical Software, 26(3), pp. 1–22. https://www.jstatsoft.org/article/view/v027i03.