Algoritmo simple para la detección de valores atípicos en línea de una serie temporal genérica

Estoy trabajando con una gran cantidad de series de tiempo. Estas series temporales son básicamente mediciones de red que se realizan cada 10 minutos, y algunas de ellas son periódicas (es decir, el ancho de banda), mientras que otras no (es decir, la cantidad de tráfico de enrutamiento).

Me gustaría un algoritmo simple para hacer una "detección de valores atípicos" en línea. Básicamente, quiero mantener en la memoria (o en el disco) todos los datos históricos para cada serie de tiempo, y quiero detectar cualquier valor atípico en un escenario en vivo (cada vez que se captura una nueva muestra). ¿Cuál es la mejor manera de lograr estos resultados?

Actualmente estoy usando un promedio móvil para eliminar algo de ruido, pero ¿qué sigue? Cosas simples como la desviación estándar, loco, ... contra todo el conjunto de datos no funciona bien (no puedo asumir que las series temporales son estacionarias), y me gustaría algo más "preciso", idealmente un cuadro negro como:

double outlier_detection (vector doble *, valor doble);

donde vector es la matriz de doble que contiene los datos históricos, y el valor de retorno es el puntaje de anomalía para el nuevo "valor" de la muestra.

Aquí hay una función R simple que encontrará valores atípicos de series temporales (y opcionalmente los mostrará en una gráfica). Manejará series temporales estacionales y no estacionales. La idea básica es encontrar estimaciones sólidas de la tendencia y los componentes estacionales y restarlos. Luego encuentre valores atípicos en los residuos. La prueba para valores atípicos residuales es la misma que para el diagrama de caja estándar: se supone que los puntos superiores a 1.5 IQR por encima o por debajo de los cuartiles superior e inferior son valores atípicos. El número de IQR por encima / debajo de estos umbrales se devuelve como un "puntaje" atípico. Por lo tanto, la puntuación puede ser cualquier número positivo y será cero para los no atípicos.

Me doy cuenta de que no está implementando esto en R, pero a menudo encuentro que una función R es un buen lugar para comenzar. Entonces la tarea es traducir esto a cualquier idioma que se requiera.

tsoutliers <- function(x,plot=FALSE)
{
    x <- as.ts(x)
    if(frequency(x)>1)
        resid <- stl(x,s.window="periodic",robust=TRUE)$time.series[,3]
    else
    {
        tt <- 1:length(x)
        resid <- residuals(loess(x ~ tt))
    }
    resid.q <- quantile(resid,prob=c(0.25,0.75))
    iqr <- diff(resid.q)
    limits <- resid.q + 1.5*iqr*c(-1,1)
    score <- abs(pmin((resid-limits[1])/iqr,0) + pmax((resid - limits[2])/iqr,0))
    if(plot)
    {
        plot(x)
        x2 <- ts(rep(NA,length(x)))
        x2[score>0] <- x[score>0]
        tsp(x2) <- tsp(x)
        points(x2,pch=19,col="red")
        return(invisible(score))
    }
    else
        return(score)
}

Una buena solución tendrá varios ingredientes, que incluyen:

• Use una ventana móvil resistente y suave para eliminar la no estacionariedad.

• Vuelva a expresar los datos originales para que los residuos con respecto a la suavidad se distribuyan aproximadamente simétricamente. Dada la naturaleza de sus datos, es probable que sus raíces cuadradas o logaritmos den residuales simétricos.

• Aplique métodos de gráfico de control, o al menos pensamiento de gráfico de control, a los residuos.

En lo que respecta al último, el pensamiento de la tabla de control muestra que los umbrales "convencionales" como 2 SD o 1.5 veces el IQR más allá de los cuartiles funcionan mal porque activan demasiadas señales falsas fuera de control. Las personas generalmente usan 3 SD en el trabajo de la tabla de control, por lo que 2.5 (o incluso 3) veces el IQR más allá de los cuartiles sería un buen punto de partida.

