DATA SCIENCE REPORT SERIES

ANOMALY DETECTION AND OUTLIER ANALYSIS

Youssouph Cissokho1 , Soufiane Fadel1 , Richard Millson1 , Razieh Pourhasan1 , Patrick Boily1,2,3

Abstract
With the advent of automatic data collection, it is now possible to store and process large troves of data.
There are technical issues associated to massive data sets, such as the speed and efficiency of analytical
methods, but there are also problems related to the detection of anomalous observations and the analysis
of outliers.

Extreme and irregular values behave very differently from the majority of observations. For instance, they
can represent criminal attacks, fraud attempts, targeted attacks, or data collection errors. As a result,
anomaly detection and outlier analysis play a crucial role in cybersecurity, quality control, etc. [1,3,4]. The
(potentially) heavy human price and technical consequences related to the presence of such observations
go a long way towards explaining why the topic has attracted attention in recent years.

This report contains a review of various detection methods, with particular attention paid to both supervised
and unsupervised methods, as well as an application to time series data and a project suggestion
(comparative analysis of various algorithms applied to 5 real-world datasets).

Keywords
Anomaly detection, outlier analysis.
Funding Acknowledgement
This report was funded by a University of Ottawa grant to develop teaching material in French (2019-2020). It was subsequently
translated into English to produce the current document.
1
Department of Mathematics and Statistics, University of Ottawa, Ottawa
2
Data Action Lab, Ottawa
3
Idlewyld Analytics and Consulting Services, Wakefield, Canada
Email: pboily@[Link]

Contents 4 Anomalies in High-Dimensional Datasets 16

4.1 Definitions and Challenges . . . . . . . . . . . . . . . . . . . 16
1 Introduction 2 4.2 Projection-Based Methods . . . . . . . . . . . . . . . . . . . 16
1.1 Basic Notions and Overview . . . . . . . . . . . . . . . . . . . 2 4.3 Ensembles Methods . . . . . . . . . . . . . . . . . . . . . . . 18
1.2 Anomaly Detection as a Statistical Learning Problem . . . . . 4 4.4 Subspace Methods . . . . . . . . . . . . . . . . . . . . . . . 20
1.3 Suggested References . . . . . . . . . . . . . . . . . . . . . . 9
5 Applications to Time Series 20
1.4 Structure and Organization . . . . . . . . . . . . . . . . . . . . 9
5.1 Outliers and Anomalies in Time Series . . . . . . . . . . . . 20
2 Quantitative Methods of Anomaly Detection 10 5.2 R Package: tsoutliers . . . . . . . . . . . . . . . . . . 21
2.1 Distance-Based Methods . . . . . . . . . . . . . . . . . . . . 10 5.3 R Package: anomalize . . . . . . . . . . . . . . . . . . . 21
2.2 Density-Based Methods . . . . . . . . . . . . . . . . . . . . . 11 5.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3 Qualitative Methods of Anomaly Detection 14 6 Project 23
3.1 Definitions and Challenges . . . . . . . . . . . . . . . . . . . 14
3.2 Review of Two Methods . . . . . . . . . . . . . . . . . . . . . 15
DATA SCIENCE REPORT SERIES ANOMALY DETECTION AND OUTLIER ANALYSIS

1. Introduction shops, for instance, a purchase which is recorded at 3AM

local time does not necessarily raise a red flag anymore.
Isaac Asimov, the prolific American author, once wrote that When anomalies are actually associated to malicious
The most exciting phrase to hear [...], the one activities, they are more than often disguised in order to
that heralds the most discoveries, is not “Eu- blend in with normal observations, which obviously com-
reka!” but “That’s funny...”. plicates the detection process.

However, anomalous observations are not only harbingers Numerous methods exist to identify anomalous observa-
of great scientific discoveries – unexpected observations tions; none of them are foolproof and judgement must
can spoil analyses or be indicative of the presence of issues be used. Methods that employ graphical aids (such as box-
related to data collection or data processing. plots, scatterplots, scatterplot matrices, and 2D tours) to
Either way, it becomes imperative for decision-makers identify outliers are particularly easy to implement, but
and analysts to establish anomaly detection protocols, and a low-dimensional setting is usually required for ease of
to identify strategies to deal with such observations. interpretability.
Analytical methods also exist (using Cooke’s or Maha-
1.1 Basic Notions and Overview lanobis’ distances, say), but in general some additional level
Outlying observations are data points which are atypi- of analysis must be performed, especially when trying to
cal in comparison to the unit’s remaining features (within- identify influential points (cf. leverage).
unit), or in comparison to the measurements for other units
(between-units), or as part of a collective subset of observa- With small datasets, anomaly detection can be conducted
tions. Outliers are thus observations which are dissimilar on a case-by-case basis, but with large datasets, the tempta-
to other cases or which contradict known dependencies tion to use automated detection/removal is strong – care
or rules.1 must be exercised before the analyst decides to go down
Observations could be anomalous in one context, but that route.3
not in another. Consider, for instance, an adult male who In the early stages of anomaly detection, simple data
is 6-foot tall. Such a man would fall in the 86th percentile analyses (such as descriptive statistics, 1- and 2-way tables,
among Canadian males [23], which, while on the tall side, and traditional visualisations) may be performed to help
is not unusual; in Bolivia, however, the same man would identify anomalous observations, or to obtain insights about
land in the 99.9th percentile [23], which would mark him the data, which could eventually lead to modifications of
as extremely tall and quite dissimilar to the rest of the pop- the analysis plan.
ulation.2
Outlier Tests
A common mistake that analysts make when dealing with How are outliers actually detected? Most methods come
outlying observations is to remove them from the dataset in one of two flavours: supervised and unsupervised (we
without carefully studying whether they are influential will discuss those in detail in later sections).
data points, that is, observations whose absence leads to Supervised methods use a historical record of labeled
markedly different analysis results. (that is to say, previously identified) anomalous observa-
When influential observations are identified, remedial tions to build a predictive classification or regression
measures (such as data transformation strategies) may need model which estimates the probability that a unit is anoma-
to be applied to minimize any undue effect. Note that lous; domain expertise is required to tag the data. Since
outliers may be influential, and influential data points may anomalies are typically infrequent, these models often also
be outliers, but the conditions are neither necessary nor have to accommodate the rare occurrence problem.4
sufficient. Unsupervised methods, on the other hand, use no pre-
viously labeled information or data, and try to determine if
Anomaly Detection an observation is an outlying one solely by comparing its
By definition, anomalies are infrequent and typically shrouded behaviour to that of the other observations.
in uncertainty due to their relatively low numbers, which
3
makes it difficult to differentiate them from banal noise or This stems partly from the fact that once the “anomalous” observations
have been removed from the dataset, previously “regular” observations
data collection errors.
can become anomalous in turn in the smaller dataset; it is not clear when
Furthermore, the boundary between normal and de- that runaway train will stop.
viant observations is usually fuzzy; with the advent of e- 4
Supervised models are built to minimize a cost function; in default
settings, it is often the case that the mis-classification cost is assumed to
1
Outlying observations may be anomalous along any of the individual be symmetrical, which can lead to technically correct but useless solutions.
variables, or in combination. For instance, the vast majority (99.999+%) of air passengers emphatically
2
Anomaly detection points towards interesting questions for analysts do not bring weapons with them on flights; a model that predicts that no
and subject matter experts: in this case, why is there such a large discrep- passenger is attempting to smuggle a weapon on board a flight would be
ancy in the two populations? 99.999+% accurate, but it would miss the point completely.

[Link], [Link], [Link], [Link], [Link] (2020) Page 2 of 25

DATA SCIENCE REPORT SERIES ANOMALY DETECTION AND OUTLIER ANALYSIS

– the Tietjen-Moore test, which is used to find a

specific number of outliers;
– the generalized extreme studentized deviate
test, if the number of outliers is unknown;
– the chi-square test, when outliers affect the
goodness-of-fit, as well as
– DBSCAN and other clustering-based outlier de-
tection methods.

Visual Outlier Detection

The following three (simple) examples illustrate the princi-
ples underlying visual outlier and anomaly detection.
Example 1. On a specific day, the height of several plants
Figure 1. Tukey’s boxplot test; suspected outliers are marked by in a nursery are measured. The records also show each
white disks, outliers by black disks. plant’s age (the number of weeks since the seed has been
planted).
Histograms of the data are shown in Figure 2 (age on
The following traditional methods and tests of outlier de-
the left, height on the middle).
tection fall into this category:5
Very little can be said about the data at that stage: the
Perhaps the most commonly-used test is Tukey’s box- age of the plants (controlled by the nursery staff) seems
plot test; for normally distributed data, regular ob- to be somewhat haphazard, as does the response variable
servations typically lie between the inner fences (height). A scatter plot of the data (rightmost chart in Fig-
ure 2), however, reveals that growth is strongly correlated
Q 1 − 1.5(Q 3 − Q 1 ) and Q 3 + 1.5(Q 3 − Q 1 ). with age during the early period of a plant’s life for the
observations in the dataset; points clutter around a lin-
Suspected outliers lie between the inner fences and
ear trend. One point (in yellow) is easily identified as an
their respective outer fences
outlier. There are (at least) two possibilities: either that
Q 1 − 3(Q 3 − Q 1 ) and Q 3 + 3(Q 3 − Q 1 ). measurement was botched or mis-entered in the database
(representing an invalid entry), or that one specimen has ex-
Points beyond the outer fences are identified as out- perienced unusual growth (outlier). Either way, the analyst
liers (Q 1 and Q 3 represent the data’s 1st and 3rd quar- has to investigate further.
tile, respectively; see Figure 1).
Example 2. A government department has 11 service points
The Grubbs test is another univariate test, which
in a jurisdiction. Service statistics are recorded: the monthly
takes into consideration the number of observations
average arrival rates per teller and monthly average service
in the dataset. Let x i be the value of feature X for
rates per teller for each service point are available.
the i th unit, 1 ≤ i ≤ N , let (x, s x ) be the mean and
A scatter plot of the service rate per teller ( y axis)
standard deviation of feature X , let α be the desired
against the arrival rate per teller (x axis), with linear re-
significance level, and let T (α, N ) be the critical value
gression trend, is shown in the leftmost chart in Figure 3.
of the Student t-distribution at significance α/2N .
The trend is seen to inch upwards with increasing x values.
Then, the i th unit is an outlier along feature X if
A similar chart, but with the left-most point removed
from consideration, is shown in the middle chart of Figure 3.
v
s x (N − 1) t T 2 (α, N )
|x i − x| ≥ p . The trend still slopes upward, but the fit is significantly im-
N N − 2 + T 2 (α, N )
proved, suggesting that the removed observation is unduly
Other common tests include: influential (or anomalous) – a better understanding of the
– the Dixon Q test, which is used in the experi- relationship between arrivals and services is afforded if it
mental sciences to find outliers in (extremely) is set aside.
small datasets – it is of dubious validity; Any attempt to fit that data point into the model must
– the Mahalanobis distance, which is linked to take this information into consideration. Note, however,
the leverage of an observation (a measure of that influential observations depend on the analysis that is
influence), can also be used to find multi-dimen- ultimately being conducted – a point may be influential for
sional outliers, when all relationships are linear one analysis, but not for another.
(or nearly linear); Example 3. Measurements of the length of the appendage
5 of a certain species of insect have been made on 71 individ-
Note that normality of the underlying data is an assumption for most
tests; how robust these tests are against departures from this assumption uals. Descriptive statistics have been computed; the results
depends on the situation. are shown in Figure 4.

