Algoritmo simples para detecção on-line de uma série temporal genérica

Estou trabalhando com uma grande quantidade de séries temporais. Essas séries temporais são basicamente medidas de rede que ocorrem a cada 10 minutos, e algumas são periódicas (ou seja, a largura de banda), enquanto outras não são (ou seja, a quantidade de tráfego de roteamento).

Eu gostaria de um algoritmo simples para fazer uma "detecção de outlier" on-line. Basicamente, quero manter na memória (ou em disco) todos os dados históricos de cada série temporal e quero detectar qualquer erro externo em um cenário ativo (cada vez que uma nova amostra é capturada). Qual é a melhor maneira de alcançar esses resultados?

Atualmente, estou usando uma média móvel para remover algum ruído, mas e depois? Coisas simples, como desvio padrão, louco, ... contra todo o conjunto de dados não funcionam bem (não posso assumir que as séries temporais sejam estacionárias) e gostaria de algo mais "preciso", idealmente uma caixa preta como:

double outlier_detection (vetor duplo *, valor duplo);

em que vetor é a matriz double que contém os dados históricos e o valor de retorno é a pontuação de anomalia para o novo "valor" da amostra.

Aqui está uma função R simples que encontrará discrepantes de séries temporais (e, opcionalmente, mostra-os em um gráfico). Ele lidará com séries temporais sazonais e não sazonais. A idéia básica é encontrar estimativas robustas da tendência e dos componentes sazonais e subtraí-los. Em seguida, encontre valores discrepantes nos resíduos. O teste para outliers residuais é o mesmo que para o boxplot padrão - pontos superiores a 1,5IQR acima ou abaixo dos quartis superior e inferior são considerados outliers. O número de IQRs acima / abaixo desses limites é retornado como uma "pontuação" externa. Portanto, a pontuação pode ser qualquer número positivo e será zero para não-discrepantes.

Sei que você não está implementando isso em R, mas geralmente acho uma função R um bom lugar para começar. A tarefa é traduzir isso para qualquer idioma necessário.

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)
}

Uma boa solução terá vários ingredientes, incluindo:

• Use uma janela resistente e em movimento suave para remover a não estacionaridade.

• Expresse novamente os dados originais para que os resíduos com relação ao liso sejam distribuídos aproximadamente simetricamente. Dada a natureza dos seus dados, é provável que suas raízes quadradas ou logaritmos gerem resíduos simétricos.

• Aplique os métodos do gráfico de controle, ou pelo menos o pensamento do gráfico de controle, aos resíduos.

Quanto ao último, o pensamento da carta de controle mostra que os limites "convencionais", como 2 SD ou 1,5 vezes o IQR além dos quartis, funcionam mal porque desencadeiam muitos sinais falsos fora de controle. As pessoas costumam usar 3 DP no trabalho da carta de controle, sendo 2,5 (ou até 3) vezes o IQR além dos quartis seria um bom ponto de partida.

Descrevi mais ou menos a natureza da solução de Rob Hyndman, acrescentando a ela dois pontos principais: a necessidade potencial de reexprimir os dados e a sabedoria de ser mais conservador ao sinalizar valores extremos. Não tenho certeza se Loess é bom para um detector on-line, porque não funciona bem nos terminais. Em vez disso, você pode usar algo tão simples quanto um filtro mediano em movimento (como na suavização resistente de Tukey). Se os outliers não aparecerem em rajadas, você poderá usar uma janela estreita (5 pontos de dados, talvez, que serão quebrados apenas com uma explosão de 3 ou mais outliers em um grupo de 5).

Depois de realizar a análise para determinar uma boa reexpressão dos dados, é improvável que você precise alterar a reexpressão. Portanto, o seu detector on-line realmente precisa apenas fazer referência aos valores mais recentes (a janela mais recente) porque não usará os dados anteriores. Se você tem realmente séries temporais muito longas, analise a autocorrelação e a sazonalidade (como flutuações diárias ou semanais recorrentes) para melhorar o procedimento.

(Esta resposta respondeu a uma pergunta duplicada (agora fechada) em Detectando eventos pendentes , que apresentou alguns dados em forma gráfica.)