He esbozado más o menos la naturaleza de la solución de Rob Hyndman al mismo tiempo que le agrego dos puntos principales: la necesidad potencial de volver a expresar los datos y la sabiduría de ser más conservador al señalar un valor atípico. Sin embargo, no estoy seguro de que Loess sea bueno para un detector en línea, porque no funciona bien en los puntos finales. En su lugar, puede usar algo tan simple como un filtro de mediana móvil (como en el suavizado resistente de Tukey). Si los valores atípicos no vienen en ráfagas, puede usar una ventana estrecha (5 puntos de datos, tal vez, que se descompondrán solo con una explosión de 3 o más valores atípicos dentro de un grupo de 5).

Una vez que haya realizado el análisis para determinar una buena reexpresión de los datos, es poco probable que necesite cambiar la reexpresión. Por lo tanto, su detector en línea realmente solo necesita hacer referencia a los valores más recientes (la última ventana) porque no utilizará los datos anteriores. Si tiene series de tiempo realmente largas, podría ir más allá para analizar la autocorrelación y la estacionalidad (como las fluctuaciones diarias o semanales recurrentes) para mejorar el procedimiento.

(Esta respuesta respondió a una pregunta duplicada (ahora cerrada) en Detección de eventos pendientes , que presentó algunos datos en forma gráfica).

La detección de valores atípicos depende de la naturaleza de los datos y de lo que esté dispuesto a asumir sobre ellos. Los métodos de uso general se basan en estadísticas sólidas. El espíritu de este enfoque es caracterizar la mayor parte de los datos de una manera que no esté influenciada por valores atípicos y luego señalar cualquier valor individual que no se ajuste a esa caracterización.

Debido a que esta es una serie temporal, agrega la complicación de la necesidad de (re) detectar valores atípicos de manera continua. Si esto se va a hacer a medida que se desarrolla la serie, entonces solo se nos permite usar datos más antiguos para la detección, ¡no datos futuros! Además, como protección contra las muchas pruebas repetidas, nos gustaría utilizar un método que tenga una tasa muy baja de falsos positivos.

Estas consideraciones sugieren ejecutar una prueba atípica de ventana móvil simple y robusta sobre los datos . Hay muchas posibilidades, pero una simple, fácil de entender e implementar, se basa en un MAD en ejecución: desviación media absoluta de la mediana. Esta es una medida de variación muy sólida dentro de los datos, similar a una desviación estándar. Un pico periférico sería de varios MAD o más que la mediana.

R$x = (1,2,\ldots,n)$$n=1150$$y$

# Parameters to tune to the circumstances:
window <- 30
threshold <- 5

# An upper threshold ("ut") calculation based on the MAD:
library(zoo) # rollapply()
ut <- function(x) {m = median(x); median(x) + threshold * median(abs(x - m))}
z <- rollapply(zoo(y), window, ut, align="right")
z <- c(rep(z[1], window-1), z) # Use z[1] throughout the initial period
outliers <- y > z

# Graph the data, show the ut() cutoffs, and mark the outliers:
plot(x, y, type="l", lwd=2, col="#E00000", ylim=c(0, 20000))
lines(x, z, col="Gray")
points(x[outliers], y[outliers], pch=19)

Aplicado a un conjunto de datos como la curva roja ilustrada en la pregunta, produce este resultado:

Los datos se muestran en rojo, la ventana de 30 días de mediana + 5 * umbrales MAD en gris, y los valores atípicos, que son simplemente esos valores de datos por encima de la curva gris, en negro.

(El umbral solo se puede calcular comenzando al final de la ventana inicial. Para todos los datos dentro de esta ventana inicial, se usa el primer umbral: es por eso que la curva gris es plana entre x = 0 yx = 30).

Los efectos de cambiar los parámetros son (a) aumentar el valor de windowtenderá a suavizar la curva gris y (b) aumentar thresholdaumentará la curva gris. Sabiendo esto, uno puede tomar un segmento inicial de los datos e identificar rápidamente los valores de los parámetros que mejor separan los picos periféricos del resto de los datos. Aplique estos valores de parámetros para verificar el resto de los datos. Si una gráfica muestra que el método está empeorando con el tiempo, eso significa que la naturaleza de los datos está cambiando y los parámetros pueden necesitar un nuevo ajuste.

