{'fraction_non_alphanumeric': 0.0331668914239895, 'fraction_numerical': 0.013886357837956045, 'mean_word_length': 4.356891839125342, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Anomaly Detection by Robust Statistics

Real data often contain anomalous cases, also known as outliers. These may spoil the resulting analysis but they may also contain valuable information. In either case, the ability to detect such anomalies is essential. A useful tool for this purpose is robust statistics, which aims to detect the outliers by first fitting the majority of the data and then flagging data points that deviate from it. We present an overview of several robust methods and the resulting graphical outlier detection tools. We discuss robust procedures for univariate, low-dimensional, and high-dimensional data, such as estimating location and scatter, linear regression, principal component analysis, classification, clustering, and functional data analysis. Also the challenging new topic of cellwise outliers is introduced.


INTRODUCTION
In real data sets it often happens that some cases behave differently from the majority. Such data points are called anomalies in machine learning, and outliers in statistics. Outliers may be caused by errors, but they could also have been recorded under exceptional circumstances, or belong to another population. It is very important to be able to detect anomalous cases, which may (a) have a harmful effect on the conclusions drawn from the data, or (b) contain valuable nuggets of information.
In practice one often tries to detect outliers using diagnostics starting from a classical fitting method. However, classical methods can be affected by outliers so strongly that the resulting fitted model may not allow to detect the deviating observations. This is called the masking effect. Additionally, some good data points might even appear to be outliers, which is known as swamping. To avoid these effects, the goal of robust statistics is to find a fit which is close to the fit we would have found without the outliers. We can then identify the outliers by their large 'deviation' (e.g. its distance or residual) from that robust fit.
First we describe some robust procedures for detecting anomalies in univariate location and scale, as well as in multivariate data and in the linear regression setting. For more details on this part see [1][2][3] . Next we discuss principal component analysis (PCA) and some available robust methods for classification, clustering, and functional data analysis. Finally we introduce the emerging research topic of detecting cellwise anomalies.

ESTIMATING UNIVARIATE LOCATION AND SCALE
As an example of univariate data, suppose we have five measurements of a length: 6.27, 6.34, 6.25, 6.31, 6.28 (1) and we want to estimate its true value. For this, one usually computes the sample mean x i which in this case equalsx = (6.27 + 6.34 + 6.25 + 6.31 + 6.28)/5 = 6.29 . Let us now suppose that the fourth measurement has been recorded wrongly and the data become 6.27, 6.34, 6.25, 63.1, 6.28 .
In this case we obtainx = 17.65, which is far off. Alternatively, we could also compute the median of these data. For this we sort the observations in (2) from smallest to largest: 6.25 6.27 6.28 6.34 63.10 .
The median is the middle value, here yielding 6.28, which is still reasonable. We say that the median is more robust against an outlier.
More generally, the location-scale model states that the n univariate observations x i are independent and identically distributed (i.i.d.) with distribution function F ((x−µ)/σ) where F is known. Typically F is the standard gaussian distribution function Φ. We then want to find estimates for the unknown center µ and the unknown scale parameter σ.
The classical estimate of location is the mean. As we saw above, the mean is very sensitive to aberrant values among the n observations. In general, replacing even a single observation by a very large value can change the mean completely. We say that the breakdown value 4,5 of the sample mean is 1/n, so it becomes 0% for large n. In general, the breakdown value is the smallest proportion of observations in the data set that need to be replaced to carry the estimate arbitrarily far away. A breakdown value of 0% is thus the worst possible. See 6 for precise definitions and extensions. The robustness of an estimator is also measured by its influence function 7 which measures the effect of a single outlier. The influence function of the mean is unbounded, which again illustrates that the mean is not robust.
For the general definition of the median, we denote the ith ordered observation as x (i) .
The median is defined as x ((n+1)/2) if n is odd and (x (n/2) + x (n/2+1) )/2 if n is even. Its breakdown value is about 50%, meaning that the median can resist almost 50% of outliers. This is the best possible breakdown value since the clean data need to be in the majority.
The situation for the scale parameter σ is similar. The classical estimator is the standard . Since a single outlier can already make s arbitrarily large, its breakdown value is 0%. For instance, for the clean data (1) above we have s = 0.035, whereas for the data (2) with the outlier we obtain s = 25.41 ! A robust measure of scale is the median of all absolute deviations from the median (MAD), given by The constant 1.4826 is a correction factor which makes the MAD consistent at gaussian distributions. The MAD of (2) is the same as that of (1), namely 0.044. We can also use the Q n estimator 8 , defined as with k = h 2 ≈ n 2 /4 and h = n 2 + 1. Here . . . rounds down to the nearest integer. This scale estimator is thus the first quartile of all pairwise distances between two data points.
The breakdown value of both the MAD and the Q n estimator is 50%.
Also the (normalized) interquartile range (IQR) can be used, given by IQR = 0.7413(Q 3 − Q 1 ) where Q 1 = x n/4 is the first quartile of the data and Q 3 = x 3n/4 is the third quartile.
The IQR has a simple expression but its breakdown value is only 25%, so it is less robust than the MAD and Q n .
The robustness of the median (and the MAD) comes at a price: at the gaussian model it is less efficient than the mean. Many robust procedures have been proposed that strike a balance between robustness and efficiency, such as location M-estimators 9 . They are defined implicitly as the solution of the equation for a given real function ψ . The denominatorσ is an initial robust scale estimate such as Q n . A solutionμ to (4) can be found by an iterative algorithm, starting from the initial location estimateμ (0) = median i (x i ). Popular choices for ψ are the Huber function ψ(x) = x min(1, c/|x|) and Tukey's bisquare function ψ(x) = x(1 − (x/c) 2 ) 2 I(|x| c). These M-estimators contain a tuning parameter c which needs to be chosen in advance. Also M-estimators for the scale parameter σ exist.
People often use rules to detect outliers. The classical rule is based on the z-scores of the observations, given by where s is the standard deviation of the data. More precisely, the rule flags x i as outlying if |z i | exceeds 2.5, say. But in the above example (2) (1), which equals 1.41. The z-score of the outlier is small because it subtracts the nonrobust mean (which was drawn toward the outlier) and because it divides by the nonrobust standard deviation (which the outlier has made much larger than in the clean data). Plugging in robust estimators of location and scale such as the median and the MAD yields the robust scores which yield a much more reliable outlier detection tool. Indeed, in the contaminated example Also Tukey's boxplot is often used to pinpoint possible outliers. In this plot a box is drawn from the first quartile Q 1 of the data to the third quartile Q 3 . Points outside the interval [Q 1 − 1.5 IQR, Q 3 + 1.5 IQR], called the fence, are traditionally marked as outliers.
Note that the boxplot assumes symmetry, since we add the same amount to Q 3 as what we subtract from Q 1 . At asymmetric distributions the usual boxplot typically flags many regular data points as outliers. The skewness-adjusted boxplot 10 corrects for this by using a robust measure of skewness 11 in determining the fence.

MULTIVARIATE LOCATION AND COVARIANCE ESTI-MATION
From now on we assume that the data are d-dimensional and are stored in an n × d data matrix X = (x 1 , . . . , x n ) T with x i = (x i1 , . . . , x id ) T the ith data point. Classical measures of location and scatter are given by the empirical meanx = 1 n n i=1 x i and the empirical covariance matrix S X = n i=1 (x i −x)(x i −x) T /(n − 1). As in the univariate case, both classical estimators have a breakdown value of 0%, that is, a small fraction of outliers can completely ruin them.
As an illustration we consider a bivariate dataset (from page 59 in 2 ) containing the logarithms of body weight and brain weight of 28 animal species, with scatterplot in Figure 1.
Any point x has a so-called Mahalanobis distance (or 'generalized distance') to the meanμ =x , taking the covariance matrixΣ = S X into account. The MD is constant on ellipsoids. The so-called 97.5% tolerance ellipsoid is given by M D(x) χ 2 d,0.975 where χ 2 d,0.975 is the 0.975 quantile of the chi-squared distribution with d degrees of freedom. In this bivariate example d = 2, and the resulting ellipse is drawn in red. We see that it is inflated in the direction of the three outliers 6, 16, and 26 which are dinosaurs having low brain weight and high body weight. As a result these data points fall near the boundary of the tolerance ellipse, i.e. their M D(x i ) are not very high. Alternatively we can compute robust estimates of location and scatter (covariance), for instance by the Minimum Covariance Determinant (MCD) method 12,13 . The MCD looks for those h observations in the data set (where the number h is given by the user) whose classical covariance matrix has the lowest possible determinant. The MCD estimate of locationμ is then the average of these h points, whereas the MCD estimate of scatterΣ is their covariance matrix, multiplied by a consistency factor. (By default this is then followed by a reweighting step to improve efficiency at gaussian data.) Instead of Mahalanobis distances we can then compute robust distances, again given by (7) but now with the robust estimatesμ andΣ .
This yields the robust tolerance ellipse shown in blue in Figure 1. This ellipse exposes the three dinosaurs, and we see two species near the upper boundary, 17 (rhesus monkey) and 14 (human). In dimension d = 4 or higher it becomes infeasible to visualize the tolerance ellipsoid, but we still have the distances. The distance-distance plot (DD-plot) in Figure 2 shows the robust distance RD(x i ) of each data point versus its classical Mahalanobis distance MD(x i ) .
The horizontal and vertical cutoff lines are at χ 2 d,0.975 and the dashed line is where classical and robust distances coincide. We see that the RD(x i ) flag all the outliers in this dataset, while the M D(x i ) don't. For a dataset in which they are very similar we can trust classical statistical methods, but when they differ much (like here) the DD-plot detects the outlying data points. This does not imply we should somehow delete them, but rather that they should be investigated and understood. Outliers are not necessarily 'errors': they can also correspond to unusual circumstances or be members of a different population.
The MCD estimator, as well as its weighted version, has a bounded influence function and breakdown value (n − h + 1)/n, hence the number h determines the robustness of the estimator. The MCD has its highest possible breakdown value when h = (n + p + 1)/2 .
When a large proportion of contamination is expected, h should thus be chosen close to 0.5n.
Otherwise an intermediate value for h, such as 0.75n, is recommended to obtain a higher finite-sample efficiency. Reference 14 gives a more detailed overview of the MCD estimator and its properties.
The computation of the MCD estimator is non-trivial and naively requires an exhaustive investigation of all h-subsets out of n. Fortunately a much faster algorithm was constructed, called FastMCD 15 . It starts by randomly drawing many p + 1 observations from the data set. Based on these subsets, h-subsets are obtained by means of so-called C-steps (see 15 for details). More recently an even faster algorithm called DetMCD was devised 16 which carries out a deterministic computation instead of random sampling.
The MCD assumes that n > d, so there must be more data points than dimensions, and it works best when n > 5d. When there are more than, say, 20 dimensions and/or d n other methods are needed. One is to compute robust principal components as described in a section below. Another is to use the minimum regularized covariance determinant (MRCD) Many other robust estimators of location and scatter have been presented in the literature.
The first such estimator was proposed by Stahel 18 and Donoho 19 (see also 21 ). They defined the so-called Stahel-Donoho outlyingness of a data point x i as where the maximum is over all directions (i.e., all d-dimensional unit length vectors u), and x T j u is the projection of x j on the direction u. In each direction this uses the robust z-scores (6). Recently a version of (8)

LINEAR REGRESSION
The multiple linear regression model assumes that there are d 'explanatory' x-variables as well as a response variable y which can be approximated by a linear combination of the x-variables. More precisely, the model says that for all data points (x i , y i ) it holds that where the errors ε i are assumed to be independent and identically distributed with zero mean and constant variance σ 2 . Applying a regression estimator to the data yields d + 1 regression coefficients, combined asβ = (β 0 , . . . ,β d ) T . The residual r i of case i is defined as the difference between the observed response y i and its estimated valueŷ i .
The classical least squares (LS) method to estimate β minimizes the sum of the squared residuals. It is popular because it allows to compute the regression estimates explicitly, and it is optimal if the errors have a gaussian distribution. Unfortunately LS is extremely sensitive to outliers, i.e. data points that do not obey the linear pattern formed by the majority of the data.
For instance, Figure 3 shows the Hertzsprung-Russell diagram of the star cluster CYG OB1, containing 47 stars. The x-coordinate of each star is the logarithm of its surface temperature, and the y-coordinate is the logarithm of its light intensity. Most of the stars belong to the so-called main sequence, whereas 11, 20, 30, 34 are giant stars and 7 is intermediate.
The least squares line is shown in red, and has a negative slope although the main sequence slopes upward. It has been pulled away by the leverage exerted by the four giant stars. As an unfortunate side effect, the giant stars do not have larger absolute residuals than some of the main sequence stars, so only looking at residuals would not allow to detect them.
where (r 2 ) (1) (r 2 ) (2) · · · (r 2 ) (n) are the ordered squared residuals. (They are first squared, and then ordered.) By not adding all the squared residuals the LTS makes it possible to fit the majority of the data, whereas the outliers can have large residuals. In We can then identify outliers by their large standardized LTS residuals r i /σ LTS . We can also use the standardized LTS residuals to assign a weight to every observation. The weighted LS estimator with these LTS weights inherits the nice robustness properties of LTS, but is more efficient and yields all the usual inferential output such as t-statistics, F-statistics, an R 2 statistic, and the corresponding p-values. Alternatively, inference for LTS can be based on the fast robust bootstrap proposed in 28,29 . In most situations we have more than one explanatory variable, and for dimension d = 3 and higher it is no longer possible to perceive the linear patterns by eye. It is in those cases that robust regression becomes the most useful. To flag and interpret the outliers we can use the outlier map of 30 which plots the standardized LTS residuals versus robust distances (7) based on (for instance) the MCD estimator applied to the x-variables only. Figure 4   The main benefit of the outlier map is when the data has more dimensions. For instance, the stackloss data 31 is a benchmark data set with 21 points with d = 3 explanatory variables, an intercept term and a response variable y i . We cannot easily interpret such 4-dimensional data, but we can still look at the outlier map in the right panel of Figure 5. We see that 4 is a vertical outlier, 1, 3, and 21 are bad leverage points, and 2 is a good leverage point. Note that the left panel of Figure 5 does not flag any of these points because it uses the classical LS residuals and the classical distances M D(x i ), both of which tend to mask atypical points. A first group of robust PCA methods is obtained by replacing the classical covariance matrix by a robust covariance estimator, such as the weighted MCD estimator or MMestimators 40,41 . Unfortunately the use of these covariance estimators is limited to small to moderate dimensions since they are not defined when d n.
A second approach to robust PCA uses Projection Pursuit techniques. These methods maximize a robust measure of spread to obtain consecutive directions on which the data points are projected, see 42,43 .
The ROBPCA 44 approach is a hybrid, which combines ideas of projection pursuit and robust covariance estimation. The projection pursuit part is used for the initial dimension reduction. Some ideas based on the MCD estimator are then applied to this lower-dimensional data space. In order to diagnose outliers we can draw an outlier map for PCA 44 , similar to the outlier map for regression in the previous section. A stylized example of such a PCA outlier map is shown in the right panel of Figure 6, which corresponds to the three-dimensional data in the left panel which is fitted by two principal components. On the vertical axis of the PCA outlier map we find the orthogonal distance of each data point to the PCA subspace. This is just the Euclidean distance of the data point to its projection. The orthogonal distance is highest for the points 3, 4, and 5 in the example. On the horizontal axis we see the score distance of each data point, which is just the robust distance (7) of its projection relative to all the projected data points. The score distance is rather high for the points 1, 2, 4, and 5 in the figure.
By combining both distance measures the outlier map allows to distinguish between four types of data points. Regular observations have both a small orthogonal distance and a small score distance. Points with a high score distance but a small orthogonal distance, such as points 1 and 2 in Figure 6, are called good leverage points as they can improve the accuracy of the fitted PCA subspace. Orthogonal outliers have a large orthogonal distance but a small score distance, like point 3. Bad leverage points have both a large orthogonal distance and a large score distance, like points 4 and 5. They lie far from the space spanned by the robust principal components, and after projection on that space they lie far from most of the other projected data. They are called 'bad' because they typically they have a large influence on classical PCA, as the eigenvectors will be tilted towards them.

OTHER MODELS Classification
The goal of classification, also known as discriminant analysis or supervised learning, is to obtain rules that describe the separation between known groups G j of d-dimensional data points, with an eye toward assigning new data points to one of the groups. We write p j for the membership probability, i.e. the probability for any observation to come from G j .
For low-dimensional data, a popular classification rule results from maximizing the Bayes posterior probability. At gaussian distributions this yields quadratic discriminant analysis (QDA), i.e. choosing the j for which x has the highest quadratic score d Q j (x) given by When all the covariance matrices are assumed to be equal, these scores can be simplified to where Σ is the common covariance matrix, yielding linear discriminant analysis (LDA).
Robust classification rules can be obtained by replacing the classical covariance matrices by robust alternatives such as the MCD estimator or S-estimators, as in [50][51][52][53] .
When the data are high-dimensional, this approach cannot be applied because the robust covariance estimators are no longer computable. One approach is to first apply robust PCA to the entire data set. Alternatively, one can also apply a PCA method to each group separately.
This is the idea behind the SIMCA (Soft Independent Modeling of Class Analogy) method 54 .
A robustification of SIMCA is obtained by first applying robust PCA to each group, and then constructing a classification rule for new observations based on their orthogonal distance to each subspace and their score distance within each subspace 55 .
An SVM classifier with an unbounded kernel, e.g. a linear kernel, is not robust and suffers the same problems as traditional linear classifiers. But when a bounded kernel is used, the resulting non-linear SVM classification handles outliers quite well 48 . As an alternative, one can apply KROBPCA combined with LDA to the scores 49 .

Clustering
Cluster analysis (also known as unsupervised learning) is an important methodology when handling large data sets. It searches for homogeneous groups in the data, which afterwards may be analyzed separately. Partitioning (non-hierarchical) clustering methods search for the best clustering in k groups.

Other applications
Robust statistics has many other uses apart from outlier detection. For instance, robust estimation can be used in automated settings such as computer vision 71,72 . Another aspect is statistical inference, such as the construction of robust hypothesis tests, p-values, confidence intervals, and model selection (e.g. variable selection in regression). This aspect is studied in 3,7 and in other works they reference.

DETECTING OUTLYING CELLS
Until recently people have always considered outliers to be cases (data points), i.e. rows of the n × d data matrix X. But recently the realization has grown that this paradigm is no longer sufficient for the high-dimensional data sets we are often faced with nowadays. Typically most data cells (entries) in a row are regular and a few cells are anomalous. The first paper to formulate the cellwise paradigm was 73 , which showed how such outlying cells propagate in computations. In more than a few dimensions, even a small percentage of outlying cells can spoil a large percentage of rows. This is fatal for rowwise robust methods, which require at least 50% of the rows to be clean.
Detecting cellwise outliers is a hard problem, since the outlyingness of a cell depends on the relation of its column to the other columns of the data, and on the values of the other cells in its row (some of which may be outlying themselves). The DetectDeviatingCells 74 algorithm addresses these issues, and apart from flagging cells it also provides a graphical output called a cellmap. As an example we consider the mortality by age for males in France from 1816 to 2010, obtained from http://www.mortality.org . Each row corresponds to the mortalities in a given calendar year. The left panel in Figure 8 shows the result of the ROBPCA method described in the section on principal components. Outlying rows are shown in black and regular rows in yellow. The analysis was carried out on the data set with the individual years and the individual ages, but as this resolution would be too high to fit on the page we have combined the cells into 5 × 5 blocks afterward. The combination of some black rows with some yellow ones has led to gray blocks. We can see that there were outlying rows in the early years, the most recent years, and during two periods in between. Note that a black row doesn't provide information about its cells.
By contrast, the result of DetectDeviatingCells in the right panel in Figure 8 identifies a lot more information. Cells with higher values than predicted are shown in red, and those with lower values in blue, after which the colors were averaged in the 5 × 5 blocks.
The outlying early years saw a high infant mortality. During the Prussian war and both world wars there was a higher mortality among young adult men. And in recent years mortality among middle-aged and older men has decreased substantially, perhaps due to medical advances. We also return to the glass data from the section on PCA. The top panel in Figure 9 shows the rows detected by the ROBPCA method. The lower panel is the cell map obtained by DetectDeviatingCells on this 180 × 750 dataset. After the analysis, the cells were again grouped in 5 × 5 blocks. We now see clearly which parts of each spectrum are higher/lower than predicted. The wavelengths of these deviating cells reveal the chemical elements responsible.
Ideally, after runningDetectDeviatingCells the user can look at the deviating cells and whether their values are higher or lower than predicted, and make sense of what is going on.
This may lead to a better understanding of the data pattern, to changes in the way the data are collected/measured, to dropping certain rows or columns, to transforming variables, to changing the model, and so on. If the data set is too large for visual inspection of the results, or the analysis is automated, the deviating cells can be set to missing after which the data set is treated by a method appropriate for incomplete data. A good rowwise robust method of this type is 75 .

SOFTWARE AVAILABILITY
All the examples in this paper were produced with the free software R 76 . The publicly available CRAN package robustbase contains Qn, covMcd, ltsReg, and lmrob, whereas rrcov has many robust covariance estimators, robust principal components, and robust LDA and QDA classification. ROBPCA and its extensions are available in rospca and robust SIMCA in rrcovHD. Robust clustering can be performed with the cluster and tclust packages. The package mrfDepth contains tools for functional data and cellWise provides cellwise outlier detection and cellmaps.
Matlab functions for many of these methods are available in the LIBRA toolbox 77,78 , which can be downloaded from http://wis.kuleuven.be/stat/robust .
The MCD and LTS methods are also built into S-PLUS as well as SAS (version 11 or higher) and SAS/IML (version 7 or higher).

CONCLUSIONS
We have surveyed the utility of robust statistical methods and their algorithms for detecting anomalous data. These methods were illustrated on real data, in frameworks ranging from covariance matrices, the linear regression model and principal component analysis, with references to methods for many other tasks such as supervised and unsupervised classification as well as the analysis of functional data. For high-dimensional data, sparse and regularized robust methods were developed recently.
We have described methods to detect anomalous cases (rowwise outliers) but also newer work on the detection of anomalous data cells (cellwise outliers). An important topic for future research is to further improve the efficiency of the robust methodologies, in terms of both predictive accuracy and computational cost.