Daphne: A tool for anomaly detection

Title: Daphne: A tool for Anomaly Detection Author: Carlos García Ling Advisor: Daniel Selva Valero Department: Mechanical and Aerospace Engineering Academic year: 2017/18 Degree in Mathematics Universitat Polit`ecnica de Catalunya Facultat de Matem`atiques i Estad´ıstica Degree in Mathematics Bachelor’s Degree Thesis Daphne: A tool for Anomaly Detection Carlos Garc´ıa Ling June 2018 Supervised by Daniel Selva Valero Department of Mechanical and Aerospace Engineering Cornell University, Ithaca, New York Acknowledgements I would like to thank my advisor, Daniel Selva, for helping me develop this project. I also thank Antoni Vir´os, for introducing me to the Daphne project and welcoming me to the Ithaca. I would like to thank all the sources of funding that have allowed me to do this project in the USA: the TFG scholarship from CFIS and Funcaci´o Cellex, the MOBINT scholarship form Generalitat de Catalunya, and Cornell University. Abstract This thesis presents a versatile tool for anomaly detection, which is presented as a new functionality for Daphne’s cognitive assistant. This tool allows the user to analyze any kind of time series data, univariate or multivariate, and search for non conforming patterns, i.e. anomalies. Several anomaly detection algorithms have been integrated to analyze either a concrete variable, or the multivariate series as a whole. Multiple algorithms allow the user to search for different kind of anomalies in the data, accounting on data properties (seasonality, correlation), or on problem requirements (isolated anomalies, clustered anomalies) Complementary features to anomaly detection were also implemented, as an attempt to help the user choose a specific algorithm or understand its output. These features include data properties extraction (Seasonality, Correlation), checking agreement between different algorithms and diagnosing the novel anomalies using a database provided by the user. To enhance user experience an interactive interface has been implemented. The main components of the interface are a plot and different functionalities. The plot displays the problem’s variables and the outputs from the anomaly detection algorithms, and also allows the user to select and obtain information about the data. The functionalities allow the user to: run and configure anomaly detection methods; upload data; execute the complementary features. The current implementation of the tool can be found at the following link: https://www.selva-research.com/anomalies/ Keywords Anomaly Detection, Anomaly diagnose, Cognitive Assistant, Unsupervised Anomaly Detection, Time Series Analysis, Multivariate Analysis 1 Contents 1 Introduction 4 2 Related work 5 2.1 Anomaly Detection Tools ................................... 5 2.2 Anomaly Detection Algorithms ................................ 6 3 Methodologies 7 3.1 Statistical Methods ...................................... 7 3.1.1 Proposed Method ................................... 7 3.2 ARIMA Model Based ..................................... 9 3.2.1 ARIMA Models .................................... 10 3.2.2 Seasonality and exogenous variables ......................... 11 3.2.3 Application to Anomaly Detection .......................... 11 3.3 Adaptive Kernel density-based ................................ 13 3.4 Isolation Forest ........................................ 16 4 Daphne integration 20 4.1 User Interface ......................................... 20 4.1.1 User Menu ....................................... 21 4.1.2 Anomaly Plot ..................................... 21 4.1.3 Question Bar ..................................... 22 4.1.4 Functionalities ..................................... 22 4.2 Algorithms questions ..................................... 23 4.2.1 Anomalies Spotter ................................... 23 4.2.2 Agreement Between algorithms ............................ 24 2 4.2.3 Anomaly isolation ................................... 25 4.3 Data Related questions .................................... 25 4.3.1 Correlation ....................................... 25 4.3.2 Seasonality ....................................... 27 4.4 Anomaly Diagnose ....................................... 28 4.4.1 Anomaly Database .................................. 28 4.4.2 Diagnose algorithm .................................. 28 5 Practical examples 30 5.1 Traffic data .......................................... 30 5.2 Satellite data ......................................... 32 6 Future Work 36 6.1 Algorithm additions ...................................... 36 6.2 Anomalies Diagnose ...................................... 36 6.3 Online integration ....................................... 37 7 Conclusions 38 References 39 3 Daphne - Anomaly Detection 1. Introduction In machine learning, anomaly detection refers to the identification of patterns of data that do not adjust to the expected behaviour [1]. This patterns are referred to as anomalies or outliers. Another definition for anomaly is: an observation which deviates enough from the others to suspect that it was generated by a different mechanism [2]. According to this, the anomalies can be an indicator that the system analyzed is not behaving in the expected way. As a result, detecting these patterns has multiple real life applications, by spotting crashes or malfunctions early, and preventing their consequences. One of these applications is fraud detection [3], this includes any kind of undesirable behaviours such as bank fraud or system intrusion. As an example, financial institutions monitor credit card transactions. Applying anomaly detection algorithms, they can spot changes on the normal transactions pattern [4], and this can warn them that the credit card might have been stolen, or that an illicit use of the card is taking place. Another broadly used branch of fraud detection is computer intrusion [5]. Companies survey data, such as account logins or information requests, to guarantee that confidential information is not being compromised. In many systems considering the behaviour of different variables can help detect and predict faults or overloads. As an example, in mechanical systems, detecting anomalies in bearings can help to predict their remaining useful life [6]. In satellite control systems, addressing control wheels [7] and sensor [8] malfunctions, can prevent catastrophic failures. In satellite power subsystems, anomalies must be monitored to guarantee overall system health [9], and each particular component’s [10]. In distribution systems, like natural gas networks [11], detecting anomalies in consumption might warn for a failure like a leak on the distribution system. In epidemiology, anomaly detection techniques are used to control different factors related with infectious deceases [12]. Anomalies in this factors might indicate events like the start of an outbreak. Early detection and diagnose of these incidents are fundamental for a fast action against them, preventing important consequences in global health. Anomaly detection techniques can be even applied to computer vision data. In particular, video surveillance data has been fed to these algorithms in order to detect infractions in crowded scenes [13] like sneaking into the subway [14]. This thesis presents a tool for anomaly detection based on Daphne cognitive assistant [15], motivated by the diversity of applications of anomaly detection and the lack of a free versatile and open source tool for this problem. Daphne’s anomaly detection tool presents an interface to analyze any kind of time series data, by performing anomaly detection, extracting data properties and understanding and diagnosing the detected anomalies. The structure of this thesis starts with a related work section 2, where other anomaly detection tools are presented, as well as the current approaches on anomaly detection methods. Sections 3and 4explain in detail the Anomaly detection algorithms used and the functionalities implemented respectively. To show the potential of the tool some practical examples are developed in 5. Section 6presents future work in order to improve the current tool. To end Section 7presents conclusions and remarks on the work done. 4 •θiare the moving average coefficients, and they model the influence of the previous errors calculating the series value on the next term. Both αiand θicoefficients are found iteratively by maximizing the likelihood that these parameters explain the series. The parameters are specific for each combination of p,d,q, parameters. In order to choose the adequate parameters first dis estimated by differentiating the series until this is stationary. To check the stationarity of the series the Augmented Dickey-Fuller test is performed [49]. Once dparameter is fixed, a grid search is performed to select pand q. The parameters are chosen attending to the AIC metric, which represents a trade-off between the complexity of the ARIMA model, and the goodness of fit. 3.2.2 Seasonality and exogenous variables To introduce further complexity the model can be accommodated to include the seasonal trend and also the influence of other variables. The resulting model expression is obtained by adding the seasonality terms to the model as follows: (1 − p X i=1 αiLi)(1 − P X i=1 AiLi·s)(1 −L)d(1 −Ls)Dxn= (1 + q X i=1 θiLi)(1 + Q X i=1 TiLi·s)n(6) In the previous equation (6), the new terms added correspond to the seasonality modeling. The s variable indicates the seasonality lag, that represents the period of the seasonality. The resulting model is denominated Seasonal Autoregressive Integrated Moving Average (SARIMA). The parameters P,Dand Q, are the equivalent of their correspondent lowercase parameters, they indicate analogous properties but for the seasonal properties. These are estimated similarly to p,dand q. The parameters Aiand Tiare fitted in an iterative process which maximizes the likelihood of the parameters. The fitting of the SARIMA model parameters is done simultaneously (parameters α,θ,A,T). In practice the statsmodels implementation of this method for python was used Also exogenous variables can similarly be added to the model (SARIMAX). The relations implied in the model are linear, in consequence the variables added must be correlated to the principal time series data analyzed. 3.2.3 Application to Anomaly Detection The main use of this model is to be used as a forecasting tool. Although it can be also used as an anomaly detection tool In Algorithm 2the process to implement anomaly detection is explained in detail. It consists on two main parts which are: building the model, and predicting a confidence interval for each point. The classification of the points on anomalous or non anomalous is quite simple, the anomalous data will lay outside the confidence interval while the non anomalous are the ones located in it. 11 Daphne - Anomaly Detection Algorithm 2 SARIMAXAD(X,α,S,ex) Input: X:input data, α: sensibility parameter, S: seasonality parameter, ex: exogenous variables Output: Detected Anomalies Anomalies = [] {initialize anomalies detected} [p,d,q,P,D,Q] = Find Model Parameters(X,S,ex) SARIMAXModel = SARIMAX(X,S,ex, [p,d,q,P,D,Q]) Fit Model(SARIMAXModel) for all (t,y)in Xdo CIy= Confidence interval(SARIMAXModel,t,α) if not yin CIythen Add (t,y)to Anomalies end if end for return Anomalies (a) SARIMAX Method applied to traffic data setting α= 0.10 (b) SARIMAX Method applied to traffic data setting α= 0.01 Figure 2: The shaded part of the plots represent the confidence intervals. A higher α,2a, rises the sensitivity by narrowing the confidence interval and detecting more anomalies as a result. On the other hand, lower αvalues, 2b, mean in a wider confidence interval and only detecting very anomalous data 12 The anomaly detection algorithm uses αin order determine the amplitude of the confidence interval, the meaning of this parameter is the uncertainty of the prediction, this means than when for example an α= 0.01 parameter is chosen, there is a 1% uncertainty of the variable not being in the on the output confidence interval. Translated into anomaly detection terms, a higher value of the parameter means that the method will be more sensitive to anomalies, detecting many but also generating false alarms, while a lower parameter will only detect extreme anomalies, usually overlooking less relevant anomalies (Figure 2). The rest of the arguments of the algorithm are the seasonal properties (s) and exogenous variables (ex), adding this instances, increases the specificity of the model to by providing the different patterns and relations between the data, and as a result adapting better to the data analyzed. 3.3 Adaptive Kernel density-based Density based algorithms is considered as a part of the nearest neighbor based anomaly detection techniques. The approach is quite simple, being the main assumption that normal data is located in higher density regions, while anomalous lays in low density regions [1]. The fist step for this kind of methods is to compute a density for each point of the dataset, to do so a distance must be defined. As a clarification, all the distances calculated in this section refer to the euclidean distances, although, this method is valid for any distance that verifies the definition. The original algorithm developed in this section is presented on [35], it consists on two main steps, which are: •Calculate the local density for each point •Calculate a Local Anomaly Score (LAS) First of all, data must be normalized in order to avoid underestimate or overestimate any variables, due differences in the magnitudes of the values. In this case, the normalization performed on the data consists on standardizing the variables, establishing zero mean and unit variance. The density is computed using a Parzen window estimate, or Kernel density Estimate (KDE). Given a sample of msamples obtain from the probability density p(x), the estimator is formulated as it follows: ˆp(x) = 1 m m X i=1 h−nK(x−xi h) (7) The term K(·) represents the kernel function, and h a width parameter that controls the smoothness of the estimator. In this approach, the kernel function considered is the Gaussian Kernel. As a result, the expression for the local density results on: ρ(xi) = 1 m−1X j∈{1,2,...,m}\{i} exp−d(xi−xj) ri2(8) 13 Daphne - Anomaly Detection The parameter ri, is the kernel width, which is not a fixed parameter, instead it is determined for each point. This determines how to weight the terms accounting on their distance to xi. Wider rimeans that, fixed a distance, the same point will have a greater weight. In anomaly detection it is useful to fix the width narrower in low density regions and wider in high density regions, in order to help highlighting the differences between the normal and the anomalous data. It will be computed as follows for each datapoint: ri=c[dk−max +dk−min +−dk(xi)] (9) In the equation above, crepresents the overall smoothing parameter, which is recommended to be between 0.5 and 1, in the case of density estimation. The term d(xi) is the average distance of the knearest neighbors. The terms dk−max and dk−min represent respectively the maximum and minimum of the set {d(xi)}. Lastly the term takes a small value, ensures the non zero value of the width parameters, guarantees existence of ρ(xi). Figure 3: Using an anomaly threshold based only in the density of the points, point p1might be always flagged as an anomaly, to classify p2as anomalous this will imply to classify normal instances from C1as anomalies too [1] Once the densities have been computed, the output of the algorithm is computed. The LAS is computed as quotient between the density of the point considered, and the density of the neighboring points. LAS(xi) = log1 kPj∈kNN(xi)ρ(xj) ρ(xi)(10) The LAS overcome the limitations of a density defined for the globally, by allowing to detect anomalies 14 in a contextual way, avoiding the miss classification problems showed in Figure 3 To classify a point as an anomaly its density should be significantly lower than its neighbors. The LAS gives a measure to this concept, as it is defined in equation 10, a higher LAS means the densities of the neighboring points are higher than the analyzed point’s, and, as a result, the data is characterized as more anomalous. The values of the LAS are not comprised between any predetermined values, as a result, fixing a threshold a priory might not be the most adequate option, although values lower than 0 are not consider as anomalies, as this means that the density of the point is greater than the ones surrounding it. The pseudocode of the implemented method is found in algorithm 3 Algorithm 3 AdaptiveKDB(X,k,c) Input: X:input data, k: number of nearest neighbors, c: smoothing parameter Output: Anomaly Score AnomalyScore = [] {initialize anomalies detected} D= Compute Distance Matrix(X) for all xiin Xdo KNNi= Compute Nearest Neighbors(D(xi)) dk(xi) = Compute Mean(D(KNNi,xi)) end for dmax = max(dk(xi)) dmin = min(dk(xi)) for all xiin Xdo ri= Compute kernel radius(dk(xi)){Use equation 9} ρi= Compute density(D(xi), ri){Use equation 8} end for for all xiin Xdo AnomalyScore[i] = Compute Anomaly Score(ρi,ρ){Use equation 10} end for return AnomalyScore In Algorithm 3the matrix Drepresents the distance matrix, in which each element D(i,j)represents the distance between elements i, and j(D(i,j)=d(xi,xj)). The computation of Dmatrix is the slowest part of the algorithm, and increases with complexity O(n2). To avoid this issue, in problems with big volumes of data, a subset of the data can be used to train the model. If there is a subset of the data that represents the normal behaviour of the data, it can be useful to consider this as the training set. In the implementation done, a random sampling is performed as determined by the user to reduce execution time when working with big data sets. Although online implementations are not considered yet, the algorithm can be used for online anomaly detection. In this case a training set must be fixed, and also special attention must be payed to equation (9), as the d(xi) might exceed dk−max , as a result ri, must be limited to a minimum value (c[dk−min +]) 15 Daphne - Anomaly Detection 3.4 Isolation Forest The last Anomaly Detection Algorithms is the Isolation Forest [50]. It is focused on multivariate data sets, but does not take into account the time dimension. The strategy of this algorithm is to randomly divide the data, by recursively partitioning it, using what is denominated as isolation trees. Once many of these trees are obtained, the data points are evaluated through these trees to obtain an Anomaly Score. The Isolation Forest methodology it is based in a property of the data called mass [51], assuming that normal data lays in regions with bigger mass, while anomalous data lay in regions with lower mass. One of the main advantages of using this property, is the ability to efficiently detect clustered anomalies. This algorithm is significantly different to standard anomaly detection algorithms, as instead of profiling the normal data, it directly the isolates anomalies. To clarify, some elements of the problem are defined. The input of the algorithm consists on a data set X={x0,x1...xn}, with each of the data points defined by mattributes: x= (xq1,xq2, ...xqm), being Q={q1,q2, ...qm}the data attributes . The algorithm processes the data and assigns each point an anomaly score AS(xi)∈[0, 1], 1 meaning that the data is very anomalous, and 0 meaning the data is not anomalous at all. This algorithm performs two different steps in order to obtain the anomaly score: •Create the Isolation Forest (Training stage) •Evaluate the Data through the Isolation Forest (Testing stage) (a) (b) Figure 4: Differences regarding number of partitions in anomalous data 4a regular data 4b [50] The training stage is responsible of creating the model, by recognizing the data distribution, and in particular, the regions where the data is more ”isolated”. This translates in the idea that when the data 16 is randomly recursively partitioned, the number of partitions to arrive to one that contains only a single instance is small(Figure 4a), in contrast, generic data requires a greater number to reach isolation(Figure 4b). These partitions are denominated Isolation Trees, which can be defined recursively as follows: Definition 3.1. Isolation Tree Given Ta node of an Isolation tree. Tcan be either an external node (exNode) with no child or and internal node formed by a test and two child nodes (Tl,Tr). A test is defined by an attribute q∈Qa split value p, the test divides the elements of Tin the nodes: Tl={x∈T|xq<p} and Tl={x∈T|xq>=p} Algorithm 4 iTree(X,e,l) Input: X:input data, e:current tree height, l-height limit Output: an iTree if e≥lor |X| ≤ 1then return exNode{Size :|X|} else randomly select q∈Q max: maximum value of attribute q in X min: maximum value of attribute q in X randomly choose p∈[min,max] Split Xin Xl={x∈X|xq<p}and Xr={x∈X|xq≥p} return InNode( Left :iTree(Xl,e+ 1, l), Right :iTree(Xr,e+ 1, l), SplitAttribute :q, SplitValue :p) end if The training phase consists on computing a reasonable number tof iTrees, in order to identify the isolation patterns of the data. By computing a big number of isolation trees and averaging them, the normal behaviour of the data can be profiled. The pseudocode (Algorithm 4), shows how a iTree is computed. To implement this function, two classes (InNode and ExNode) must be defined as stated in definition 3.1. In contrast with other anomaly detection algorithms, in which more data means a better fit, isolation trees work better when keeping the training data small, as large datasets reduce the capacity of isolating the data. As a result, in the training phase the iTrees are not defined using the whole data X, instead, a sub sample of the data X0⊂X||X0|=ψ, is obtained to fit each iTree. The subsampling size parameter ψ is fixed to a number, fixing ψ= 256 generally provides a good performance of the algorithm [50]. Once the iForest is computes, the testing stage is performed. In the testing stage the average path length is computed for each x∈X, in order to obtain the Anomaly Score. Definition 3.2. Given an instance xand an iTree T, the path length of xreferred to T(h(x)) is the number of edges that must be gone through in order to reach an external node. 17 Daphne - Anomaly Detection Figure 5: Convergence of the average path length[50] When increasing the number of generated isolation trees, the average number of partitions converges to a concrete value as shown in Figure 5. The interest of the algorithms lays in identifying elements with shorter paths, as a result, there is no point in calculating in detail the path of elements with the longest paths. When defining the iTree, a limit of height was fixed in order to avoid this computations. This height is defined as the average path of an iTree (log2(ψ)), and represents the lparameter from Algorithm 4. In order to estimate the longitude remaining when an external node is reached, a term based (c(s)) in the size of the external node (s) is added to the path length of the node. This term is computed based in equation 11 that represents the average path length of an unsuccessful search in a Binary Search Tree (BST) [52] of size n. The algorithm to compute the path length is detailed in Algorithm 5 c(n)=2H(n−1) −(2(n−1)/n) (11) In equation 11 the term H(i) represents the harmonic number which can be estimated by H(i) = ln i+γ, being γthe Euler-Mascheroni constant. Algorithm 5 pathLength(x,T,e) Input: x: a data point, T: an iTree , e: current path length (initialized to 0) Output: path length of x if Tis exNode then return e + c(T.Size) {c() defined in } end if q = T.SplitAttribute if xq<T.SplitValue then return pathLenght(x,T.left,e+ 1) else {xq≥T.SplitValue} return pathLenght(x,T.right,e+ 1) end if 18 Once the average path has been computed the anomaly score can be computed. The output score of this method is based on the average path length of a BST with size the sub sampling size: c(ψ). The score is defined by equation 12. In the Anomaly score equation the term E(h(x)) indicates the expected value of the path length of element x, which can be estimated by averaging the lengths over the trees computed E(h(x)) ≈¯ h(x) AS(x)=2−E(h(x)) c(ψ)= 2−¯ h(x) c(ψ)(12) Considering equation 12 the anomaly score can be easily interpreted. Data points xwith values AS(x)≤0.5 mean that ¯ h(x)≥c(ψ), i.e. the mean path length is greater than an the mean unsuccessful BST path, which implies that xis not characterized as anomalous. On the other hand, if xhas an anomaly score significantly greater than 0.5 this will indicate an anomalous point. Once all the stages of the algorithms have been defined, the pseudocode of iTree can be found at Algorithm 6 Algorithm 6 iForest(x,ψ,t) Input: x: a data point, ψ: subsampling size , t: number of iTrees computed Output: Anomaly Score l = Set Limit iTree Height (ψ) T = initialize iTree array for i=0 to i=t do Xsample = Subsample a ψsized data from x T = Append(T, iTree(Xsample, 0, l)) end for AnomalyScore = initialize Anomaly Score array for all yin xdo PathLenghts = initialize Path Lenghts array for all iT in Tdo PathLenghts = Append(PathLenghts, pathLenght(y,iT,0)) end for AnomalyScore = Append(AnomalyScore, Compute Anomaly Score (PathLenghts)){Use equation 12} end for return AnomalySore 19 Daphne - Anomaly Detection 4. Daphne integration The main objective of this project is to integrate these algorithms into one tool which is easy to use, that provides the user an understanding on why the data is classified as anomalous, and aids the diagnose of these anomalies. The integration has been done in JavaScript and in Python, the logic has been integrated in Python using the Django framework, while all the user interface was done using the Vue framework based on JavaScript. 4.1 User Interface Figure 6: DAPHNE Anomaly Detection Feature Interface. 1 User Menu, 2 Anomaly Plot, 3 Question Bar and 4 Functionalities The user interface uses Daphne’s [15] design and disposition as a template. The interface is divided in four parts that can be seen in Figure 6which are: •1. User Menu: Displays the available features •2. Anomaly Plot: Graphically represents the data and the anomalies detected, it also allows the user to select different regions or to obtain information about a concrete instance of the data. •3. Question Bar: Allow the user to ask questions to Daphne, either by writing it in the bar, or by speaking it to a microphone. 20 4.3.2 Seasonality The notion of seasonality used in this section is based on the developed in subsection 3.2.2, as the main objective of this question is to be able to adequate choose the seasonality parameter in a SARIMA model. To determine seasonality auto correlation function (ACF) is analyzed [54]. The ACF is defined as follows: ACF(t) = E[(X−µ)(Lt−µ)] σ2(15) The previous equation defines the ACF given a fixed lag t. The terms that appear the equation, µ and σ2, correspond respectively to the mean and variance of variable X. In theory, a process in which seasonality is present, the seasonality term appears as a single spike on the ACF, in practice, a local extrem appears in the function, this can be seen in Figure 10. Figure 10: Sample of a ACF plot: The local maximums at lags 12, 24, 36, ... suggest a 12 lag seasonality The method implemented to asses if there is seasonality in the data, and consists on analyzing the ACF searching local extreme. The peaks of the function are classified according to their relative amplitude, if the amplitude of the peak is higher than 0.25 times the lag analyzed, the lag is showed as possible for being a seasonality corresponding to it and if it is higher than 0.40 times is displayed as a lag with high probability of containing a seasonality. It must be also noted that is common, when there is seasonality, to find the multiples of the values also as likely seasonalities. This lags actually reinforce the idea of the first one having a seasonality. 27 Daphne - Anomaly Detection The seasonality question analyzes the lags until a maximum fixed by the user. The respond consists on a list of seasonalities ordered by likelyhood. An example of response is found in Figure 11 4.4 Anomaly Diagnose Anomaly detection aside, the main objective of Daphne is to help the user understand the anomalies detected by the implemented algorithms. To do so an anomaly diagnose feature has been implemented. This feature is based on a database that the user must provide. It works by comparing the anomalies detected to the anomalies present on the data base and returns a possible diagnose that helps the user to understand what is happening, and gives a path of actuation in order to fix the anomaly 4.4.1 Anomaly Database A database must be provided to Daphne in order to infer from it a consistent diagnose. The current algorithm needs to be fed a data set of anomalies in which the appear the value of the variables when the anomaly took place, the date and time when the it took place in a the timestamp column, the diagnose of the anomaly in a column denominated anomalyType and the action taken on the anomaly in an anomalyAction column. This database is uploaded through a functionality called Anomaly Database Loader, which extracts the information to feed the diagnose algorithm. The anomaly, database, with size lwill be denoted as ADB ={a1,a2, ..., al} 4.4.2 Diagnose algorithm In order to classify the anomalies the algorithm simply finds the anomaly that is most similar to the novel anomaly. To do so, the following steps are followed: •Normalization of the data: both the data from the anomaly database, and the detected anomalies are normalized. The normalization is done attending to the problem data. Both mean µand standard deviation σare extracted from the data, after, the data is normalized: ˆx=x−µ σ •Distance calculation: for every novel anomaly na the distance to all the anomalies of the database is computed d(na,ai) •Nearest neighbor selection: the anomaly ajwith minimum distance d(na,aj) = minai∈ADB d(na,ai), is selected and taken as reference for the novel anomaly. •Diagnose: once a reference anomaly is selected, the novel anomaly is classified attending to the type from aj, and the action suggested is also the one that corresponds to the reference anomaly. The implementation also considers the case when the variables of the data provided and the anomaly database don’t coincide, in this case the algorithm is executed taking into account only the common 28 variables, but shows a warning to the user. Once the algorithm has obtained the most similar anomaly on the database Daphne responds as shown in Figure 18 29 Daphne - Anomaly Detection 5. Practical examples Two datasets were used to test Daphne’s anomaly detection tool. The fist was extracted from the Numenta Anomaly Benchmark [44], and the other cosists on real data from the battery system of a cubesat [55] 5.1 Traffic data This sample dataset contains real traffic data from the Twin Cities Metro area in Minnesota. The original dataset contains 7 variables, but, due to incomplete data, only the variables that represent occupancy and speed have been chosen to the analysis. Also data have been re sampled to one hour frequency due to inconsistent time steps in the original data. As a result the data analyzed consists on a time series of 5 variables, 2 occupancy variables and 3 speed variables, with a sampling frequency of one hour. The fist step to the analysis consist on run the anomaly detection methods. Both univariate methods were executed to all variables with the same parameters. In the case of the Windowed Statistics the window parameter was fixed to w= 100. This number has been chosen intuitively in order to provide information from the previous 4 days. The SARIMA based algorithm was executed taking into account a 24 hour seasonality, this can be induced from an a priori knowledge of the data (the traffic patterns show a daily seasonality), or asking daphne about it. When this question is asked to any of the variables on the dataset, a similar response to the one displayed in Figure 11 is obtained, always with high probability in lag 24, as expected. Figure 11: Answer about seasonality of variable occupancy t4013. The seasonalities are displayed from most importance to less. Having 48 as a probable value reinfoces the idea that there is a seasonality in lag 24. A high probability in lag 12 also makes sense, the acf plot that corresponds to this variable is the one in Figure 10, that shows that 12 is a value to consider. 30 Figure 12: Adaptive Kernel Density Based method’s Anomaly Score. When Daphne is asked to classify the tree most anomalous points, a threshold is returned, and the points classified as anomalies are marked in the anomaly plot with a purple dot Multivariate methods were executed by keeping the default parameters, except in the case of the isolation forest, in which the sampling parameters was modified to consider enough training points. When analyzing the anomaly score of the Adaptive Kernel Density Based method, three clear anomalies are spotted (Figure 12). Because of this, Daphne was asked to classify the three most anomalous points using the score of the multivariate methods for both algorithms. Figure 13: When both methods are asked to classify the 3 most anomalous points there is a correspondence of 100%, as the anomalies detected are the same. All the anomalies are marked as coincident and as related according to definitions 4.1 and 4.2, respectively 31 Daphne - Anomaly Detection Also Daphne was asked whether there is agreement between both multivariate methods. The response returned (Figure 13) shows there is complete agreement between the methods, providing higher certainty about these points being anomalies. Once there is certainty enough that the marked data points are in fact anomalies, Daphne can be used to spot the variables that trigger, or are related to these anomalies. In particular, when the points in the surrounding the anomalies are selected, the user can ask which are the variables with the higher number of spotted anomalies on the selection (Figure 14). Once these variables are known the plot of the variable and the anomaly score can be overlapped to understand the behaviour of the system. Figure 14: When asked to show the most anomalous variables, Daphne returns the number of anomalies spotted in the selection, as well as the variable, and the anomaly detection method associated to the anomalies. If the speed t4013 is displayed, it is clear that the anomalies spotted are driven by this variable This procedure is an example of how, Daphne can be used to understand the behaviour of the data, spot the anomalies on the data, and understand the behaviour of the system in this anomalies. 5.2 Satellite data The original housekeeping data from the cubesat contained 14 variables, but, in order to simplyfy the analysis of the data only three of the variables were kept. The variables discharded had to do with the orientation of different parts of the satellite, such as the solar panels or the chasis. The battery level (mV), photo voltaic current (mA) and total System current (mA), were the remaining variables considered. 32 (a) SARIMA Method applied to total system current (b) SARIMA Method applied to battery Figure 15: Both methods detect anomalies o the data when this constant steps don’t adjust to the seasonal data pattern. In the case of the battery 15a also the introduced anomalies are detected, due to the seasonal model, the these anomalies induce false alarms As the datasets available on the internet correspond to normal functioning of the satellite mostly didn’t contained anomalous data. Analyzing the data some irregularities might be spotted at first glance, some sections of the data appear suspiciously constant, which might indicate overflow of data. This anomalies are not detected by the Windowed Statistic methodology, but can be detected using a SARIMA based algorithm (Figure 15). Figure 16: Windowed Statistic method applied to Battery variable. Using this method, only the introduced anomalies are detected, but not the data overflow anomalies. 33 Daphne - Anomaly Detection Although these overflow data anomalies are detected, the values do not correspond to anomalous behaviour of the system. In order to test the diagnose function of Daphne, two anomalies were introduced, by setting anomalous values at the battery variable, one corresponding to a low battery level and another corresponding to an overcharge on the batteries. The anomalies were successfully detected by the Windowed Statistic (Figure 16), SARIMA based (Figure 15b), and appear as peaks in the anomaly score of the Kernel Adaptive Density Based Algorithm (Figure 17) Because of the coincidence of a lot of the values of the data due to low resolution, the performance of the iForest algorithm is quite poor. When the 2 most anomalous points according to the iForest are compared with the 2 according to the Adaptive Kernel Density Based, only one anomaly coincides (Figure 13). Figure 17: Adaptive Kernel Density Based Algorithm Anomaly scores and Agreement with iForest. Two spikes appear on the anomaly score, corresponding to the anomalies introduced. To test the anomaly diagnose feature a sample database was uploaded, containing anomalies corresponding to low battery levels and overcharge of the battery. These anomalies are compared with the detected ones to diagnose them. The response of this diagnose is shown in Figure 18. 34 Figure 18: Anomaly Diagnose Response. A Nearest Neighbors algorithm matches the novel anomalies with the ones on the database, and similar anomalies are displayed to the user. The procedure followed shows how can inconsistencies in seasonality can be found, and how the anomaly diagnose feature results useful into understanding the anomalies, and taking decisions regarding the actuation that must be followed 35 Daphne - Anomaly Detection 6. Future Work The current version of DAPHNE’s Anomaly Detection feature is the first step in building a versatile tool to detect irregularities on any kind of data. This section is dedicated to present the possibilities of improvement for the current implementation. This might be done by adding new anomaly detection algorithms specific for other kinds of anomalies, by improving current functionalities such as the Anomaly Diagnose, or by adapting the tool for new needs such as Online Anomaly detection. 6.1 Algorithm additions Daphne’s anomaly detection tool mainly relies on unsupervised algorithms in order to detect anomalies on the input data. Deepening in this kind of methods, other algorithms might be useful to implement, in order to address different problems in a more specific way. When dealing with systems in which accumulative effects are relevant, is useful to use anomaly metrics as the cumulative sum (CUSUM), that do not focus that much on the deviation of the data from normal instances, but on the sum of these deviations over time. They can be applied to computer traffic data to detect flooding attacks [56], or to monitor power systems [10]. When the system variables are highly correlated, this CUSUM approach, can be performed together with Principal Components Analysis (PCA). This has been successfully applied to early detect the collapse of a hospital’s emergency department [57]. Another interesting direction to improve Daphne’s anomaly detection tool might be to consider the temporal dimension on multivariate methods. The current implemented approaches in Daphne, only consider the general distribution of the data, but no the temporal scope. To do so, Dynamic Time Warping (DTW) similarity measure [58] might be interesting to introduce with anomaly detection techniques [33]. This measure evaluates the similarity in shape between two sequences, and can be applied to multivariate sequences. DTW has been already applied to detect anomalies in Satellite Telemetry Data [59]. Numenta’s Anomaly detection tool [18] is based in HTM networks that mimic the processes that take place in the human brain [60]. The algorithm developed for anomaly detection [20], has been proven to significantly overpass other state of the art algorithms [61]. This algorithm is developed for univariate series, and it can be interesting to add it to the list of implemented methods. 6.2 Anomalies Diagnose The current approach consists on an empirical model, also called data driven models, these do not depend on specific introduced expert knowledge, but, instead only on the patterns of the data provided [62]. The anomaly diagnose algorithm implemented is quite simplistic and may induce error, as it only relays on the nearest neighbors of the novel anomaly. In order to improve the reliability of the model more complex classification algorithms might be interesting to implement. Combining Kernel Fisher Discriminant Analysis and the k-nearest neighbors can be used to improve the performance of fault diagnosis [63]. Also, using associated rules, might be an 36 [63] Zhi Bo Zhu and Zhi Huan Song. A novel fault diagnosis system using pattern classification on kernel FDA subspace. Expert Systems with Applications, 38(6):6895–6905, 2011. [64] Zhiwei Gao, C Cecati, and S X Ding. A Survey of Fault Diagnosis and Fault-Tolerant Techniques Part I: Fault Diagnosis. IEEE Transactions On Industrial Electronics, 62(6):3768 – 3774, 2015. [65] M Orchard and G Vachtsevanos. A Particle Filtering Approach for On-Line Fault Diagnosis and Failure Prognosis. Measurement And Control, 31(3-4):1–18, 2007. [66] Zirije Hasani. Robust Anomaly Detection Algorithms for Real-time Big Data. 2017 6th MEDITERRANEAN CONFERENCE ON EMBEDDED COMPUTING, (June):11–15, 2017. [67] Shen Yin, Xiangping Zhu, and Student Member. Intelligent Particle Filter and Its Application to Fault Detection of Nonlinear System. IEEE transactions on industrial electronics, 62(6):3852–3861, 2015. [68] Laura Rettig, Mourad Khayati, Philippe Cudre-Mauroux, and Michal Piorkowski. Online anomaly detection over Big Data streams. 2015 IEEE International Conference on Big Data (Big Data), pages 1113–1122, 2015. [69] Li Dong, Shulin LIU, and Hongli ZHANG. A method of anomaly detection and fault diagnosis with online adaptive learning under small training samples. Pattern Recognition, 64(May 2016):374–385, 2017. 43