Observe lo poco que supone este método sobre los datos: no tienen que distribuirse normalmente; no necesitan exhibir ninguna periodicidad; ni siquiera tienen que ser no negativos. Todo lo que supone es que los datos se comportan de manera razonablemente similar a lo largo del tiempo y que los picos periféricos son visiblemente más altos que el resto de los datos.

Si a alguien le gustaría experimentar (o comparar alguna otra solución con la que se ofrece aquí), aquí está el código que utilicé para producir datos como los que se muestran en la pregunta.

n.length <- 1150
cycle.a <- 11
cycle.b <- 365/12
amp.a <- 800
amp.b <- 8000

set.seed(17)
x <- 1:n.length
baseline <- (1/2) * amp.a * (1 + sin(x * 2*pi / cycle.a)) * rgamma(n.length, 40, scale=1/40)
peaks <- rbinom(n.length, 1,  exp(2*(-1 + sin(((1 + x/2)^(1/5) / (1 + n.length/2)^(1/5))*x * 2*pi / cycle.b))*cycle.b))
y <- peaks * rgamma(n.length, 20, scale=amp.b/20) + baseline

Si le preocupan las suposiciones con un enfoque particular, un enfoque es capacitar a varios alumnos en diferentes señales, luego usar métodos de conjunto y agregar los "votos" de sus alumnos para hacer una clasificación atípica.

Por cierto, esto podría valer la pena leerlo o leerlo, ya que hace referencia a algunos enfoques del problema.

• Detección de valores atípicos en línea a través de flujos de datos

Supongo que el modelo sofisticado de series temporales no funcionará para usted debido al tiempo que lleva detectar los valores atípicos utilizando esta metodología. Por lo tanto, aquí hay una solución alternativa:

1. Primero establezca una línea de base de patrones de tráfico 'normales' para un año basado en el análisis manual de datos históricos que tenga en cuenta la hora del día, día de la semana vs fin de semana, mes del año, etc.

2. Use esta línea base junto con algún mecanismo simple (por ejemplo, promedio móvil sugerido por Carlos) para detectar valores atípicos.

Es posible que también desee revisar la literatura de control de procesos estadísticos para algunas ideas.

Ajuste estacionalmente los datos para que un día normal parezca más plano. Puede tomar la muestra de hoy a las 5:00 pm y restar o dividir el promedio de los 30 días anteriores a las 5:00 pm. Luego mire más allá de N desviaciones estándar (medidas usando datos preajustados) para valores atípicos. Esto podría hacerse por separado para las "estaciones" semanales y diarias.

Una alternativa al enfoque esbozado por Rob Hyndman sería utilizar el pronóstico de Holt-Winters . Las bandas de confianza derivadas de Holt-Winters se pueden usar para detectar valores atípicos. Aquí hay un documento que describe cómo usar Holt-Winters para "Detección de comportamiento aberrante en series de tiempo para monitoreo de red". Una implementación para RRDTool se puede encontrar aquí .

El análisis espectral detecta la periodicidad en series de tiempo estacionarias. El enfoque de dominio de frecuencia basado en la estimación de densidad espectral es un enfoque que recomendaría como su primer paso.

Si durante ciertos períodos la irregularidad significa un pico mucho más alto que el típico para ese período, entonces la serie con tales irregularidades no sería estacionaria y el análisis espectral no sería apropiado. Pero suponiendo que haya identificado el período que tiene las irregularidades, debería poder determinar aproximadamente cuál sería la altura máxima normal y luego establecer un umbral en algún nivel por encima de ese promedio para designar los casos irregulares.

Como se trata de datos de series temporales, un filtro exponencial simple http://en.wikipedia.org/wiki/Exponential_smoothing suavizará los datos. Es un filtro muy bueno ya que no necesita acumular puntos de datos antiguos. Comparar todos los recién suavizada valor de datos con su unsmoothed valor. Una vez que la desviación excede un cierto umbral predefinido (dependiendo de lo que creas que es un valor atípico en tus datos), entonces tu valor atípico se puede detectar fácilmente.