A detecção de outlier depende da natureza dos dados e do que você deseja assumir sobre eles. Métodos de uso geral confiam em estatísticas robustas. O espírito dessa abordagem é caracterizar a maior parte dos dados de uma maneira que não seja influenciada por discrepantes e, em seguida, apontar para quaisquer valores individuais que não se encaixam nessa caracterização.

Por se tratar de uma série temporal, ela acrescenta a complicação de precisar (re) detectar outliers continuamente. Se isso for feito à medida que a série se desenrola, temos permissão apenas para usar dados mais antigos para a detecção, não dados futuros! Além disso, como proteção contra muitos testes repetidos, gostaríamos de usar um método com uma taxa de falsos positivos muito baixa.

Essas considerações sugerem a execução de um teste externo simples e robusto da janela móvel sobre os dados . Existem muitas possibilidades, mas uma simples, fácil de entender e facilmente implementada é baseada em um MAD: desvio absoluto médio da mediana. Essa é uma medida de variação fortemente robusta dentro dos dados, semelhante a um desvio padrão. Um pico externo seria vários MADs ou mais maior que a mediana.

Ainda há alguns ajustes a serem feitos : quanto de um desvio em relação à maior parte dos dados deve ser considerado distante e a que distância do tempo se deve olhar? Vamos deixá-los como parâmetros para experimentação. Aqui está uma Rimplementação aplicada aos dados (com para emular os dados) com os valores correspondentes :n = 1150 $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 um conjunto de dados como a curva vermelha ilustrada na pergunta, ele produz este resultado:

Os dados são mostrados em vermelho, a janela de 30 dias dos limiares medianos + 5 * MAD em cinza e os valores extremos - que são simplesmente aqueles valores acima da curva cinza - em preto.

(O limite só pode ser calculado a partir do final da janela inicial. Para todos os dados nessa janela inicial, o primeiro limite é usado: é por isso que a curva cinza é plana entre x = 0 ex = 30).

Os efeitos da alteração dos parâmetros são (a) aumentar o valor de windowtenderá a suavizar a curva cinza e (b) aumentar thresholdaumentará a curva cinza. Sabendo disso, pode-se pegar um segmento inicial dos dados e identificar rapidamente os valores dos parâmetros que melhor segregam os picos externos do restante dos dados. Aplique esses valores de parâmetro para verificar o restante dos dados. Se um gráfico mostrar que o método está piorando com o tempo, isso significa que a natureza dos dados está mudando e os parâmetros podem precisar de um novo ajuste.

Observe o quão pouco esse método assume sobre os dados: eles não precisam ser normalmente distribuídos; eles não precisam exibir periodicidade; eles nem precisam ser negativos. Tudo pressupõe que os dados se comportam de maneiras razoavelmente semelhantes ao longo do tempo e que os picos externos são visivelmente mais altos que o restante dos dados.

Se alguém quiser experimentar (ou comparar alguma outra solução à oferecida aqui), eis o código que usei para produzir dados como os mostrados na pergunta.

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

Se você está preocupado com suposições com qualquer abordagem específica, uma abordagem é treinar vários alunos em sinais diferentes, use métodos de conjunto e agregue os "votos" de seus alunos para fazer a classificação externa.

BTW, isso pode valer a pena ser lido ou lido, pois faz referência a algumas abordagens do problema.

• Detecção de outlier on-line sobre fluxos de dados

Suponho que o modelo sofisticado de séries temporais não funcione para você por causa do tempo necessário para detectar discrepâncias usando essa metodologia. Portanto, aqui está uma solução alternativa:

1. Primeiro, estabeleça padrões de tráfego 'normais' de linha de base para um ano com base na análise manual de dados históricos, que representam a hora do dia, o dia da semana versus o final de semana, o mês do ano etc.

2. Use esta linha de base juntamente com algum mecanismo simples (por exemplo, média móvel sugerida por Carlos) para detectar valores extremos.

Você também pode revisar a literatura de controle estatístico de processos para obter algumas idéias.

Ajuste os dados sazonalmente para que um dia normal pareça mais plano. Você pode coletar a amostra das 17:00 de hoje e subtrair ou dividir a média dos 30 dias anteriores às 17:00. Em seguida, procure N desvios padrão (medidos usando dados pré-ajustados) em busca de outliers. Isso pode ser feito separadamente para as "estações" semanais e diárias.