[Link], [Link], [Link], [Link], [Link] (2020) Page 3 of 25

DATA SCIENCE REPORT SERIES ANOMALY DETECTION AND OUTLIER ANALYSIS

Figure 2. Summary visualisations for an (artificial) plant dataset: age distribution (left), height distribution (middle), height vs. age,
with linear trend (right).

Figure 3. Visualisations for an (artificial) service point dataset: trend for 11 service points (left), trend for 10 service points (middle),
influential observations (right).

Analysts who are well-versed in statistical methods might 1.2 Anomaly Detection as a Statistical Learning Problem
recognize the tell-tale signs that the distribution of ap- Fraudulent behaviour is not always easily identifiable, even
pendage lengths is likely to be asymmetrical (since the after the fact. Credit card fraudsters, for instance, will try
skewness is non-negligible) and to have a “fat” tail (due to to disguise their transactions as regular and banal, rather
the kurtosis being commensurate with the mean and the than as outlandish; to fool human observers into confusing
standard deviation, the range being so much larger than what is merely plausible with what is probable (or at least,
the interquartile range, and the maximum value being so not improbable).
much larger than the third quartile).
The mode, minimum, and first quartile values belong At its most basic level, anomaly detection is a problem in
to individuals without appendages, so there appears to be applied probability: if I denotes what is known about the
at least two sub-groups in the population (perhaps split dataset (behaviour of individual observations, behaviour of
along the lines of juveniles/adults, or males/females). The observations as a group, anomalous/normal verdict for a
maximum value has already been seen to be quite large number of similar observations, etc.), is
compared to the rest of the observations, which at first
P(obs. is anomalous | I) > P(obs. is normal | I)?
suggests that it might belong to an outlier.
The histogram of the measurements, however, shows Anomaly detection models usually assume stationarity for
that there are 3 individuals with very long appendages (see normal observations, which is to say, that the underlying
right-most chart in Figure 4): it now becomes plausible for mechanism that generates data does not change in a sub-
these anomalous entries to belong to individuals from a stantial manner over time, or, if it does, that its rate of
different species altogether who were erroneously added change or cyclicity is known.
to the dataset. This does not, of course, constitute a proof
of such an error, but it raises the possibility, which is often For time series data, this means that it may be necessary to
the best that an analyst can do in the absence of subject first perform trend and seasonality extraction.
matter expertise.
Example 4. Supply chains play a crucial role in the trans-
This traditional approach to anomaly detection fails for portation of goods from one part of the world to another. As
high-dimensional datasets, however, and a fundamentally the saying goes, “a given chain is only as strong as its weak-
different approach is advocated. est link” – in a multi-modal context, comparing the various

[Link], [Link], [Link], [Link], [Link] (2020) Page 4 of 25

DATA SCIENCE REPORT SERIES ANOMALY DETECTION AND OUTLIER ANALYSIS

Figure 4. Summary and visualisation for an (artificial) appendage length dataset: descriptive statistics (left), appendage length
distribution (right).

transportation segments is far from an obvious endeavour: The supply chain under investigation has Shanghai as the
if shipments departing Shanghai in February 2013 took two point of origin of shipments, with Toronto as the final des-
more days, on average, to arrive in Vancouver than those tination; the containers enter the country either through
departing in July 2017, can it be said with any certainty Vancouver or Prince Rupert. Containers leave their point of
that the shipping process has improved in the intervening origin by boat, arrive and dwell in either of the two ports
years? Are February departures always slower to cross the before reaching their final destination by rail.
Pacific Ocean? For each of the three segments (Marine Transit, Port
The seasonal variability of performance is relevant to Dwell, Rail Transit), the data consists of the monthly em-
supply chain monitoring; the ability to quantify and account pirical distribution of transit times, built from sub-samples
for the severity of its impact on the data is thus of great (assumed to be randomly selected and fully representative)
interest. of all containers entering the appropriate segment.
One way to tackle this problem is to produce an index
Each segment’s performance is measured using fluidity
to track container transit times. This index should depict
indicators, which are computed using various statistics
the reliability and the variability of transit times but in
of the transit/dwelling time distributions for each of the
such a way as to be able to allow for performance compari-
supply chain segments, such as:
son between differing time periods.
Reliability Indicator (RI) – the ratio of the 95th percentile
To simplify the discussion, assume that the ultimate goal to the 5th percentile of transit/dwelling times (a high
is to compare quarterly and/or monthly performance data, RI indicates high volatility, whereas a low RI (≈ 1)
irrespective of the transit season, in order to determine how indicates a reliable corridor);
well the network is performing on the Shanghai → Port Buffer Index (BI) – the ratio of the positive difference be-
Metro Vancouver/Prince Rupert → Toronto corridor, say. tween the 95th percentile and the mean, to the mean.
A small BI (≈ 0) indicates only slight variability in
the upper (longer) transit/dwelling times; a large
BI indicates that the variability of the longer tran-
sit/dwelling times is high, and that outliers might be
found in that domain;
Coefficient of Variation (CV) – the ratio of the standard
deviation of transit/dwelling times to the mean tran-
sit/dwelling time.

The time series of monthly indicators (which are derived

from the monthly transit/dwelling time distributions in
each segment) are then decomposed into their

trend;
seasonal component (seasonality, trading-day, moving-
Figure 5. Multi-modal supply chain. holiday), and

[Link], [Link], [Link], [Link], [Link] (2020) Page 5 of 25

DATA SCIENCE REPORT SERIES ANOMALY DETECTION AND OUTLIER ANALYSIS

The multiplicative and pseudo-additive models are defined

in similar ways (consult [27–31] for details).8

The data decomposition/preparation process is illustrated

with the 40-month time series of marine transit CVs from
2010-2013, whose values are shown in Figure 7. The size of
the peaks and troughs seems fairly constant with respect to
the changing trend; the SAS implementation of X12 agrees
with that assessment and suggests the additive decomposi-
tion model, with no need for further data transformations.

The diagnostic plots are shown in Figure 8: the CV se-

Figure 6. Conceptual time series decomposition; potential
ries is prior-adjusted from the beginning until OCT2010
anomalous behaviour should be searched for in the irregular
component. after the detection of a level shift. The SI (Seasonal Irreg-
ular) chart shows that there are more than one irregular
component which exhibits volatility. The adjusted series
irregular component. is shown below in Figure 9; the trend and irregular com-
ponents are also shown separately for readability. It is on
The trend and the seasonal components provide the ex-
the irregular component that detection anomaly would be
pected behaviour of the indicator time series;6 the irreg-
conducted.
ular component arise as a consequence of supply chain
volatility. A high irregular component at a given time
point indicates a poor performance against expectations for
that month, which is to say, an anomalous observation. The last example shows the importance of domain under-
standing and data preparation to the anomaly detection
In general, the decomposition follows a model which is process. Given that the vast majority of observations in a
multiplicative; general problem are typically "normal", another conceptu-
additive, or ally important approach is to view anomaly detection as
pseudo-additive. a rare occurrence learning classification problem or as a
novelty detection data stream problem (these problems
The choice of a model is driven by data behaviour and will be tackled in other data science reports of this series).
choice of assumptions; the X12 model automates some of Either way, while there a number of strategies that use
the aspects of the decomposition, but manual intervention regular classification/clustering algorithms for anomaly de-
and diagnostics are still required.7 The additive model, for tection, they are rarely successful unless they are adapted
instance, assumes that: or modified for the anomaly detection context.
1. the seasonal component S t and the irregular compo-
nent I t are independent of the trend Tt ; Basic Concepts
2. the seasonal component S t remains stable from year A generic system (such as the monthly transit times from
to year; and the previous subsection, say) may be realized in normal
P12 states or in abnormal states. Normality, perhaps counter-
3. there is no seasonal fluctuation: j=1 S t+ j = 0.
intuitively, is not confined to finding the most likely state,
Mathematically, the model is expressed as: however, as infrequently occurring states could still be nor-
mal or plausible under some interpretation of the system.
Ot = Tt + S t + I t
As the authors of [12] see it, a system’s states are the re-
All components share the same dimensions and units. After sults of processes or behaviours that follow certain natural
seasonality adjustment,the seasonality adjusted series is: rules and broad principles; the observations are a manifes-
tation of these states. Data, in general, allows for inferences
SA t = Ot − S t = Tt + I t to be made about the underlying processes, which can then
6
Before carrying out seasonal adjustment, it is important to identify be tested or invalidated by the collection of additional data.
and pre-adjust for structural breaks (using the Chow test, for instance), as When the inputs are perturbed, the corresponding outputs
their presence can give rise to severe distortions in the estimation of the
Trend and Seasonal effects. Seasonal breaks occur when the usual seasonal 8
The simplest way to determine whether to use multiplicative or ad-
activity level of a particular time reporting unit changes in subsequent ditive decomposition is by graphing the time series. If the size of the
years. Trend breaks occurs when the trend in a data series is lowered or seasonal variation increases/decreases over time, multiplicative decompo-
raised for a prolonged period, either temporarily or permanently. Sources sition should be used. On the other hand, if the seasonal variation seems to
of these breaks may come from changes in government policies, strike be constant over time, additive model should be used. A pseudo-additive
actions, exceptional events, inclement weather, etc. model should be used when the data exhibits the characteristics of the
7
X12 is implemented in SAS and R, among other platforms. multiplicative series, but parameter values are close to zero.

[Link], [Link], [Link], [Link], [Link] (2020) Page 6 of 25

DATA SCIENCE REPORT SERIES ANOMALY DETECTION AND OUTLIER ANALYSIS

Figure 7. Monthly marine transit CVs and estimation summary.

Figure 8. Diagnostic plots. Note that the analysis of a time series starts with estimation of the effects of festivals and trading days.
These pre-calculated estimates are then used for prior adjustment of the series. The prior adjusted original series is subsequently
analyzed using the seasonal adjustment.

are likely to be perturbed as well; if anomalies arise from

perturbed processes, being able to identify when the pro-
cess is abnormal, that is to say, being able to capture the
various normal and abnormal processes, may lead to useful
anomaly detection.

Any supervised anomaly detection algorithm requires a

training set of historical labeled data (which may be costly
to obtain) on which to build the prediction model, and a As discussed previously, the rare occurrence problem makes
testing set on which to evaluate the model’s performance optimizing for maximum accuracy
in terms of True Positives (TP, detected anomalies that
actually arise from process abnormalities); True Negatives TN + TP
a=
(TN, predicted normal observations that indeed arise from TN + TP + FN + FP
normal processes); False Positives (FP, detected anomalies
a losing strategy; instead, algorithms attempt to minimize
corresponding to regular processes), and False Negatives
the FP rate and the FN rate under the assumption that the
(FN, predicted normal observations that are in fact the
cost of making a false negative error could be substantially
product of an abnormal process).
higher than the cost of making a false positive error.

