Simple algorithm for online outlier detection of a generic time series

Question

I am working with a large amount of time series. These time series are basically network measurements coming every 10 minutes, and some of them are periodic (i.e. the bandwidth), while some other aren't (i.e. the amount of routing traffic).

I would like a simple algorithm for doing an online "outlier detection". Basically, I want to keep in memory (or on disk) the whole historical data for each time series, and I want to detect any outlier in a live scenario (each time a new sample is captured). What is the best way to achieve these results?

I'm currently using a moving average in order to remove some noise, but then what next? Simple things like standard deviation, mad, ... against the whole data set doesn't work well (I can't assume the time series are stationary), and I would like something more "accurate", ideally a black box like:

double outlier_detection(double* vector, double value);

where vector is the array of double containing the historical data, and the return value is the anomaly score for the new sample "value" .

Answer 1

Here is a simple R function that will find time series outliers (and optionally show them in a plot). It will handle seasonal and non-seasonal time series. The basic idea is to find robust estimates of the trend and seasonal components and subtract them. Then find outliers in the residuals. The test for residual outliers is the same as for the standard boxplot -- points greater than 1.5IQR above or below the upper and lower quartiles are assumed outliers. The number of IQRs above/below these thresholds is returned as an outlier "score". So the score can be any positive number, and will be zero for non-outliers.

I realise you are not implementing this in R, but I often find an R function a good place to start. Then the task is to translate this into whatever language is required.

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

Answer 2

A good solution will have several ingredients, including:

• Use a resistant, moving window smooth to remove nonstationarity.

• Re-express the original data so that the residuals with respect to the smooth are approximately symmetrically distributed. Given the nature of your data, it's likely that their square roots or logarithms would give symmetric residuals.

• Apply control chart methods, or at least control chart thinking, to the residuals.

As far as that last one goes, control chart thinking shows that "conventional" thresholds like 2 SD or 1.5 times the IQR beyond the quartiles work poorly because they trigger too many false out-of-control signals. People usually use 3 SD in control chart work, whence 2.5 (or even 3) times the IQR beyond the quartiles would be a good starting point.

I have more or less outlined the nature of Rob Hyndman's solution while adding to it two major points: the potential need to re-express the data and the wisdom of being more conservative in signaling an outlier. I'm not sure that Loess is good for an online detector, though, because it doesn't work well at the endpoints. You might instead use something as simple as a moving median filter (as in Tukey's resistant smoothing). If outliers don't come in bursts, you can use a narrow window (5 data points, perhaps, which will break down only with a burst of 3 or more outliers within a group of 5).

Once you have performed the analysis to determine a good re-expression of the data, it's unlikely you'll need to change the re-expression. Therefore, your online detector really only needs to reference the most recent values (the latest window) because it won't use the earlier data at all. If you have really long time series you could go further to analyze autocorrelation and seasonality (such as recurring daily or weekly fluctuations) to improve the procedure.

Answer 3

(This answer responded to a duplicate (now closed) question at Detecting outstanding events, which presented some data in graphical form.)

---

Outlier detection depends on the nature of the data and on what you are willing to assume about them. General-purpose methods rely on robust statistics. The spirit of this approach is to characterize the bulk of the data in a way that is not influenced by any outliers and then point to any individual values that do not fit within that characterization.

Because this is a time series, it adds the complication of needing to (re)detect outliers on an ongoing basis. If this is to be done as the series unfolds, then we are allowed only to use older data for the detection, not future data! Moreover, as protection against the many repeated tests, we would want to use a method that has a very low false positive rate.

These considerations suggest running a simple, robust moving window outlier test over the data. There are many possibilities, but one simple, easily understood and easily implemented one is based on a running MAD: median absolute deviation from the median. This is a strongly robust measure of variation within the data, akin to a standard deviation. An outlying peak would be several MADs or more greater than the median.

There is still some tuning to be done: how much of a deviation from the bulk of the data should be considered outlying and how far back in time should one look? Let's leave these as parameters for experimentation. Here's an R implementation applied to data $x = (1,2,\ldots,n)$ (with $n=1150$ to emulate the data) with corresponding values $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)

Applied to a dataset like the red curve illustrated in the question, it produces this result:

The data are shown in red, the 30-day window of median+5*MAD thresholds in gray, and the outliers--which are simply those data values above the gray curve--in black.