Uma alternativa à abordagem descrita por Rob Hyndman seria usar a Holt-Winters Forecasting . As faixas de confiança derivadas de Holt-Winters podem ser usadas para detectar discrepâncias. Aqui está um artigo que descreve como usar Holt-Winters para "Detecção de comportamento aberrante em séries temporais para monitoramento de rede". Uma implementação para o RRDTool pode ser encontrada aqui .

A análise espectral detecta periodicidade em séries temporais estacionárias. A abordagem no domínio da frequência baseada na estimativa da densidade espectral é uma abordagem que eu recomendaria como seu primeiro passo.

Se durante certos períodos a irregularidade significa um pico muito mais alto do que o normal para esse período, as séries com tais irregularidades não seriam estacionárias e a anlise espectral não seria apropriada. Mas, supondo que você tenha identificado o período que apresenta as irregularidades, você poderá determinar aproximadamente qual seria a altura normal do pico e, em seguida, poderá definir um limite em algum nível acima dessa média para designar os casos irregulares.

Como se trata de dados de séries temporais, um simples filtro exponencial http://en.wikipedia.org/wiki/Exponential_smoothing suavizará os dados. É um filtro muito bom, pois você não precisa acumular pontos de dados antigos. Compare cada valor de dados recém- suavizado com seu valor não suavizado . Depois que o desvio exceder um determinado limite predefinido (dependendo do que você acredita que seja um erro de informação), ele pode ser facilmente detectado.

No IC, faça o seguinte para uma amostra de 16 bits em tempo real (acredito que isso seja encontrado em algum lugar aqui <Explicação - https://dsp.stackexchange.com/questions/378/what-is-the-best-first-order -iir-aproximação-para-um-filtro-de-média-móvel >)

#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;   
}

Você pode usar o desvio padrão das últimas N medidas (você precisa escolher um N adequado). Uma boa pontuação de anomalia seria quantos desvios padrão uma medida é da média móvel.

o que faço é agrupar as medidas por hora e dia da semana e comparar os desvios padrão disso. Ainda não corrige coisas como feriados e sazonalidade de verão / inverno, mas é correta na maioria das vezes.

A desvantagem é que você realmente precisa coletar um ano ou mais de dados para ter o suficiente para que o stddev comece a fazer sentido.

Sugiro o esquema abaixo, que deve ser implementado em um dia ou mais:

Treinamento

• Colete o maior número possível de amostras na memória
• Remova discrepantes óbvios usando o desvio padrão para cada atributo
• Calcular e armazenar a matriz de correlação e também a média de cada atributo
• Calcule e armazene as distâncias de Mahalanobis de todas as suas amostras

Cálculo de "outicidade":

Para a amostra única da qual você deseja conhecer sua "aparência extravagante":

• Recuperar as médias, matriz de covariância e distância s de Mahalanobis do treinamento
• Calcule a distância de Mahalanobis "d" para sua amostra
• Retorne o percentil em que "d" cai (usando as distâncias de Mahalanobis do treinamento)

Essa será sua pontuação outlier: 100% é uma outlier extrema.

No caso em que é preciso calcular os valores discrepantes rapidamente, pode-se usar a ideia de Rob Hyndman e Mahito Sugiyama ( https://github.com/BorgwardtLab/sampling-outlier-detection , library (spoutlier), função qsp) para calcular os valores extremos da seguinte forma:

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)
}

A detecção de anomalias requer a construção de uma equação que descreva a expectativa. A detecção de intervenção está disponível em um cenário não causal e causal. Se alguém tem uma série de preditores como preço, as coisas podem ficar um pouco complicadas. Outras respostas aqui não parecem levar em consideração a causa atribuível atribuível a séries de preditores especificadas pelo usuário, como preço, e, portanto, podem ter falhas. A quantidade vendida pode depender do preço, talvez preços anteriores e talvez quantidade vendida no passado. A base para a detecção de anomalias (pulsos, pulsos sazonais, mudanças de nível e tendências da hora local) encontra-se em https://pdfs.semanticscholar.org/09c4/ba8dd3cc88289caf18d71e8985bdd11ad21c.pdf