[Link], [Link], [Link], [Link], [Link] (2020) Page 7 of 25

DATA SCIENCE REPORT SERIES ANOMALY DETECTION AND OUTLIER ANALYSIS

Figure 9. Adjusted time series components plots.

Assume that for a testing set with d = FN + TP true outliers, Another supervised approach is to estimate the relative
an anomaly detection algorithm identifies m = FP + TP abnormality of various observations: it is usually quite
suspicious observations, of which n = TP are known to difficult to estimate the probability that an observation x1
be true outliers. Performance evaluation in this context is is anomalous with any certainty, but it might be possible
often measured using: to determine that it is more likely to be anomalous than
Precision – the proportion of true outliers among the sus- another observation x2 , say (denoted by x1  x2 ).
picious observations This paradigm allows the suspicious observations to be
ranked; let ki ∈ {1, . . . , m} be the rank of the i th true outlier,
n TP i ∈ {1, . . . , n}, in the sorted list of suspicious observations
p= = ;
m FP + TP
x1  xk1  · · ·  xki  · · · xkn  xm ;
when most of the points identified by the algorithm
are true outliers, p ≈ 1; the rank power of the algorithm is
Recall – the proportion of true outliers detected by the n(n + 1)
RP = Pn .
algorithm 2 i=1 ki
n TP
r= = ; When the d actual anomalies are ranked in (or near) the
d FN + TP
when most of the true outliers are identified by the top d suspicious observations, RP ≈ 1.
algorithm, r ≈ 1; Rank power is well-defined only when m ≥ d; as with
most performance evaluation metrics, a single raw number
F1 −Score – the harmonic mean of the algorithm’s precision is meaningless – it is in comparison with the performance
and its recall of other algorithms that it is most useful.
2pr 2TP
F1 = = ;
p+r 2TP + FP + FN On the unsupervised front, where anomalous/normal la-
bels are not known or used, if anomalies are those obser-
one drawback of precision, recall, and the F1 −score is
vations that are dissimilar to other observations, and if
that they do not incorporate TN in the evaluation pro-
clusters represent groupings of similar observations, then
cess, but this is unlikely to be problematic as regular
observations that do not naturally fit into a cluster could
observations that are correctly seen as unsuspicious
be potential anomalies (see Figure 10).
are not usually the observations of interest.9 ;
There are a number of challenges, not the least of which
Example 5. Consider a test dataset with 5000 observations, being that most clustering algorithms do not recognize po-
100 of which are anomalous. An algorithm which predicts tential outliers (DBSCAN is a happy exception) and that
all observations to be anomalous would score a = p = 0.02, some appropriate measure of similarity/dissimilarity of ob-
r = 1, and F1 ≈ 0.04, whereas an algorithm that detects 10 servations has to be agreed upon (different measures could
of the true outliers would score r = 0.1 (the other values lead to different cluster assignments).
would change according to the TN and FN counts).
9 Finally, it is worth mentioning that the definitions of terms
Nevertheless, the analyst for whom the full picture is important might
want to further evaluate the algorithm with the help of the Matthews like normal and anomalous are kept purposely vague, to
TN
Correlation Coefficient [32] or the specificity s = FP+TN . allow for flexibility.

[Link], [Link], [Link], [Link], [Link] (2020) Page 8 of 25

DATA SCIENCE REPORT SERIES ANOMALY DETECTION AND OUTLIER ANALYSIS

Figure 10. Clusters of customers (red, green, blue) and potential anomalies/outliers (grey) in an artificial dataset.

1.3 Suggested References 1.4 Structure and Organization

The main references that were consulted in the preparation In this report, we aim to provide a better understanding
fo this report are: of anomaly detection, outlier techniques and some of the
field’s challenges. We also provide some step-by step appli-
Aggarwal, C.C. [2017], Outlier Analysis (2nd ed.), cations of the techniques over real-life examples.
Springer [1]
Mehrotra, K.G., Mohan, C.K., Huang, H. [2017], Ano- The purpose of Section 2 is to provide a comprehensive
maly Detection Principles and Algorithms, Springer and structured overview of different methods of anomaly
[12] detection and outliers analysis in the quantitative case;
Other good survey documents include Chandola, Banerjee, the qualitative case is tackled in Section 3. In these sec-
and Kumar’s Outlier detection: a survey [21], and Hodge tions, particular attention is paid to supervised and unsu-
and Austin’s A survey of outlier detection methodologies [22]. pervised methods, including distance-based and density-
based methods.
Specific methods and approaches are the focus of other Section 4 is dedicated to approaches for large data
papers: [2,8,41] (high-dimensional data), [7] (DOBIN), [9] sets, known as HDLSS (high dimension low sample size)
(outlier ensembles), [17,24] (isolation forest), [18,25] (DB- where the sample size n is smaller than the dimension p.
SCAN), [39] (LOF), [37, 38, 40, 42, 43] (subspace method), Outlier detection in HDLSS dataset is even more challen-
[35] (time series data). gin, mostly due to the curse of dimensionality. Feature
bagging, ensemble methods, and various dimension re-
On the practical side of things, we would be remiss if we duction methods are also discussed.
didn’t mention Arora’s An Awesome Tutorial to Learn Outlier Finally, in Sections 5 and 6, we provide detailed prac-
Detection in Python using PyOD Library [13]; note that there tical, real-life examples of anomaly detection in stock ex-
is a plethora of quality tutorials for anomaly detection in change data, and suggest project work for airline data and
the programming language of your choice online. fatal driving collision data.

[Link], [Link], [Link], [Link], [Link] (2020) Page 9 of 25

DATA SCIENCE REPORT SERIES ANOMALY DETECTION AND OUTLIER ANALYSIS

2. Quantitative Methods of Anomaly Detection If f X = f Y (or f X = f Y almost everywhere in the continuous

case, that is, except over a countable set), then
Quantitative methods are divided into distance- and density- Z
based methods. XÆ Æ
f x f Y = 1 or f X f Y dz = 1
Ω Ω
2.1 Distance-Based Methods
In order to determine whether an observation is anomalous and H(X , Y ) = 0. The fact that H(X , Y ) ∈ [0,
p1] is a
or not, it must be compared to a set of other observations consequence of Cauchy’s inequality, with f X∗ = f X and
f Y∗ = f Y :
p
(anomalies are relative, not absolute). In the distance-
based context, one natural way to compare observations is Z Z
to consider their distance from one another, with increasing Æ
0≤ f X f Y dz = f X∗ f Y∗ dz
distance from the others being increasingly suggestive of
Ω Ω
anomalous status. Z 1/2 Z 1/2
This approach works both in continuous and discrete ≤ | f X∗ |2 dz | f Y∗ |2 dz
cases, as long as a distance function or a pre-computed ta- Ω Ω
ble of pair-wise distances between observations is given. Z 1/2 Z 1/2
The choice of which sets of points to use in this compar- = f X dz f Y dz = 1;
ison distinguishes the different distance-based algorithms. Ω Ω

(a similar argument holds for discrete random variables).

This discussion begins with the introduction of some nota-
tion. Let D ⊂ Rn be an n-dimensional data set, p, q ∈ D, Recall that the covariance matrices ΣX and ΣY are n × n-
P ⊂ D be a subset of D, and d : D×D → R gives the distance matrices whose (i, j)-th entries are the covariance between
between p and q, written d(p, q). the i-th and j-th positions of X and Y , respectively. Given
An anomaly detection algorithm provides a function a collection of identically distributed samples, these covari-
a : D → R that describes how anomalous a given point is. ance matrices can be estimated.
This induces an ordering on the points of D: if a(p) < a(q) We can also consider a single point p to represent prob-
for p, q ∈ D, then p is less anomalous than q. ability distribution. In that case, the Hellinger distance
It could be necessary to define a threshold beyond which between that point and any other distribution with mean
a point is considered anomalous; if α ∈ R is such a thresh- µ and covariance matrix Σ can be studied using the frame-
old, then any p ∈ D is absolutely anomalous if a(p) > α. work above, using the Mahalanobis distance:
Æ
Similarity Measures M (p) = (p − µ)> Σ−1 (p − µ).
A similarity measure is a real-valued function that de-
scribes the similarity between two objects. A common con- Alternatively, if p and q are drawn from the same distribu-
struction is to define the similarity w between two points tion with covariance Σ, then the Mahalanobis distance is a
p, q as dissimilarity measure between p and q:
Æ
1 d M (p, q) = (p − q)> Σ−1 (p − q).
w(p, q) = , for some distance d,
1 + d(p, q)
Now, if Σ is diagonal, then
so that w → 1 as d → 0, and w → 0 as d → ∞. v
n
(pi − qi )2
uX
A similarity measure can also be constructed between d M (p, q) =
t
,
probability distributions. Let X and Y be two n-dimensional i=1
σ2i
random vectors of (possibly) different distribution with
probability mass/density functions (p.m.f./p.d.f.) f X and where σ2i is the variance along the i-th dimension. If Σ is the
f Y , respectively. Let Ω be their shared domain. For discrete identity matrix, then we recover the Euclidean distance
random variables, the Hellinger distance is defined by v
uX n
d2 (p, q) = (pi − qi )2 .
t
‚ Œ1/2
XÆ
i=1
H(X , Y ) = 1 − f X (z) f Y (z) ;
z∈Ω When using the Euclidean distance in an anomaly detection
context, a linear normalization is usually applied to each
for continuous random variables, it is defined by dimension so that each entry lies in the hypercube [−1, 1]n .
 Z 1/2
Æ
H(X , Y ) = 1 − f X (z) f Y (z) dz .
Ω

[Link], [Link], [Link], [Link], [Link] (2020) Page 10 of 25

DATA SCIENCE REPORT SERIES ANOMALY DETECTION AND OUTLIER ANALYSIS

The Minkowski distance of order p is a generalization of Distance-Based Approaches

the Euclidean distance: All these distance functions can be used to create basic
anomaly detection algorithms (the ideas can also be ex-
‚ n
X
Œ1/p tended to more complex algorithms).
d p (p, q) = |pi − qi | p

i=1 Given some distance function d, dataset D, and integers

k, ν ≤ |D|, the distance to all points anomaly detection
For p = 2 we recover the Euclidean distance d2 , for p = 1
algorithm considers each point p in D and adds the distance
the Manhattan distance
from p to every other point in D, i.e.
n X
a(p) =
X
d1 (p, q) = |pi − qi |, d(q, p).
q6=p∈D
i=1
The ν points with largest values for a are then said to be
and for p = ∞ the supremum distance anomalous according to a. This approach often selects
n the most extreme observations as anomalous, which may
d∞ (p, q) = max |pi − qi |. be of limited use in practice.
i=1
The distance to nearest neighbour algorithm defines
The Minkowski distance d p is only actually a distance func-
a(p) = min d(q, p),
tion (i.e., a metric) when p ≥ 1, but an exception is made q6=p∈D
for with a similar definition for the ν anomalous points.
n
d−∞ (p, q) = min |pi − qi | The average distance to k nearest neighbours and
i=1
median distance to k nearest neighbours are defined sim-
to fall within the same framework. ilarly.