(The threshold can only be computed beginning at the end of the initial window. For all data within this initial window, the first threshold is used: that's why the gray curve is flat between x=0 and x=30.)

The effects of changing the parameters are (a) increasing the value of window will tend to smooth out the gray curve and (b) increasing threshold will raise the gray curve. Knowing this, one can take an initial segment of the data and quickly identify values of the parameters that best segregate the outlying peaks from the rest of the data. Apply these parameter values to checking the rest of the data. If a plot shows the method is worsening over time, that means the nature of the data are changing and the parameters might need re-tuning.

Notice how little this method assumes about the data: they do not have to be normally distributed; they do not need to exhibit any periodicity; they don't even have to be non-negative. All it assumes is that the data behave in reasonably similar ways over time and that the outlying peaks are visibly higher than the rest of the data.

---

If anyone would like to experiment (or compare some other solution to the one offered here), here is the code I used to produce data like those shown in the question.

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

Answer 4

If you're worried about assumptions with any particular approach, one approach is to train a number of learners on different signals, then use ensemble methods and aggregate over the "votes" from your learners to make the outlier classification.

BTW, this might be worth reading or skimming since it references a few approaches to the problem.

• Online outlier detection over data streams

Answer 5

I am guessing sophisticated time series model will not work for you because of the time it takes to detect outliers using this methodology. Therefore, here is a workaround:

1. First establish a baseline 'normal' traffic patterns for a year based on manual analysis of historical data which accounts for time of the day, weekday vs weekend, month of the year etc.

2. Use this baseline along with some simple mechanism (e.g., moving average suggested by Carlos) to detect outliers.

You may also want to review the statistical process control literature for some ideas.

Answer 6

Seasonally adjust the data such that a normal day looks closer to flat. You could take today's 5:00pm sample and subtract or divide out the average of the previous 30 days at 5:00pm. Then look past N standard deviations (measured using pre-adjusted data) for outliers. This could be done separately for weekly and daily "seasons."

Answer 7

An alternative to the approach outlined by Rob Hyndman would be to use Holt-Winters Forecasting . The confidence bands derived from Holt-Winters can be used to detect outliers. Here is a paper that describes how to use Holt-Winters for "Aberrant Behavior Detection in Time Series for Network Monitoring". An implementation for RRDTool can be found here.

Answer 8

Spectral analysis detects periodicity in stationary time series. The frequency domain approach based on spectral density estimation is an approach I would recommend as your first step.

If for certain periods irregularity means a much higher peak than is typical for that period then the series with such irregularities would not be stationary and spectral anlsysis would not be appropriate. But assuming you have identified the period that has the irregularities you should be able to determine approximately what the normal peak height would be and then can set a threshold at some level above that average to designate the irregular cases.

Answer 9

Since it is a time series data, a simple exponential filter http://en.wikipedia.org/wiki/Exponential_smoothing will smoothen the data. It is a very good filter since you don't need to accumulate old data points. Compare every newly smoothed data value with its unsmoothed value. Once the deviation exceeds a certain predefined threshold (depending on what you believe an outlier in your data is), then your outlier can be easily detected.

In C I will do the following for a real-time 16 bit sample (I believe this is found somewhere here < Explanation - https://dsp.stackexchange.com/questions/378/what-is-the-best-first-order-iir-approximation-to-a-moving-average-filter>)

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

Answer 10

You could use the standard deviation of the last N measurements (you have to pick a suitable N). A good anomaly score would be how many standard deviations a measurement is from the moving average.

Answer 11

what I do is group the measurements by hour and day of week and compare standard deviations of that. Still doesn't correct for things like holidays and summer/winter seasonality but its correct most of the time.

The downside is that you really need to collect a year or so of data to have enough so that stddev starts making sense.

Answer 12

anomaly detection requires the construction of an equation which describes the expectation. Intervention Detection is available in both a non-causal and causal setting . If one has a predictor series like price then things can get a little complicated. Other responses here don't seem to take into account assignable cause attributable to user specified predictor series like price and thus might be flawed. Quantity sold may well depend on price , perhaps previous prices and perhaps quantity sold in the past. The basis for the anomaly detection ( pulses,seasonal pulses, level shifts and local time trends ) is found in https://pdfs.semanticscholar.org/09c4/ba8dd3cc88289caf18d71e8985bdd11ad21c.pdf

Answer 13

For the case where one has to compute the outliers quickly, one could use the idea of Rob Hyndman and Mahito Sugiyama ( https://github.com/BorgwardtLab/sampling-outlier-detection , library(spoutlier), function qsp ) to compute the outliers as follows:

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