En CI hará lo siguiente para una muestra de 16 bits en tiempo real (creo que esto se encuentra en algún lugar aquí <Explicación - https://dsp.stackexchange.com/questions/378/what-is-the-best-first-order -iir-aproximación-a-un-filtro-promedio-móvil >)

#define BITS2 2     //< This is roughly = log2( 1 / alpha ), depending on how smooth you want your data to be

short Simple_Exp_Filter(int new_sample) 
{static int filtered_sample = 0;
long local_sample = sample << 16; /*We assume it is a 16 bit sample */
filtered_sample += (local_sample - filtered_sample) >> BITS2;   
return (short) ((filtered_sample+0x8000) >> 16); //< Round by adding .5 and truncating.   
}


int main()
{
newly_arrived = function_receive_new_sample();
filtered_sample = Simple_Exp_Filter(newly_arrived);
if (abs(newly_arrived - filtered_sample)/newly_arrived > THRESHOLD)
    {
    //AN OUTLIER HAS BEEN FOUND
    }
 return 0;   
}

Puede usar la desviación estándar de las últimas N mediciones (debe elegir un N adecuado). Un buen puntaje de anomalía sería cuántas desviaciones estándar es una medición del promedio móvil.

lo que hago es agrupar las mediciones por hora y día de la semana y comparar las desviaciones estándar de eso. Todavía no se corrige para cosas como vacaciones y temporada de verano / invierno, pero es correcto la mayor parte del tiempo.

La desventaja es que realmente necesita recopilar aproximadamente un año de datos para tener suficiente para que stddev comience a tener sentido.

Sugiero el siguiente esquema, que debería ser implementable en un día más o menos:

Formación

• Recolecta tantas muestras como puedas guardar en la memoria
• Eliminar valores atípicos obvios utilizando la desviación estándar para cada atributo
• Calcule y almacene la matriz de correlación y también la media de cada atributo
• Calcule y almacene las distancias de Mahalanobis de todas sus muestras

Cálculo de "valores atípicos":

Para la muestra única de la que desea saber su "atípica":

• Recupere las medias, la matriz de covarianza y las distancias de Mahalanobis del entrenamiento
• Calcule la distancia de Mahalanobis "d" para su muestra
• Devuelve el percentil en el que cae "d" (usando las distancias de Mahalanobis desde el entrenamiento)

Ese será su puntaje atípico: 100% es un valor atípico extremo.

Para el caso en que uno tiene que calcular los valores atípicos rápidamente, uno podría usar la idea de Rob Hyndman y Mahito Sugiyama ( https://github.com/BorgwardtLab/sampling-outlier-detection , library (spoutlier), function qsp) para calcular los valores atípicos de la siguiente manera:

library(spoutlier)
rapidtsoutliers <- function(x,plot=FALSE,seed=123)
{
    set.seed(seed)
    x <- as.numeric(x)
    tt <- 1:length(x)
    qspscore <- qsp(x)
    limit <- quantile(qspscore,prob=c(0.95))
    score <- pmax((qspscore - limit),0)
    if(plot)
    {
        plot(x,type="l")
        x2 <- ts(rep(NA,length(x)))
        x2[score>0] <- x[score>0]
        tsp(x2) <- tsp(x)
        points(x2,pch=19,col="red")
        return(invisible(score))
    }
    else
        return(score)
}

La detección de anomalías requiere la construcción de una ecuación que describa la expectativa. La detección de intervención está disponible en un entorno no causal y causal. Si uno tiene una serie de predictores como el precio, las cosas pueden complicarse un poco. Otras respuestas aquí no parecen tener en cuenta la causa asignable atribuible a series de predictores especificadas por el usuario como el precio y, por lo tanto, pueden ser defectuosas. La cantidad vendida puede depender del precio, quizás precios anteriores y tal vez la cantidad vendida en el pasado. La base para la detección de anomalías (pulsos, pulsos estacionales, cambios de nivel y tendencias de hora local) se encuentra en https://pdfs.semanticscholar.org/09c4/ba8dd3cc88289caf18d71e8985bdd11ad21c.pdf