The Jaccard similarity of two datasets P and Q, is defined 2.2 Density-Based Methods
as the size of their intersection divided by the size of their Density-based approaches, on the other hand, view points
union as anomalous if they occur in low density regions.

|P ∩ Q| |P ∩ Q| Local Outlier Factor

J(P, Q) = = The Local Outlier Factor (LOF) algorithm was proposed
|P ∪ Q| |P| + |Q| − |P ∩ Q|
in 2000 by [39] (a summary can be found in Section 6.4.2
Their Jaccard distance is then taken to be 1 − J(P, Q). of [12]). LOF works by measuring the local deviation of
This definition can be extended to compare binary vec- each point in a dataset from its k nearest neighbours, with
tors (i.e. vectors with entries in {0, 1}) of the same length. a point said to be anomalous if this deviation is large.
Given two binary vectors p and q of length n, consider an A local k−region around a point p is defined as the k
arbitrary set D of size n. Then p and q can be viewed as nearest neighbours of p. The density of points in each of
subsets of D: if pi = 1 then p is said to contain the i-th their respective local k−neighbourhoods is estimated, and
element of D, while if pi = 0 then it does not. Viewing compared to the density of the local k−neighbourhoods of
p and q in this way allows us to compute their Jaccard the point within their own k−neighbourhood.
similarity, and thus their Jaccard distance. This can then be used to identify outliers that inhabit re-
gions of lower density than their neighbours, as p would be
Finally, let p, q 6= 0. Recall that p · q = kpkkqk cos θ , where in Figure 11. The formal procedure is shown in Algorithm 1.
θ is the angle between p and q. The cosine similarity be-
tween p and q is the cosine of θ , which can be computed
as Pn
p·q i=1 pi qi
cos θ = = qP .
kpkkqk n
qP n
p2 q2
i=1 i i=1 i
p
This value ranges between 1 and −1, with 1 attained when
p = q, −1 when p = −q, and 0 when p and q are perpen-
dicular.

Armed with these concepts, we can now explore distance-

(and eventually density-) based methods for anomaly de- Figure 11. For k = 2, p is an outlier as it has lower density
tection. than its neighbours.

[Link], [Link], [Link], [Link], [Link] (2020) Page 11 of 25

DATA SCIENCE REPORT SERIES ANOMALY DETECTION AND OUTLIER ANALYSIS

Algorithm 1: Local Outlier Factor (LOF) DBSCAN

Density-Based Spatial Clustering of Applications with Noise
1 Input: dataset D, point p ∈ D, integer k for
(DBSCAN) was proposed in 1996 by [18] (a summary can
number of nearest neighbours to consider,
be found in Section 4.1.5 of [12]). As its name suggests, it
distance function d
is a density-based clustering algorithm that groups nearby
2 Compute the distance between all points in D
points together and labels points that do not fall in the
3 for p ∈ D do
clusters as anomalies.
4 for q ∈ D \ {p} do
Hierarchical DBSCAN (HDBSCAN) [25] was introduced
5 Compute d(p, q)
in 2013. It notably removes the problem of choosing the pa-
6 end
rameter for the radius of a neighbourhood by considering all
7 Order D by increasing distance from p
possible radii. Further documentation can be found at [26].
8 Set dk (p) = d(p, qk )
9 end
In DBSCAN,
10 Find the k nearest neighbours of p
11 Set Nk (p) = {q ∈ D \ {p} : d(p, q) ≤ dk (p)} a point p is a core point if there are a minimum
12 Define the reachability distance number m of points within distance r of p;
dreach (p, q) = max{dk (q), d(p, q)} a point q is a border point if it is not itself a core
13 Define the average reachability distance point but is within distance r of one, and
dreach (p,q)
P
q∈N (p)
dreach (p) = k
|Nk (p)|
a point o is an outlier if it is neither a core nor a
14 Define the local reachability density border point.
€ Š−1
`k (p) = dreach (p) DBSCAN considers each point in the dataset individually. If
`k (q)
P
q∈Nk (p) `k (p)
that point is an outlier, then it is added to a list of outliers.
15 Compute the local outlier factor ak (p) = |Nk (p)| Otherwise if it is a core point, then its r-neighbourhood
16 Output: LOF ak (p) forms the beginning of a new cluster. Each point in this
r-neighbourhood is then considered in turn, with the r-
neighbourhoods of other core points contained in the neigh-
LOF is able to identify local outliers, but selecting a thresh- bourhood being added to the cluster.
old beyond which a point is considered an outlier is difficult. This expansion repeats until all points have been exam-
ined. During this step points that were previously labelled
LOF introduces the idea of a reachability distance, which as outliers may be updated as they become border points
improves the stability of results within clusters/regions: in this new cluster. This process continues until every point
within a local k−region around p, it is simply the maximal has either been assigned to a cluster or labelled as an outlier
distance to its k−neighbours; outside of that region, it is (see Algorithm 2 and Figure 13).
the actual distance from p. While DBSCAN’s dual use as a clustering algorithm may
In Figure 12 (with k = 3), for instance, the points seem irrelevant in the outlier detection setting, its ability to
q1 , q2 , q3 all have the same reachability distance from p succesfully identify clusters is crucial to being able to label
as they are all 3-neighbours of p, that is, the remaining points as outliers.
dreach (p, q1 ) = dreach (p, q2 ) = dreach (p, q3 ) = d(p, q3 ).
On the one hand, in DBSCAN the number of clusters does
The point q4 , on the other hand, has dreach (p, q4 ) = d(p, q4 ) not need to be known beforehand (unlike in k−means and
as it is not a k-neighbour of p. other clustering algorithms) and clusters of arbitrary shape
can be detected.
Furthermore, when using HDBSCAN, only the parame-
q4 ter for the minimum cluster size m is required, which can
q3

p q1

o q1 p q2
q2

Figure 13. For minimum neighbourhood size m = 2 and

Figure 12. The region of uniform reachability distance this fixed radius r, o is an outlier, p a core point, and q1
around p for k = 3. and q2 are border points.

[Link], [Link], [Link], [Link], [Link] (2020) Page 12 of 25

DATA SCIENCE REPORT SERIES ANOMALY DETECTION AND OUTLIER ANALYSIS

Algorithm 2: DBSCAN cluster. This issue is overcome in HDBSCAN.

1 Input: dataset D, distance function d, Isolation Forest
neighbourhood radius r > 0, minimum number of The previously discussed approaches first construct models
points to be considered a cluster m ∈ N of what normal points look like, and then identify points
2 Clusters = {} that do not fit this model. The Isolation Forest algorithm
3 Outliers = {} [17] introduced in 2008 instead tries to explicitly identify
4 for p ∈ D do outliers under the assumptions that there are few outliers
5 if p ∈ Outliers ∪ (∪C∈Clusters C) then and that these outliers have very different attributes com-
6 continue pared to normal points. Doing so allows the use of sampling
7 end techniques that increase algorithmic speed while decreasing
8 Set N (p) = {q ∈ D : d(p, q) ≤ r} memory requirements.
9 if |N (p)| < m then The Isolation Forest algorithm tries to isolate anoma-
10 Add p to Outliers lous points. It does this by randomly selecting an attribute
11 continue and then randomly selecting a split value between that at-
12 end tribute’s min and max values. This recursively partitions
13 else the points until every point is isolated in its own partition.
14 Cluster = N (p)
15 for q ∈ Cluster \ {p} do Recursive partitioning yields a binary tree called an Iso-
16 if q ∈ Outliers then lation Tree. The root of this tree is the entire dataset; each
17 Remove q from Outliers node is a subset of the observations, and each branch cor-
18 end responds to one of the generated partitions. The leaf nodes
19 else if q ∈ ∪C∈Clusters C then are singleton sets containing a single isolated point. Each
20 continue point is then assigned a score derived from how deep in
21 end the tree its singleton partition appears (see Figure 14 and
22 Set N (q) = {q0 ∈ D : d(q, q0 ) ≤ r} Algorithm 3).
23 if |N (q)| ≥ m then As points that are shallower in the tree were easier to
24 Cluster = Cluster ∪ N (q) separate from the rest, these are the likely outliers. Since
25 end only shallow points are of interest, once the height of the
26 end tree has reached a given threshold (the expected height of
27 end a random binary tree, say), further construction of the tree
28 Add Cluster to Clusters can be stopped to decrease computational cost.
29 end Additionally, instead of building a tree from the entire
30 return Outliers dataset, a tree can be constructed from a subset. The loca-
31 Output: a list of outliers tion of any point within this smaller tree can then be esti-
mated, again saving computational and memory resources.
These two improvements are detailed in the original pa-
per [17].
be set fairly intuitively, which is not the case for the param-
eters in general clustering algorithms: if the elements of Once a number of Isolation Trees have been randomly gen-
D are n−dimensional, take m ≥ n + 1 (larger values of m erated (Isolation Forest), a score can be computed for each
allow for better noise identification). point. This is done by searching each tree for the location
of a given point and noting the path length required to
On the other hand, DBSCAN is not deterministic, as border reach it. Once a point’s path length in each tree has been
points can be assigned to different clusters depending on computed, the average path length is taken to be its score.
the order in which core points are considered (this does not
affect its use as an anomaly detection algorithm, however).
In high dimensions, the ability of any distance function
based on Euclidean distance to distinguish near and distant
points diminishes due to the Curse of Dimensionality; thus
in high dimension spaces, it become ineffective (as do other
clustering algorithms).
Finally, DBSCAN cannot handle differences in local den-
sities as the radius of a neighbourhood r is fixed; this could
lead to sparser clusters being labelled as outliers, or to out- Figure 14. A partitioning constructed during Isolation
liers surrounding a denser cluster being included in the Tree generation.

[Link], [Link], [Link], [Link], [Link] (2020) Page 13 of 25

DATA SCIENCE REPORT SERIES ANOMALY DETECTION AND OUTLIER ANALYSIS

Algorithm 3: Recursive Isolation Tree Construc- Algorithm 4: Isolation Forest

tion: iTree(D) 1 Input: dataset D, integer t number of Isolation
1 Input: dataset D Trees
2 if |D| ≤ 1 then 2 Forest = {}
3 return {} 3 for i = 1 to t do
4 end 4 Tree = iTree(D)
5 else 5 Add Tree to Forest
6 Let A be a list of attributes in D 6 end
7 Randomly select an attribute A ∈ A 7 for p ∈ D do
8 Randomly sample a point s from 8 PathLengths = {}
[minq∈D A(q), maxq∈D A(q)] 9 for Tree in Forest do
9 Return 10 Find the path length ` from the root of Tree
to node {p}
= iTree({q ∈ D : A(q) ≤ s})

LeftChild 11 Add ` to PathLengths
Node RightChild = iTree({q ∈ D : A(q) > s})
12 end
NodeValue = D

`∈PathLengths `
P
13 AveragePathLength = t
10 end AveragePathLength
11 Output: Binary tree with node values that are 14 Set a(p) = 2− c(|D|)

subsets of D 15 end
16 Output: Anomaly score a(p) ∈ [0, 1] for each
p∈D
It can be desirable to construct a normalized anomaly score
that is independent of the size of the dataset. In order to
do this, the expected path length of a random point in an 3. Qualitative Methods of Anomaly Detection
Isolation Tree (i.e. binary tree) must be estimated. With
n = |D|, it can be shown that the expected length is New challenges are presented by non-numerical variables.

2(n − 1)
c(n) = 2H(n − 1) − , 3.1 Definitions and Challenges
n Categorical Variables
where H(n − 1) is the (n − 1)th harmonic number, which A categorical variable (or qualitative variable) is one whose
can be approximated by ln(n − 1) + 0.577; c(n) is then used levels are measured on a nominal scale; examples include
to normalize the final anomaly score a(p) for p ∈ D, which an object’s colour, the mother tongue of an individual, her
is given by favourite meal, and so forth.
The central tendency of the values of a categorical
average path length to p in the Isolation Trees variable is usually given by its mode; measures of spread
log2 a(p) = − .
c(n) are harder to define consistently (the proportion of levels
with more than a certain percentage of the observations
Thus defined, a(p) ∈ [0, 1], with a(p) ≈ 1 suggesting p is above a given threshold could be used as rough gauge, but
an anomaly, a(p) ≤ 0.5 suggesting p is a normal point; difficulties with this approach are readily apparent).
if all points receive a score around 0.5, this suggests that We often associate qualitative feature to numerical val-
there are no anomalies present. ues, but with the caveat that these should not be interpreted
as numerals; if we use the code “red” = 1 and “blond” = 2
Isolation Forests have small time and memory requirements; to represent hair colour, for instance, we obviously cannot
can handle high dimensional data, and do not need observa- conclude that “blond” > “red”, even though 2 > 1.
tions to have been labeled anomalies in the training set, but A categorical variable that has exactly two levels is
the anomaly score assigned to a given point can have high called a dichotomous feature (or a binairy variable); those
variance over multiple runs of the algorithm. The authors with more than two levels are called polytomous vari-
of [24] propose some solutions. ables.

Challenges of Anomaly Detection with Categorical Data

In general, density-based schemes are more powerful than Representing categorical variables with numerical features
distance-based schemes when a dataset contains patterns can lead to traps; consequently, using anomaly detection
with diverse characteristics, but less effective when the methods based on distance metrics or on density is not
patterns are of comparable densities with the outliers [44]. recommended in the qualitiative context, unless they have
first been modified appropriately.

[Link], [Link], [Link], [Link], [Link] (2020) Page 14 of 25

DATA SCIENCE REPORT SERIES ANOMALY DETECTION AND OUTLIER ANALYSIS

3.2 Review of Two Methods Algorithm 5: AVF

We present two of the categorical methods below.
1 Inputs: dataset D (n observations, m features),
AVF Algorithm number of anomalous observations k
The Attribute Value Frequency (AVF) algorithm offers a 2 while i ≤ n do
fast and simple way to detect outlying observations in cate- 3 j=1
gorical data, which minimizes the amount of data analyses, 4 AVFscore(xi ) = f (x i, j )
without having to create or search through various combi- 5 while j ≤ m do
nations feature levels (which increase the search time). 6 AVFscore(xi ) = AVFscore(xi ) + f (x i, j );
j = j+1
Intuitively, outlying observations are points which occur rel- 7 end
atively infrequently in the (categorical) dataset; an “ideal” 8 AVFscore(xi ) =Mean(AVFscore(xi ))
anomalous point is one for which each feature value is 9 i = i+1
extremely anomalous (or relatively infrequent). 10 end
11 Outputs: k observations with smallest AVF scores
The rarity of an attribute level can be measured by sum-
ming the number of times the corresponding feature takes
that value in the dataset. where the entropy of the subset D \ O is the sum of the
Let’s say that there are n observations in the dataset: entropy of each of feature on D \ O:
{xi }, i = 1, . . . , n, and that each observation is a collection
of m features. We write H(D \ O) = H(X 1 ; D \ O) + · · · + H(X m ; D \ O)
xi = (x i,1 , · · · , x i,` , · · · , x i,m ), et X
where x i,` `th feature’s level. Using the reasoning presented H(X ` ; D \ O) = − p(z` ) log p(z` ),
above, the AVF score (shown below) is a good tool to z` ∈S(X ` ;D\O)

determine whether xi should be considered an outlier or where S(X ` ; D \ O) is the set of levels that the `th feature
not: takes in D \ O.
m
1 X
AVFscore(xi ) = f (x i,` ),
m `=1 The "greedyAlg1" algorithm solves the optimization prob-
where f (x i,` ) is the number of observations xi for which lem as follows:
the `th feature takes on the level x i,` . A low AVF score
1. The set of outlying and/or anomalous observations
indicates that the observation is more likely to be an outlier.
OS is initially set to be empty, and all observations of
D \ OS are identified as normal (or regular).
Since AVFscore(xi ) is essentially a sum of m positive num-
bers, it is minimized when each of the sum’s term is mini- 2. Compute H(D \ OS).
mized, individually. Thus, the “ideal” anomalous observa-
tion described above minimizes the AVF score; the minimal 3. Scan the dataset in order to select a candidate anoma-
score is reached when each of the observation’s features’ lous observation: every normal observation x is tem-
levels occurs only once in the dataset. porarily taken out of D \ OS to create a subset Dx0 ,
whose entropy H(Dx0 ) is also computed.
As shown by the AVF pseudocode (see Algorithm 5), once
the AVF score is calculated for all points, the k outliers 4. The observation z which provides the maximal en-
returned by the algorithms are the k observations with the tropy impact, i.e. the one that minimizes
smallest AVF scores (the algorithm’s complexity is O (nm)).
H(D \ OS) − H(Dx0 ), x ∈ D \ OS,
Greedy Algorithm
The greedy algorithm “greedyAlg1” is an algorithm which is added to OS.
identifies the set OS of candidate anomalous observations
in an efficient manner. 5. Repeat steps 2-4 another k − 1 times to obtain a set
OS of k candidate anomalous observations.
The mathematical formulation of the problem is simple
– given a dataset D and a number k of anomalous observa- You can find more details in the in the source article [46];
tions to identify, we solve the optimization problem an interesting detail is that the complexity of the algorithm
is expected to be O (nmp), which implies that it is scalable.
OS = arg min{H(D \ O)}, subject to |O| = k,
O⊆D

[Link], [Link], [Link], [Link], [Link] (2020) Page 15 of 25

DATA SCIENCE REPORT SERIES ANOMALY DETECTION AND OUTLIER ANALYSIS

4. Anomalies in High-Dimensional Datasets 4.2 Projection-Based Methods

Nowadays, it is common to deal with very large data sets
Anomaly detection is a broad field of study that has been known as HDLSS (high dimension, low sample size),
applied to a large number of areas. Nowadays, many real which can contain hundreds of variables (or even more).
datasets are very large; in some scenarios, the observa- As a result, the curse of dimensionality affects the effi-
tions may contain hundreds or thousands of features (or ciency of conventional anomaly/outlier detection methods.
dimensions).
Many classical methods use proximity (distance) con- One solution to this problem is to reduce the dimension-
cepts for anomaly detection (see Section 2 for a sample of ality of the dataset while preserving its essential character-
such methods) and can only be applied in cases where the istics. Such projecion-based methods
sample size n is larger than the dimension p (n > p).
principal component analysis,
The management of high-dimensional data (n < p) of- linear discriminant analysis,
fers specific difficulties: indeed, in such spaces observations feature selection, etc.
are often isolated and scattered (or sparse) and the notion In this section, We provide details on one such method:
of proximity fails to maintain its relevance. PCA.
In that case, the notion of defining significant outliers is
much more complex and not obvious: many conventional Principal Components Analysis
methods of detecting outliers are simply not efficient in the Principal components analysis (PCA) aims to find a rep-
high-dimensional context, due to this curse of dimension- resentation of the original dataset in a lower-dimensional
ality. subspace (such as a line or a plane) containing the greatest
possible variation.
The remainder of this section is organized as follows: first, PCA corresponds to an orthogonal linear transforma-
an attempt is made to define the concept and the challenges; tion of the data into a new coordinate system, such that the
then, anomaly detection techniques are discussed; finally, largest variance resulting from a scalar projection of the
we end with a detailed description of ensembles and sub- data is on the first coordinate (the first principal compo-
space methods. Our approach mainly follows those found nent), the second largest variance on the second coordinate,
in [1, 9, 10, 12, 16]. and so forth.

4.1 Definitions and Challenges PCA is used in various contexts:

As we have seen previously, an anomalous observation is
as a dimension reduction method used during the
one that deviates or behaves differently from other the
data pre-processing step;
observations in the dataset, which makes us suspect that
as a data visualization aid, and, in the scenario of
it was generated by some other mechanism [1]; such an
interest for this report,
observation would, of course, be considered to be irregular.
as an anomaly and outlier detection method.
The challenges of anomaly and outlier detection in high-
dimensional data lie in the facts that: Let the dataset be represented by a numerical, centered,
and scaled n × p matrix X = [X1 , · · · , X p ] with n observa-
the notion of distance fails to retain its relevance due tions (number of rows) and p features (number of columns).
to the curse of dimensionality (whence “the prob- The principal components can be written as linear combi-
lem of detecting outliers is like finding a needle in a nations of the variables
haystack” [16]);
every point in such datasets has a tendency to be an Yi = `>i X = `1,i X1 + · · · + ` p,i X p ; i = 1, · · · , k,
outlier, and
datasets become more sparse as the dimension of the with k ≤ p, yielding the largest variance subjet to the con-
feature space increases. straint k`i k = 1 (where k·k represents the Euclidean norm).
We can thus deduce that
The authors of [2] consider that in order to deal properly
Var (Yi ) = Var `>i X = `>i Σ`i ,


with large datasets, detection methods should:
Cov (Yi , Yk ) = Cov `>i X , `>k X = `>i Σ`k .


1. allow for effective management of sparse data issues;
2. provide interpretability of the discrepancies (i.e. how In other words, PCA finds the loadings vector `1 which
the behaviour of such observations is different); maximizes the variance of Y1 , i.e.
3. allow anomalie measurements to be compared, and
`1 = arg max `>1 X> X`1 ,

4. consider the local data behaviour to determine whether
an observation is abnormal or not. k`1 k=1

[Link], [Link], [Link], [Link], [Link] (2020) Page 16 of 25

DATA SCIENCE REPORT SERIES ANOMALY DETECTION AND OUTLIER ANALYSIS

then the loadings vector `2 (not correlated with `1 ) which Example 6. PCA is applied on a dataset of genetic ex-
maximizes the variance of Y2 , i.e. pression measurements, for n = 72 leukemia patients and
p = 7128 genes [47]. The scree plot suggests that only one
`2 = arg max `>2 X> X`2 .

k`2 k=1, `>1 `2 =0
principal component should be retained; the projection on
the first 3 principal components is also shown in Figure 15
Similarly, the loadings vector `k is not correlated with any (on the right). Some R code is given below.
of the `i , i < k, and maximizes the variance of Yk , i.e.
[Link] <-
`k = arg max `>k X> X`k .

(1) [Link]("[Link]
k`k k=1,
`>i `k =0, ∀ i<k CASI_files/DATA/leukemia_big.csv")
[Link] <- t([Link])
We solve (1) for all i < k through the Lagrangian [Link] <-
scale([Link])
L = `>k X> X`k − λk (`>k `k − 1) − w`>i `k .
[Link] <-
The critical points are found by differentiating with respect prcomp([Link])
to each of the entries of `k , λk and w, and setting the result plot([Link])
to 0. Simplifying, we obtain [Link].s <- summary([Link])
plot([Link].s$importance[3,])
X> X`k = λk `k
`>k `k = 1 et `>k `i = 0, for all i < k. There are other PCA-associated dimension reduction meth-
The loadings vector `k is thus the eigenvector of the design ods, such as the singular value decomposition, kernel PCA,
matrix X> X associated to the kth largest eigenvalue. and so forth; more details are available in [48].

The proportion of the variance which can be explained What is the link with anomaly and/or outlier detection?
by the PCA can be calculated by first noting that Once the dataset has been projected on a lower-dimensional
p p p
subspace, the curse of dimensionality is mitigated – it is on
the projected data that the traditional detection methods
X X X
Var (Yi ) = `>i Σ`i = λi .
i=1 i=1 i=1
are applied.
Note, however, that any such reduction necessarily leads
Consequently, the proportion of the total variance explained
to a loss of information, which can affect the accuracy of
by the ith principal component is
the detection procedure, especially if the presence/absence
λi of anomalies is not aligned with the dataset’s principal
0 ≤ Pp ≤1
components.
i=1 λi

The quality of the PCA results is strongly dependent on the Distance-Based Outlier Basis Using Neighbours
number of retained principal components, that is, on the Using PCA for anomaly detection is potentially problematic,
dimension k of the subspace on which the observations are however: whether an observation is anomalous or not does
projected. There are multiple ways to select the “right” k – not figure in the construction of the principal component
we will briefly present two of them. basis {PC1 , . . . , PCk } – there isn’t necessarily a correlation
between the axes of heightened variance and the presence
The proportion of the total variance explained by the first or absence of anomalies.
k principal components is given by
Pk The distance-based outlier basis using neighbours algo-
i=1 λi rithm (DOBIN) builds a basis which is better suited for the
pk = P p .
i=1 λi eventual detection of outlying observations. DOBIN’s main
One approach is to retain k principal components, where idea is to search for nearest neighbours that are in fact
k is the smallest value for which pk surpasses some pre- relatively distant from one another:
established threshold (often taken between 80% and 90%).
1. We start by building a space Y = {y` } which contains
M  n(n + 1)/2 vectors of the form
The scree plot method, on the other hand, consists in
drawing the curve given by the decreasing eigenvalues (the y` = (xi − x j ) (xi − x j ),
scree plot), and to identify the curve’s “elbows”. These
points correspond to principal components for which the where is the element-by-element Hadamard multi-
variance decreases at a slower rate with added components. plication, and for which the 1−norm
If such an elbow exists, we would retain the eigenvalues up
to it (and thus, the corresponding principal components). ky` k1 = (x 1,1 − x 2,1 )2 + · · · + (x 1,p − x 2,p )2

[Link], [Link], [Link], [Link], [Link] (2020) Page 17 of 25

DATA SCIENCE REPORT SERIES ANOMALY DETECTION AND OUTLIER ANALYSIS

Figure 15. Scree plot (left); projection on the first 3 principal components (right).

is the square of the distance between xi , x j ∈ X (the 4.3 Ensembles Methods

selection of each of the M observation pairs is made In the preceding sections, we have described various anomaly
according to a rather complex procedure which only detection algorithms whose relative performance varies
considers xi and x j if they are part of one another’s with the type of data being considered. It’s usually impossi-
k−neighbourhood, for k ∈ {k1 , . . . , k2 }); the set Y ble to come up with an algorithm that outperforms all the
thus contains points for which ky` k1 is relatively large, others.
which is to say that the observations xi are x j fairly dis- This is because a particular anomaly detection algorithm
tant from one another even if they are k−neighbours may be well adapted to a data set and may be successful
of each other; in detecting abnormal or outlier observations, but it may
2. we next build a basis {η1 , . . . , η p } ⊂ R p where each not work with other data sets whose characteristics do not
ηi is a unit vector given by a particular linear com- match the first data set.
bination of points in Y; they can be found using a
Gram-Schmidt-like procedure: The impact of such a mismatch between algorithms can be
mitigated by using ensemble methods, where the results
y`0 = y` , ` = 1, . . . , M of several algorithms are considered before making a final
y` b−1 = y` b−2 − 〈η b−1 | y` b−2 〉, ` = 1, . . . , M decision. Such an approach often provides the best results
PM and thus improves the performance of the base anomaly
`=1 y` b−1 detection algorithms [12].
ηb = P ,
M
y
`=1 ` p−1
2
We will consider two tyes of ensemble methods: sequential
for b = 1, . . . , p, ensembles (boosting) and independent ensembles,
3. and we tranform the original dataset X according to
X̂ = T (X)Θ, where T (X) normalizes each feature of Sequential Ensembles
X according to a problem-specific scheme (Min-Max Sequential ensembles requires a given algorithm (or a set of
or Median-IQR, say) and algorithms) to be applied to a dataset in a sequential man-
ntter, each time on a slightly different dataset derived from
Θ = [η1 | · · · | η p ] the previous step’s dataset based on the previous steps’ re-
sults, and so forth. At each step, the weight associated with
is a orthogonal p × p matrix.
each observation is modified according to the preceding
It is on the transformed space (which plays an analogous results using some “boosting” method (such as AdaBoost
role to the subspace projection of X in PCA) that we apply or XGBoost, for instance).
the various outlier and anomaly detection algorithms.
The final result is either some weighted combination of
The full details contain a fair number of technical complica- all preceding results, or simply the results output by the
tions; the interested reader is invited to consult the original last step in the sequence (see Algorithme 6).
documentation [7] (note that the algorithm is implemented The details are out-of-scope for this report, but can be
in R via the module dobin). studied in [49].

[Link], [Link], [Link], [Link], [Link] (2020) Page 18 of 25

DATA SCIENCE REPORT SERIES ANOMALY DETECTION AND OUTLIER ANALYSIS

Algorithm 6: SequentialEnsemble Let αi (p) represent the (normalized) anomaly score of

p ∈ D, according to algorithm Ai . If αi (p) ≈ 0, it is unlikely
1 Inputs: dataset D, base algorithms A1 , . . . , A r
that p is an anomaly according to Ai , whereas if αi (p) ≈ 1,
2 j = 1;
it is quite likely that p according to Ai .
3 while stopping criteria are not met do
The rank of p ∈ D according to Ai , on the other hand,
4 Select an algorithm A j based on the results from
is denoted by ri (p): the higher the rank (smaller number),
the preceding steps;
the higher the anomaly score vice versa. In a dataset with n
5 Create a new dataset D j from D by modifying
observations, the rank varies from 1 to n (ties are allowed).
the weight of each observation based on the
If the base detection algorithms are A1 , . . . , Am , the anomaly
results from the preceding steps;
score and the rank of an observation p ∈ D according to
6 Apply A j to D j ;
the independent ensemble method are, respectively,
7 j = j + 1;
8 end m
1 X
9 Output: anomalous observations obtained by α(p) = αi (p) and r(p) = min {ri (p)}.
m i=1 1≤i≤m
weighing the results of all previous steps
If n = m = 3, for instance, we could end up with
Algorithm 7: IndependantEnsemble α1 (p1 ) = 1.0, α1 (p2 ) = 0.9, α1 (p3 ) = 0.0;
1 Inputs: dataset D, base algorithms A1 , . . . , A r α2 (p1 ) = 1.0, α2 (p2 ) = 0.8, α2 (p3 ) = 0.0;
2 j = 1;
α3 (p1 ) = 0.1, α3 (p2 ) = 1.0, α3 (p3 ) = 0.0.
3 while stopping criteria are not met do
4 Select an algorithm A j ; Using the mean as the combination techniques, we obtain
5 Create a new dataset D j from D by (potential)
re-sampling, but independently of the α (p1 ) = 0.7, α (p2 ) = 0.9, α (p3 ) = 0.0,
preceding steps’ results;
6 Apply A j to D j ; whence
7 j = j + 1; p2  p1  p3 ,
8 end
that is, p2 is more anomalous than p1 , which is itself more
9 Output: anomalous observations obtained by
anomalous than p3 (see the notation introduced on page 8).
combining the results of all previous steps
Using the minimal rank method, we obtain

Independent Ensembles r1 (p1 ) = 1, r1 (p2 ) = 2, r1 (p3 ) = 3;

In an independent ensemble, we instead apply different al- r2 (p1 ) = 1, r2 (p2 ) = 2, r2 (p3 ) = 3;
gorithms (or different instanciations of the same algorithm)
r3 (p1 ) = 2, r3 (p2 ) = 1, r3 (p3 ) = 3,
to the dataset (or some resampled dataset).
Choices made at the data and algorithm level are inde- from which
pendent of the results obtained in previous runs (unlike in
a sequential ensemble). The results are then combined to r (p1 ) = r (p2 ) = 1, r (p3 ) = 3,
obtain more robust outliers (see Algorithm 7).
whence p1  p3 and p2  p3 , but p1 and p2 have the same
Every base anomaly detection algorithm provides an ano- anomalous levels.
maly score (or an abnormal/regular classification) for each
observation in D; observations with higher scores are con- Evidently, the results depend not only on the data set under
sidered to be more anomalous, observations with lower consideration and on the base algorithms that are used in
scores more normal. the ensemble, but also on how the results are combined.
The results are then combined using a task-specific
method in order to provide a more robust classification In the context of HDLSS data, ensemble methods can some-
of anomalous or outlying observations. times allow the analyst to mitigate some of the effects of the
curse of dimensionality by selecting fast base algorithms
Many such combination techniques used in practice: (which can be run multiple times) and focusing on building
robust relative anomaly scores.
majority vote, Another suggested approach is to use a different subset
average, of the original dataset’s features at each step, in order to
minimal rank, etc. de-correlate the base detection models.

[Link], [Link], [Link], [Link], [Link] (2020) Page 19 of 25

DATA SCIENCE REPORT SERIES ANOMALY DETECTION AND OUTLIER ANALYSIS

4.4 Subspace Methods 5. Applications to Time Series

Subspace methods have been used particularly effectively
by analysts for anomaly and outlier detection in high-dimen- In this section, we discuss outliers and anomalies in time
sional data sets [9, 15, 16]; it is often easier to find the series. A time series is a sequential set of values tracked
sought-after observations by exploring lower-dimensional over a time period. The additional structure of time series
subspaces (rather than the original set). makes detection of anomalies challenging; yet many algo-
There is thus an intrinsic interest in exploring subspaces rithms attempt the task (R provides implementations of a
in their own right [1, 10]. This approach eliminates addi- number of these approaches). In this section, we discuss
tive noise effects often found in high dimensional spaces two of the most commonly-used methods.
and leads to more robust outliers (that is, outliers which
are identified as such even when using different methods). 5.1 Outliers and Anomalies in Time Series
Outliers in time series are sudden changes in the dynamics
The problem is rather difficult to solve effectively and effi- of the data that can be temporary or permanent. These
ciently, since the potential number of subspace projections anomalous recordings are usually inconsistent with the rest
of high-dimensional data is related exponentially to the of the series and cannot be explained by standard time
number of features in the dataset. series models. If left as is, they can have a tremendous
impact on the analysis (on model selection and parameter
The Feature Bagging algorithm formalizes the idea pre- estimation, for instance).
sented at the end of the preceding sub-section; it officially They may, therefore, affect the forecasting power of the
uses the LOF algorithm of Section 2, but any fast anomaly fitted model. It is thus important to detect and treat outliers
detection algorithm can be used instead. The anomaly in time series before fitting a model.
scores and rankings from each run are aggregated as they
are in the Independent Ensemble approach. Outliers that only change the mean level of the series are
said to be deterministic. A simple procedure can be used
to detect deterministic outliers in time series: compare a
Algorithm 8: FeatureBagging
time series model with no outliers to a model that includes
1 Input: dataset D the outliers [34]. We can then estimate the effect of pro-
2 j = 1; cessing the anomalies by looking at the differences between
3 while stopping criteria are not met do the models.
4 Sample an integer r between p/2 et p − 1;
5 Randomly select r features (variables) of D in In general, detecting outliers in applied time series con-
order to create a projected dataset D̃r in the sists of determining the location, type, and magnitude of
corresponding r−dimensional sub-space; any existing outliers. There are several types of outliers:
6 Compute the LOF result for each observation in
the projected D̃r ; an additive outlier (AO) is an abrupt change for only
7 j = j + 1; one observed value – such an outlier has no effect on
8 end the subsequent observations;10
9 Output: anomaly scores given by the independent an innovational outlier (IO) is an unsual innova-
ensemble method (average, minimal rank, etc.). tion11 in the generating process that affects all later
observation – the influence of such outliers may in-
crease with the passage of time;
There are other, more sophisticated, subspace anomaly de- a level shift outlier (LS) affects the mean level of
tection methods, including: observations so that all the observations after the
outlier shift to a new level – clearly such outliers
High-dimensional Outlying Subspaces (HOS) [37]; have a permanent effect on the time series, and it is
Subspace Outlier Degree (SOD) [38]; important to detect and process them prior to building
Projected Clustering Ensembles (OutRank) [40]; any forecasting model;12
Local Selection of Subspace Projections (OUTRES) a transient change outlier (TC) is similar to a LS
[42]. but its effect is not permanent and disappears over
subsequent observations.
It should be noted that anomaly detection and outlier analy-
10
sis is still very active as an area of research, with numerous An AO is a seasonal additive outlier (SAO) when the additive outlier
challenges. The “No Free Lunch” Theorem suggests that, reappears at regular intervals.
11
The innovations in time series play the same role as errors in cross-
importantly, there is no magic method: all methods have sectional analysis (such as OLS).
strengths and limitations, and the results depend heavily 12
A LS is called a seasonal level shift outlier (SLS) when the mean
on the data. level shift occurs at regular intervals.

[Link], [Link], [Link], [Link], [Link] (2020) Page 20 of 25

DATA SCIENCE REPORT SERIES ANOMALY DETECTION AND OUTLIER ANALYSIS

Essentially, we view anomalies as outliers; regularly oc-

curing events are part of the the series trend, while rarely
occuring deviances from the trend are anomalous.

The rest of this section is devoted to application of R pack-

ages to detect outliers and anomalies.

5.2 R Package: tsoutliers

The package tsoutlier automatically detects outliers in
time series based on the procedure outlined in [35]. In this
approach, the outlier effects are estimated simultaneously
using multiple ARIMA-based regression, and the model
parameters and the outlier effects are estimated jointly.
The main interface to the automatic detection proce-
Figure 16. Outliers detection for Apple Inc. daily closing return.
dure is through the funtion tso(x, types), where x The x−axis is labelled by the date index.
is a time series object and types is a character vector indi-
cating the type of outliers to be considered by the procedure.
All five types of outliers described in the previous section The algorithm finds 14 outliers, with two being of transient
can be considered; AO, IO, LS, SLS and TC. If types is change (TC) type while the rest are additive outliers (AO).
not specified, then AO, LS, and TC are be considered by The outliers are shown in Figure 16;
default. The isolated sharp spikes represent AO, while a spike
that takes a few periods to disappear represents a TC. Note
Example 7. The DJ30 index consists of 30 stocks that are that the x−axis is labelled by the time (date) index, not
meant to reflect US market performance. Historical data the actual date. The latter can easily be extracted using the
of all stocks currently involved in the Dow Jones Industrial following code:
Average is available online on Kaggle [36]. For each of the
30 components of the index, there is one CSV file named ots=out_tso$outliers
by the stock’s symbol (e.g. AAPL for Apple Inc.). Each file [Link](type=ots$type,
provides historically adjusted market-wide data (daily, max. date=data[ots$ind,1])
5 years back).
After reading all the CSV files in R via [Link](), yields
the daily closing return for each component is sorted by type date type date
date and stored in a data frame named data. Outliers 1 TC 2015-01-28 8 AO 2017-08-01
are identified with tsoutliers’s tso() function. For 2 AO 2015-08-21 9 AO 2018-08-01
Apple Inc.’s stock, for instance, the outliers are found to be: 3 AO 2015-08-26 10 AO 2018-11-02
tso(ts(data$AAPL), 4 AO 2016-01-27 11 AO 2018-12-26
types=c("TC","AO","LS","IO","SLS")) 5 TC 2016-04-27 12 AO 2019-01-03
6 AO 2016-07-27 13 AO 2019-01-30
7 AO 2017-02-01 14 AO 2019-05-13
Outliers:
type ind time coefhat tstat 5.3 R Package: anomalize
1 TC 159 159 4.040 4.074 In this section, we will show how to use anomalize
2 AO 302 302 -6.116 -4.405 to detect anomalies in time series data. The package is
3 AO 305 305 5.735 4.130 available on CRAN, with the latest version always available
4 AO 410 410 -6.571 -4.732 on github.
5 TC 473 473 -4.521 -4.560 It is recommended to first install the package from CRAN
6 AO 536 536 6.496 4.678 (so that the dependencies are also installed locally), then
7 AO 666 666 6.098 4.392 update the package using devtools as shown below:
8 AO 791 791 6.629 4.774
9 AO 1043 1043 5.891 4.243 [Link](’anomalize’)
10 AO 1109 1109 -6.633 -4.777 library(devtools)
11 AO 1144 1144 7.042 5.072 install_github("business-science/anomalize")
12 AO 1149 1149 -9.961 -7.174 library(anomalize)
13 AO 1167 1167 6.833 4.921 library(tidyverse)
14 AO 1238 1238 -5.812 -4.186

[Link], [Link], [Link], [Link], [Link] (2020) Page 21 of 25

DATA SCIENCE REPORT SERIES ANOMALY DETECTION AND OUTLIER ANALYSIS

The anomalize() function is used to detect outliers in

a distribution with no trend or seasonality for tidy data,
and returns three columns:

remainder-l1 (lower limit for anomalies);

remainder-l2 (upper limit for anomalies, and
anomaly (Yes/No).
The first argument of the function is the “tibble” or “tbl_time“
object data; the second argument is the column target,
to which the function is applied; the third is the anomaly
detection method, either iqr geds.

The IQR method (inter-quartile range Q 3 − Q 1 ) is a gener-

alization of Tukey’s test (see Section 1). By defaults, the
limits are set at 3 times the IQR above Q 3 and below Q 1 Figure 17. Anomaly detection for Apple Inc.’s daily closing
(corresponding to α = 0.05); anything beyond those limits returns, computed with anomalize().
is considered to be an anomalous observation.
The alpha parameter can be adjusted; at α = 0.025, the
limits are 6 times the IQR above Q 3 and below Q 1 , making anomalize(remainder,
it more difficult for data to be an anomaly. Conversely, method="gesd",
α = 0.1 contracts the limits to 1.5 times the IQR above Q 3 alpha=0.05,
and below Q 1 , making it more likely that observations will max_anoms=0.2) %>%
be deemed anomalous. plot_anomaly_decomposition()
The IQR method does not depend on loops and is there-
fore fast and easily scaled, but it may not be as accurate in The output is shown in Figure 17, where the top plot dis-
detecting anomalies since the high leverage anomalies can plays the observations, the second and third plot displays
skew the centerline (median) of the IQR. the trend and seasonality components, respectively, and the
bottom plot displays the extracted data on which anomalies
The GESD method (generlized extreme studentized devi- are detected. The red markers show the anomalies found
ate test) progressively eliminates outliers using a Student’s by anomalize(). The recomposed series can also be
T −test comparing the test statistic to a critical value. Each plotted via time_recomposed():
time an outlier is removed, the test statistic is updated.
Once the test statistic drops below the critical value, all data_tb %>% time_decompose(AAPL) %>%
outliers are considered removed. anomalize(remainder) %>%
The α parameter adjusts the width of the critical values. time_recompose() %>%
plot_anomalies(time_recomposed=TRUE,
By default, α = 0.05. Because this method involves contin-
ncol=3,
uous updating via a loop, it is slower than the IQR method.
alpha_dots=0.5)
However, it tends to outperform IQR for outlier detection
and removal.
The output is shown in Figure 18. Anomalous points that
are shown in red can be extracted by the following code:
Other arguments include max_anoms (the maximum per-
centage of observations that can be identified as anomalies) anomalies=data_tb %>%
and verbose (boolean linked to the type of output). time_decompose(AAPL) %>%
anomalize(remainder) %>%
Example 8. In the previous sub-section, we used the func- time_recompose() %>%
tion tsoutliers() to detect outliers of daily closing filter(anomaly==’Yes’)
return for one of the components of DJ30, namely Apple
Inc. An alternative manner to detect anomalies in R is to The output is a time tibble with 16 rows (with date and
use anomalize(), as follows: observed as the first two variables); these are the 16
observations that anomalize() reports. Recall that in
data_tb=data %>% [Link]()
data_tb %>% previous section, tsoutliers() detected 14 outliers.
time_decompose(AAPL,
method="stl", It is interesting to compare the outlier/anomaly dates for
frequency=10, Apple Inc. data from both approaches. This can be done
trend="auto") %>% with the code on the following page:

[Link], [Link], [Link], [Link], [Link] (2020) Page 22 of 25

DATA SCIENCE REPORT SERIES ANOMALY DETECTION AND OUTLIER ANALYSIS

Figure 19. Comparing two methods of outlier/anomaly

detection for Apple Inc.’s daily closing returns.

shorter than a year, such as daily, weekly, monthly, or quar-

Figure 18. Anomaly detection for Apple Inc.’s daily closing terly while the later consists of longer term growth patterns.
returns with recomposed data; the grey portion represent the
expected (normal) trend. The first task in the anomaly detection of a time series
is thus to generate its remainders.
There are different ways to decompose a time series
odates=[Link](type=ots$type, to produce remainders: ARIMA and X12 are popular algo-
date=data[ots$ind,1]) rithms to do so.13
adates=[Link](date=anomalies$date,
observed=anomalies$observed) In general, high performance machine learning techniques
left_join(adates,odates,by="date") are not recommended for anomaly detection since the over-
fitting reduces the difference between the observed and
The list of detected outliers is found below: fitted values whereas in anomaly detection this difference
date observed type is essential to highlight the anomaly.
1 2015-01-28 5.653 TC On the other hand, seasonal decomposition performs
2 2015-08-11 -5.204 <NA> best for this task by removing the right features (i.e. sea-
3 2015-08-21 -6.116 AO sonal and trend components) while preserving the charac-
4 2015-08-26 5.735 AO teristics of anomalies in the remainders.
5 2016-01-22 5.317 <NA>
6 2016-01-27 -6.571 AO Finally, we note that there are other popular R outlier analy-
7 2016-04-27 -6.258 TC sis packages in R, such as AnomalyDetection, which
8 2016-07-27 6.496 AO uses a method similar to anomalize’s GESD. Interested
9 2017-02-01 6.098 AO readers are invited to try this package and compare their
10 2017-08-01 6.629 AO results with the functions reviewed in this section.
11 2018-08-01 5.891 AO
12 2018-11-02 -6.633 AO 6. Project
13 2018-12-26 7.042 AO
14 2019-01-03 -9.961 AO There is no substitute for practice: for the accompanying
15 2019-01-30 6.833 AO project, we ask you to test the performance and limitations
16 2019-05-13 -5.812 AO of four outlier detection algorithms (LOF, IsolationForest,
kNN, and PCA) on five datasets:
All the observations that have been detected by the original
approach (tsoutliers()) match the ones reported by the Airline dataset;
this new approach (anomalize()), but anomalize() the Distracted Driving Fatality dataset;
also detects two extra points, for which the type cannot be the House Prices dataset;
specified. The comparison is also illustrated in Figure 19. the Melbourne Temperature dataset, and
the Sale Transactions dataset.
5.4 Summary
For time series data, anomaly detection is usually performed Consult the project statement for a series of guided steps.
on time series remainders, where both the seasonal and
the trend components were removed; the former is the pres- 13
tsoutliers() uses ARIMA while anomalize() uses seasonal
ence of variations that occur at specific regular intervals decomposition.

[Link], [Link], [Link], [Link], [Link] (2020) Page 23 of 25

DATA SCIENCE REPORT SERIES ANOMALY DETECTION AND OUTLIER ANALYSIS

[19]
References Orchard, T., Woodbury, M. [1972], A Missing Infor-
[1]
mation Principle: Theory and Applications, Berkeley
Aggarwal, C.C. [2017], Outlier Analysis (2nd ed.), Symposium on Mathematical Statistics and Probability,
Springer. University of California Press.
[2]
Aggarwal, C. C. et YU, P. S. (2001). Outlier detection [20]
Torgo, L. [2017], Data Mining with R (2nd edition),
for high dimensional data. In ACM Sigmod Record, CRC Press.
pages 37–46. ACM. [21]
[3]
Chandola, V., Banerjee, A., Kumar, V. [2007], Outlier de-
Maimon, O., Rokach, L. [2010], Data Mining and tection: a survey, Technical Report TR 07-017, Depart-
Knowledge Discovery Handbook, Springer. ment of Computer Science and Engineering, University
[4]
Prasanta, G., et al. [2011], A Survey of Outlier De- of Minnesota.
tection Methods in Network Anomaly Identification. [22]
Hodge, V., Austin, J. [2004], A survey of outlier detec-
Oxford University Press, 54 (4), 570–588 tion methodologies, [Link]., 22(2):85-126.
[5]
Ranga, S.N.N.R., [2019], Outlier Detection: Tech- [23]
Height Percentile Calculator, by Age and Country,
niques and Applications: A Data Mining Perspective. on [Link]
Springer Nature Switzerland AG; 1st ed
[24]
[6]
Hariri, S., Kind, M. C., Brunner, R. J. [2018], Extended
Manish, G., et al., [2014], Outlier Detection for Tem- Isolation Forest, Computing Research Repository.
poral Data, IEEE Transactions on Knowledge and Data
[25]
Engineering 26(9), 2250–2267. Campello, R., Moulavi, D., Sander, J. [2013], Density-
[7] Based Clustering Based on Hierarchical Density Esti-
Kandanaarachchi, S., Hyndman, R.J. [2019], Dimen-
mates, Advances in Knowledge Discovery and Data
sion reduction for outlier detection using DOBIN,
Mining, Springer Berlin Heidelberg: 160–172.
Monash Business School.
[26]
[8] How HDBSCAN Works [2016] McInnes, L., Healy, J.,
Priyanga, D.T., et al., [2019], Anomaly detection in
Astels, S.
high-dimensional data.
[27]
[9] Findley, D.F., Hood, C.C., X-12-ARIMA and its Applica-
Aggarwal, C.C., Sathe, S. [2017], Outlier Ensembles,
tion to Some Italian Indicator Series, U.S. Bureau of
an Introduction, Springer.
the Census.
[10]
Aurore Archimbaud (2018). Détection non-supervisée [28]
Findley, D.F., Monsell, B.C., Bell, Otto and Chen [1998],
d’observations atypiques en contrôle de qualité : un
New Capabilities and Methods of the X-12-ARIMA Sea-
survol. Journal de la Société Française de Statistique,
sonal Adjustment Program, U.S. Bureau of the Census.
Vol. 159 No. 3 1-39
[29]
[11] An Introductory Course on Time Series Analysis, Aus-
Badr, W., 5 ways to detect outliers that every data
tralian Bureau of Statistics.
scientist should know, [Link]
[30]
[12] Seasonal Adjustment of Economic Time Series, Singa-
Mehrotra, K.G., Mohan, C.K., Huang, H. [2017],
pore Department of Statistics.
Anomaly Detection Principles and algorithms, Springer.
[31]
[13] [Link], [Link], Seasonal Adjustment Using The
Arora, L. [2019], An Awesome Tutorial to Learn Out-
X12 Procedure, SAS Institute.
lier Detection in Python using PyOD Library, on
[32]
Analitics Vidhya. Matthews Correlation Coefficient (MCC) on
[14] Wikipedia.
Santoyo, S. [2017], A Brief Overview of Outlier De-
[33]
tection Techniques, on [Link]. Principal Component Analysis on Wikipedia
[15] [34]
He, Z., Deng, S., Xu, X. [2005], A Unified Subspace Ruey S. Tsay [1988], level shifts, and variance changes
Outlier Ensemble Framework for Outlier Detection, in time series.
Advances in Web Age Information Management. [35]
Chen, C. and Liu, Lon-Mu (1993). Joint Estimation
[16]
Lazarevic, A., Kumar, V. [2005], Feature Bagging for of Model Parameters and Outlier Effects in Time Se-
Outlier Detection, ACM KDD Conference. ries. Journal of the American Statistical Association,
[17]
Fei, T.L., Ting, K.M., Zhou, Z.H. [2008], Isolation Forest, 88(421), pp. 284-297. doi: 10.2307/2290724
[36]
2008 Eighth IEEE International Conference on Data EOD data for all Dow Jones stocks
Mining: 413–422. [37]
Zhang, J., Lou, M., Ling, T. W. et Wang, H.(2004). Hos-
[18]
Ester, M., Kriegel, H.P., Sander, J., Xu, X. [1996], A miner : a system for detecting outlyting subspaces of
density-based algorithm for discovering clusters in high-dimensional data. In Proceedings of the Thirtieth
large spatial databases with noise, AAAI Press: 226– International Conference on Very Large Data Bases,
231. volume 30, pages 1265–1268. VLDB Endowment

[Link], [Link], [Link], [Link], [Link] (2020) Page 24 of 25

DATA SCIENCE REPORT SERIES ANOMALY DETECTION AND OUTLIER ANALYSIS

[38]
Zimek, A., Kriegel, H.-P., Kröger, P., Schubert, E. (2009).
Outlier detection in axis-parallel subspaces of high di-
mensional data. In Pacific-Asia Conference on Knowl-
edge Discovery and Data Mining, pages 831–838.
Springer.
[39]
Kriegel H.-P. et al. (2000) "LOF: Identifying density-
based local outliers," in Proceedings of the ACM SIG-
MOD International Conference on Management of Data
(ACM, New York), pp. 93-104
[40]
Müller, E., Assent, I., Iglesias S, P., Mulle, Y. et Bohm,
K. (2012). Outlier ranking via subspace analysis in
multiple views of the data. In IEEE 12th International
Conference on Data Mining (ICDM), pages 529-538.
IEEE.
[41]
Müller, E., Assent, I., Steinhausen, U. et Seidl T. (2008).
OutRank : ranking outliers in high dimensional data.
In ICDEW 2008, IEEE 24th International Conference
on Data Engineering Workshop, pages 600–603. IEEE
[42]
Müller, E., Schiffer, M. et Seidl, T. (2010b). Adaptive
outlierness for subspace outlier ranking. In Proceedings
of the 19th ACM International Conference on Informa-
tion and Knowledge Management, pages 1629–1632.
ACM.
[43]
Müller, E., Schiffer, M. et Seidl T. (2011). Statistical
selection of relevant subspace projections for outlier
ranking. In IEEE 27th International Conference on Data
Engineering (ICDE), pages 434–445. IEEE.
[44]
Jian T., Zhixiang C., Ada W. F., David, W C., Capa-
bilities of outlier detection schemes inlarge datasets,
framework and methodologies, Springer-Verlag Lon-
don Limited 2006
[45]
ROC curve
[46]
He, Z., Xu, X., Deng, S. [2005], A Fast Greedy Algo-
rithm for Outlier Mining.
[47]
Hastie, T., Leukemia dataset.
[48]
Leduc, O., Macfie, A., Maheshwari, A., Pelletier, M.,
Boily, P. [2019], Feature Selection and Dimension Re-
duction, Data Science Report Series, Data Action Lab
blog.
[49]
Leduc, O., Boily, P. [2019], Boosting with AdaBoost
and Gradient Boosting, Data Action Lab blog.
[50]
Lei, X. [2020], Distributed LOF: Density-Sensitive
Anomaly Detection With MapReduce, on
[Link].
[51]
Krishna, G. [2020] Performing Real-time Anomaly
Detection using AWS, on [Link].
[52]
Baron, D. [2018], Outlier Detection, XXX Winter
School of Astrophysics on Big Data in Astronomy,
GitHub repository.
[53]
Baron, D. [2016], Outlier Detection Algorithm on
Galaxy Spectra, GitHub repository.

[Link], [Link], [Link], [Link], [Link] (2020) Page 25 of 25