osquery

In osquery, you can schedule queries to be run periodically by the osqueryd daemon, populating tables that you can then query for later inspection. For investigative purposes, you can also run queries in an ad hoc fashion by using the command-line interface, osqueryi. An example query that gives you a list of all users on the system is as follows:

SELECT * FROM users;

If you wanted to locate the top five memory-hogging processes:

SELECT pid, name, resident_size FROM processes
ORDER BY resident_size DESC LIMIT 5

Although you can use osquery to monitor system reliability or compliance, one of its principal applications is to detect behavior on the system that could potentially be caused by intruders. A malicious binary will usually try to reduce its footprint on a system by getting rid of any traces it leaves in the filesystem; for example, by deleting itself after it starts execution. A common query for finding anomalous running binaries is to check for currently running processes that have a deleted executable:

SELECT * FROM processes WHERE on_disk = 0;

Suppose that this query generates some data that looks like this:

2017-06-04T18:24:17+00:00        []
2017-06-04T18:54:17+00:00        []
2017-06-04T19:24:17+00:00        ["/tmp/YBBHNCA8J0"]
2017-06-04T19:54:17+00:00        []

A very simple way to convert this data into a numerical time series is to use the length of the list as the value. It should then be clear that the third entry in this example will register as an anomaly.

Besides tapping into system state, the osquery daemon can also listen in on OS-level events such as filesystem modifications and accesses, drive mounts, process state changes, network setting changes, and more. This allows for event-based OS introspection to monitor filesystem integrity and audit processes and sockets.

You can set up osquery by creating a configuration file that defines what queries and packs the system should use.⁵ For instance, you can schedule a query to run every half hour (i.e., every 1,800 seconds) that detects deleted running binaries by putting the following query statement in the configuration:

{
 ...
  // Define a schedule of queries to run periodically:
  "deleted_running_binary": {
    "query": "SELECT * FROM processes WHERE on_disk = 0;",
    "interval": 1800
  }
  ...
}

Note that osquery can log query results as either snapshots or differentials. Differential logging can be useful in reducing the verbosity of received information, but it can also be more complex to parse. After the daemon logs this data, extracting time series metrics is simply a matter of analyzing log files or performing more SQL queries on the generated tables.

Deep packet inspection

Deep packet inspection (DPI) is the process of examining the data encapsulated in network packets, in addition to the headers and footers. This allows for the collection of signals and statistics about the network correspondence originating from the application layer. Because of this, DPI is capable of collecting signals that can help detect spam, malware, intrusions, and subtle anomalies. Real-time streaming DPI is a challenging computer science problem because of the computational requirements necessary to decrypt, disassemble, and analyze packets going through a network intersection.

Bro is one of the earliest systems that implemented a passive network monitoring framework for network intrusion detection. Bro consists of two components: an efficient event engine that extracts signals from live network traffic, and a policy engine that consumes events and policy scripts and takes the relevant action in response to different observed signals.

One thing you can use Bro for is to detect suspicious activity in web applications by inspecting the strings present in the POST body of HTTP requests. For example, you can detect SQL injections and cross-site scripting (XSS) reflection attacks by creating a profile of the POST body content for a particular web application entry point. A suspicion score can be generated by comparing the presence of certain anomalous characters (the ' character in the case of SQL injections, and the < or > script tag symbols in the case of XSS reflections) with the baseline, which can be valuable signals for detecting when a malicious actor is attacking your web application.⁷

The set of features to generate through DPI for anomaly detection is strongly dependent on the nature of the applications that operate within your network as well as the threat vectors relevant to your infrastructure. If your network does not include any outward-facing web servers, using DPI to detect XSS attacks is irrelevant. If your network contains only point-of-sale systems connected to PostgreSQL databases storing customer data, perhaps you should focus on unexpected network connections that could be indicative of an attacker pivoting in your network.

If DPI is used in an environment with Transport Layer Security/Secure Sockets Layer (TLS/SSL) in place, where the packets to be inspected are encrypted, the application performing DPI must terminate SSL. DPI essentially requires the anomaly detection system to operate as a man-in-the-middle, meaning that communication passing through the inspection point is no longer end-to-end secure. This architecture can pose a serious security and/or performance risk to your environment, especially for cases in which SSL termination and reencryption of packets is improperly implemented. You need to review and audit feature generation techniques that intercept TLS/SSL traffic very carefully before deploying them in production.

Features for network intrusion detection

The Knowledge Discovery and Data Mining Special Interest Group (SIGKDD) from the Association of Computing Machinery (ACM) holds the KDD Cup every year, posing a different challenge to participants. In 1999, the topic was “computer network intrusion detection”, in which the task was to “learn a predictive model capable of distinguishing between legitimate and illegitimate connections in a computer network.” This artificial dataset is very old and has been shown to have significant flaws, but the list of derived features provided by the dataset is a good source of example features to extract for network intrusion detection in your own environment. Staudemeyer and Omlin have used this dataset to find out which of these features are most important;⁸ their work might be useful to refer to when considering what types of features to generate for network anomaly and intrusion detection. Aggregating transactions by IP addresses, geolocation, netblocks (e.g., /16, /14), BFP prefixes, autonomous system number (ASN) information, and so on can often be good ways to distill complex network captures and generate simple count metrics for anomaly detection.⁹

ARIMA

Using the ARIMA (autoregressive integrated moving average) family of functions is a powerful and flexible way to perform forecasting on time series. Autoregressive models are a class of statistical models that have outputs that are linearly dependent on their own previous values in combination with a stochastic factor.¹² You might have heard of exponential smoothing, which can often be equated/approximated to special cases of ARIMA (e.g., Holt-Winters exponential smoothing). These operations smooth jagged line charts, using different variants of weighted moving averages to normalize the data. Seasonal variants of these operations can take periodic patterns into account, helping make more accurate forecasts. For instance, seasonal ARIMA (SARIMA) defines both a seasonal and a nonseasonal component of the ARIMA model, allowing periodic characteristics to be captured.¹³

In choosing an appropriate forecasting model, always visualize your data to identify trends, seasonalities, and cycles. If seasonality is a strong characteristic of the series, consider models with seasonal adjustments such as SARIMA and seasonal Holt-Winters methods. Forecasting methods learn characteristics of the time series by looking at previous points and making predictions about the future. In exploring the data, a useful metric to learn is the autocorrelation, which is the correlation between the series and itself at a previous point in time. A good forecast of the series can be thought of as the future points having high autocorrelation with the previous points. ARIMA uses a distributed lag model in which regressions are used to predict future values based on lagged values (an autoregressive process). Autoregressive and moving average parameters are used to tune the model, along with polynomial factor differencing—a process used to make the series stationary (i.e., having constant statistical properties over time, such as mean and variance), a condition that ARIMA requires the input series to have.

In this example, we attempt to perform anomaly detection on per-minute metrics of a host’s CPU utilization.¹⁴ The y-axis of Figure 3-2 shows the percentage CPU utilization, and the x-axis shows time.

Figure 3-2. CPU utilization over time

We can observe a clear periodic pattern in this series, with peaks in CPU utilization roughly every 2.5 hours. Using a convenient time series library for Python, PyFlux, we apply the ARIMA forecasting algorithm with autoregressive (AR) order 11, moving average (MA) order 11, and a differencing order of 0 (because the series looks stationary).¹⁵ There are some tricks to determining the AR and MA orders and the differencing order, which we will not elaborate on here. To oversimplify matters, AR and MA orders are needed to correct any residual autocorrelations that remain in the differenced series (i.e., between the time-shifted series and itself). The differencing order is a term used to make the series stationary—an already stationary series should have a differencing order of 0, a series with a constant average trend (steadily trending upward or downward) should have a differencing order of 1, and a series with a time-varying trend (a trend that changes in velocity and direction over the series) should have a differencing order of 2. Let’s plot the in-sample fit to get an idea of how the algorithm does:¹⁶

import pandas as pd
import pyflux as pf
from datetime import datetime

# Read in the training and testing dataset files
data_train_a = pd.read_csv('cpu-train-a.csv',
    parse_dates=[0], infer_datetime_format=True)
data_test_a = pd.read_csv('cpu-test-a.csv',
    parse_dates=[0], infer_datetime_format=True)

# Define the model
model_a = pf.ARIMA(data=data_train_a,
                   ar=11, ma=11, integ=0, target='cpu')

# Estimate latent variables for the model using the
# Metropolis-Hastings algorithm as the inference method
x = model_a.fit("M-H")

# Plot the fit of the ARIMA model against the data
model_a.plot_fit()

Figure 3-3 presents the result of the plot.

Figure 3-3. CPU utilization over time fitted with ARIMA model prediction

As we can observe in Figure 3-3, the results fit the observed data quite well. Next, we can do an in-sample test on the last 60 data points of the training data. The in-sample test is a validation step that treats the last subsection of the series as unknown and performs forecasting for those time steps. This process allows us to evaluate performance of the model without running tests on future/test data:

> model_a.plot_predict_is(h=60)

The in-sample prediction test (depicted in Figure 3-4) looks pretty good because it does not deviate from the original series significantly in phase and amplitude.

Figure 3-4. In-sample (training set) ARIMA prediction

Now, let’s run the actual forecasting, plotting the most recent 100 observed data points followed by the model’s 60 predicted values along with their confidence intervals:

> model_a.plot_predict(h=60, past_values=100)

Bands with a darker shade imply a higher confidence; see Figure 3-5.

Figure 3-5. Out-of-sample (test-set) ARIMA prediction

Comparing the prediction illustrated in Figure 3-5 to the actual observed points illustrated in Figure 3-6, we see that the prediction is spot-on.

Figure 3-6. Actual observed data points

To perform anomaly detection using forecasting, we compare the observed data points with a rolling prediction made periodically. For example, an arbitrary but sensible system might make a new 60-minute forecast every 30 minutes, training a new ARIMA model using the previous 24 hours of data. Comparisons between the forecast and observations can be made much more frequently (e.g., every three minutes). We can apply this method of incremental learning to almost all the algorithms that we will discuss, which allows us to approximate streaming behavior from algorithms originally designed for batch processing.

Let’s perform the same forecasting operations on another segment of the CPU utilization dataset captured at a different time:

data_train_b = pd.read_csv('cpu-train-b.csv',
    parse_dates=[0], infer_datetime_format=True)
data_test_b = pd.read_csv('cpu-test-b.csv',
    parse_dates=[0], infer_datetime_format=True)

Forecasting using the same ARIMAX model¹⁷ trained on data_train_b, the prediction is illustrated in Figure 3-7.

Figure 3-7. Out-of-sample (test-set, data_train_b) ARIMAX prediction

The observed values are, however, very different from the predictions illustrated in Figure 3-8.

Figure 3-8. Actual observed data points (data_train_b)

We see a visible anomaly that occurs a short time after our training period. Because the observed values fall within the low-confidence bands, we will raise an anomaly alert. The specific threshold conditions for how different the forecasted and observed series must be to raise an anomaly alert is something that is highly application specific but should be simple enough to implement on your own.

Artificial neural networks

Another way to perform forecasting on time series data is to use artificial neural networks. In particular, long short-term memory (LSTM) networks^18,19 are suitable for this application. LSTMs are a variant of recurrent neural networks (RNNs) that are uniquely architected to learn trends and patterns in time series input for the purposes of classification or prediction. We will not go into the theory or implementation details of neural networks here; instead, we will approach them as black boxes that can learn information from time series containing patterns that occur at unknown or irregular periodicities. We will use the Keras LSTM API, backed by TensorFlow, to perform forecasting on the same CPU utilization dataset that we used earlier.

The training methodology for our LSTM network is fairly straightforward. We first extract all continuous length-n subsequences of data from the training input, treating the last point in each subsequence as the label for the sample. In other words, we are generating n-grams from the input. For example, taking n = 3, for this raw data:

raw: [0.51, 0.29, 0.14, 1.00, 0.00, 0.13, 0.56]

we get the following n-grams:

n-grams: [[0.51, 0.29, 0.14],
          [0.29, 0.14, 1.00],
          [0.14, 1.00, 0.00],
          [1.00, 0.00, 0.13],
          [0.00, 0.13, 0.56]]

and the resulting training set is:

The model is learning to predict the third value in the sequence following the two already seen values. LSTM networks have a little more complexity built in that deals with remembering patterns and information from previous sequences, but as mentioned before, we will leave the details out. Let’s define a four-layer²⁰ LSTM network:²¹

from keras.models import Sequential
from keras.layers.recurrent import LSTM
from keras.layers.core import Dense, Activation, Dropout

# Each training data point will be length 100-1,
# since the last value in each sequence is the label
sequence_length = 100

model = Sequential()

# First LSTM layer defining the input sequence length
model.add(LSTM(input_shape=(sequence_length-1, 1),
               units=32,
               return_sequences=True))
model.add(Dropout(0.2))

# Second LSTM layer with 128 units
model.add(LSTM(units=128,
               return_sequences=True))
model.add(Dropout(0.2))

# Third LSTM layer with 100 units
model.add(LSTM(units=100,
               return_sequences=False))
model.add(Dropout(0.2))

# Densely connected output layer with the linear activation function
model.add(Dense(units=1))
model.add(Activation('linear'))

model.compile(loss='mean_squared_error', optimizer='rmsprop')

The precise architecture of the network (number of layers, size of each layer, type of layer, etc.) is arbitrarily chosen, roughly based on other LSTM networks that work well for similar problems. Notice that we are adding a Dropout(0.2) term after each hidden layer—dropout²² is a regularization technique that is commonly used to prevent neural networks from overfitting. At the end of the model definition, we make a call to the model.compile() method, which configures the learning process. We choose the rmsprop optimizer because the documentation states that it is usually a good choice for RNNs. The model fitting process will use the rmsprop optimization algorithm to minimize the loss function, which we have defined to be the mean_squared_error. There are many other tunable knobs and different architectures that will contribute to model performance, but, as usual, we opt for simplicity over accuracy.

Let’s prepare our input:

...
# Generate n-grams from the raw training data series
n_grams = []
for ix in range(len(training_data)-sequence_length):
n_grams.append(training_data[ix:ix+sequence_length])

# Normalize and shuffle the values
n_grams_arr = normalize(np.array(n_grams))
np.random.shuffle(n_grams_arr)

# Separate each sample from its label
x = n_grams_arr[:, :-1]
labels = n_grams_arr[:, −1]
...

Then, we can proceed to run the data through the model and make predictions:

...
model.fit(x,
   labels,
   batch_size=50,
   nb_epochs=3,
   validation_split=0.05)

y_pred = model.predict(x_test)
...

Figure 3-9 shows the results alongside the root-mean-squared (RMS) deviation.

We see that the prediction follows the nonanomalous observed series closely (both normalized), which hints to us that the LSTM network is indeed working well. When the anomalous observations occur, we see a large deviation between the predicted and observed series, evident in the RMS plot. Similar to the ARIMA case, such measures of deviations between predictions and observations can be used to signal when anomalies are detected. Thresholding on the observed versus predicted series divergence is a good way to abstract out the quirks of the data into a simple measure of unexpected deviations.

Figure 3-9. Observed, predicted, and RMS deviation plots of LSTM anomaly detection applied on the CPU time series

Summary

Forecasting is an intuitive method of performing anomaly detection. Especially when the time series has predictable seasonality patterns and an observable trend, models such as ARIMA can capture the data and reliably make forecasts. For more complex time series data, LSTM networks can work well. There are other methods of forecasting that utilize the same principles and achieve the same goal. Reconstructing time series data from a trained machine learning model (such as a clustering model) can be used to generate a forecast, but the validity of such an approach has been disputed in academia.²³

Note that forecasting does not typically work well for outlier detection; that is, if the training data for your model contains anomalies that you cannot easily filter out, your model will fit to both the inliers and outliers, which will make it difficult to detect future outliers. It is well suited for novelty detection, which means that the anomalies are only contained in the test data and not the training data. If the time series is highly erratic and does not follow any observable trend, or if the amplitude of fluctuations varies widely, forecasting is not likely to perform well. Forecasting works best on one-dimensional series of real-valued metrics, so if your dataset contains multidimensional feature vectors or categorical variables, you will be better off using another method of anomaly detection.

Median absolute deviation

The standard deviation of a data series is frequently used in adaptive thresholding to detect anomalies. For instance, a sensible definition of anomaly can be any point that lies more than two standard deviations from the mean. So, for a standard normal dataset with a mean of 0 and standard deviation of 1, any data points that lie between −2 and 2 will be considered regular, while a data point with the value 2.5 would be considered anomalous. This algorithm works if your data is perfectly clean, but if the data contains outliers the calculated mean and standard deviations will be skewed.

The median absolute deviation (MAD) is a commonly used alternative to the standard deviation for finding outliers in one-dimensional data. MAD is defined as the median of the absolute deviations from the series median:²⁴

import numpy as np

# Input data series
x = [1, 2, 3, 4, 5, 6]

# Calculate median absolute deviation
mad = np.median(np.abs(x - np.median(x)))

# MAD of x is 1.5

Because median is much less susceptible than mean to being influenced by outliers, MAD is a robust measure suitable for use in scenarios where the training data contains outliers.

Grubbs’ outlier test

Grubbs’ test is an algorithm that finds a single outlier in a normally distributed dataset by considering the current minimum or maximum value in the series. The algorithm is applied iteratively, removing the previously detected outlier between each iteration. Although we do not go into the details here, a common way to use Grubbs’ outlier test to detect anomalies is to calculate the Grubbs’ test statistic and Grubbs’ critical value, and mark the point as an outlier if the test statistic is greater than the critical value. This approach is only suitable for normal distributions, and can be inefficient because it only detects one anomaly in each iteration.

Summary

Statistical metric comparison is a very simple way to perform anomaly detection, and might not be considered by many to be a machine learning technique. However, it does check many of the boxes for features of an optimal anomaly detector that we discussed earlier: anomaly alerts are reproducible and easy to explain, algorithms adapt to changing trends in the data, it can be very performant because of its simplicity, and it is relatively easy to tune and maintain. Because of these properties, statistical metric comparison might be an optimal choice for some scenarios in which statistical measures can perform accurately, or for which a lower level of accuracy can be accepted. Because of their simplicity, statistical metrics have limited learning capabilities, and often perform worse than more powerful machine learning algorithms.

Elliptic envelope fitting (covariance estimate fitting)

For normally distributed datasets, elliptic envelope fitting can be a simple and elegant way to perform anomaly detection. Because anomalies are, by definition, points that do not conform to the expected distribution, it is easy for these algorithms to exclude such outliers in the training data. Thus, this method is only minimally affected by the presence of anomalies in the dataset.

The use of this method requires that you make a rather strong assumption about your data—that the inliers come from a known analytical distribution. Let’s take an example of a hypothetical dataset containing two appropriately scaled and normalized features (e.g., peak CPU utilization and start time of user-invoked processes on a host in 24 hours). Note that it is rare to find time series datasets that correspond to simple and known analytical distributions. More likely than not, this method will be suitable in anomaly detection problems with the time dimension excluded. We will synthesize this dataset by sampling a Gaussian distribution and then including a 0.01 ratio of outliers in the mixture:

import numpy as np

num_dimensions = 2
num_samples = 1000
outlier_ratio = 0.01
num_inliers = int(num_samples * (1-outlier_ratio))
num_outliers = num_samples - num_inliers

# Generate the normally distributed inliers
x = np.random.randn(num_inliers, num_dimensions)

# Add outliers sampled from a random uniform distribution
x_rand = np.random.uniform(low=-10, high=10, size=(num_outliers, num_dimensions))
x = np.r_[x, x_rand]

# Generate labels, 1 for inliers and −1 for outliers
labels = np.ones(num_samples, dtype=int)
labels[-num_outliers:] = −1

Plotting this dataset in a scatter plot (see Figure 3-10), we see that the outliers are visibly separated from the central mode cluster:

import matplotlib.pyplot as plt

plt.plot(x[:num_inliers,0], x[:num_inliers,1], 'wo', label='inliers')
plt.plot(x[-num_outliers:,0], x[-num_outliers:,1], 'ko', label='outliers')
plt.xlim(-11,11)
plt.ylim(-11,11)
plt.legend(numpoints=1)
plt.show()

Figure 3-10. Scatter plot of synthetic dataset with inlier/outlier ground truth labels

Elliptical envelope fitting does seem like a suitable option for anomaly detection given that the data looks normally distributed (as illustrated in Figure 3-10). We use the convenient sklearn.covariance.EllipticEnvelope class in the following analysis:²⁶

from sklearn.covariance import EllipticEnvelope

classifier = EllipticEnvelope(contamination=outlier_ratio)
classifier.fit(x)
y_pred = classifier.predict(x)
num_errors = sum(y_pred != labels)
print('Number of errors: {}'.format(num_errors))

> Number of errors: 0

This method performs very well on this dataset, but that is not surprising at all given the regularity of the distribution. In this example, we know the accurate value for outlier_ratio to be 0.01 because we created the dataset synthetically. This is an important parameter because it informs the classifier of the proportion of outliers it should look for. For realistic scenarios in which the outlier ratio is not known, you should make your best guess for the initial value based on your knowledge of the problem. Thereafter, you can iteratively tune the outlier_ratio upward if you are not detecting some outliers that the algorithm should have found, or tune it downward if there is a problem with false positives.

Let’s take a closer look at the decision boundary formed by this classifier, which is illustrated in Figure 3-11.

Figure 3-11. Decision boundary for elliptic envelope fitting on Gaussian synthetic data

The center mode is shaded in gray, demarcated by an elliptical decision boundary. Any points lying beyond the decision boundary of the ellipse are considered to be outliers.

We need to keep in mind that this method’s effectiveness varies across different data distributions. Let’s consider at a dataset that does not fit a regular Gaussian distribution (see Figure 3-12).

Figure 3-12. Decision boundary for elliptic envelope fitting on non-Gaussian synthetic data

There are now eight misclassifications: four outliers are now classified as inliers, and four inliers that fall just outside the decision boundary are flagged as outliers.

Applying this method in a streaming anomaly detection system is straightforward. By periodically fitting the elliptical envelope to new data, you will have a constantly updating decision boundary with which to classify incoming data points. To remove effects of drift and a continually expanding decision boundary over time, it is a good idea to retire data points after a certain amount of time. However, to ensure that seasonal and cyclical effects are covered, this sliding window of fresh data points needs to be wide enough to encapsulate information about daily or weekly patterns.

The EllipticEnvelope() function in sklearn is located in the sklearn.covariance module. The covariance of features in a dataset refers to the joint variability of the features. In other words, it is a measure of the magnitude and direction of the effect that a change in one feature has on another. The covariance is a characteristic of a dataset that we can use to describe distributions, and in turn to detect outliers that do not fit within the described distributions. Covariance estimators can be used to empirically estimate the covariance of a dataset given some training data, which is exactly how covariance-based fitting for anomaly detection works.

Robust covariance estimates²⁷ such as the Minimum Covariance Determinant (MCD) estimator will minimize the influence of training data outliers on the fitted model. We can measure the quality of a fitted model by the distance between outliers and the model’s distribution, using a distance function such as Mahalanobis distance. Compared with nonrobust estimates such as the Maximum Likelihood Estimator (MLE), MCD is able to discriminate between outliers and inliers, generating a better fit that results in inliers having small distances and outliers having large distances to the central mode of the fitted model.

The elliptic envelope fitting method makes use of robust covariance estimators to attain covariance estimates that model the distribution of the regular training data, and then classifies points that do not meet these estimates as anomalies. We’ve seen that elliptic envelope fitting works reasonably well for a two-dimensional contaminated dataset with a known Gaussian distribution, but not so well on a non-Gaussian dataset. You can apply this technique to higher-dimensional datasets as well, but your mileage may vary—elliptic envelopes work better on datasets with low dimensionality. When using it on time series data, you might find it useful in some scenarios to remove time from the feature set and just fit the model to a subset of other features. In this case, however, note that you will not be able to capture an anomaly that is statistically regular relative to the aggregate distribution, but in fact is anomalous relative to the time it appeared. For example, if some anomalous data points from the middle of the night have features that exhibit values that are not out of the ordinary for a midday data point, but are highly anomalous for nighttime measurements, the outliers might not be detected if you omit the time dimension.

One-class support vector machines

We can use a one-class SVM to detect anomalies by fitting the SVM with data belonging to only a single class. This data (which is assumed to contain no anomalies) is used to train the model, creating a decision boundary that can be used to classify future incoming data points. There is no robustness mechanism built into standard one-class SVM implementations, which means that the model training is less resilient to outliers in the dataset. As such, this method is more suitable for novelty detection than outlier detection; that is, the training data should ideally be thoroughly cleaned and contain no anomalies.

Where the one-class SVM method pulls apart from the pack is in dealing with non-Gaussian or multimodal distributions (i.e., when there is more than one “center” of regular inliers), as well as high-dimensional datasets. We will apply the one-class SVM classifier to the second dataset we used in the preceding section. Note that this dataset is not ideal for this method, because outliers comprise one percent of the data, but let’s see how much the resulting model is affected by the presence of contaminants:²⁸

from sklearn import svm

classifier = svm.OneClassSVM(nu=0.99 * outlier_ratio + 0.01,
                             kernel="rbf",
                             gamma=0.1)
classifier.fit(x)
y_pred = classifier.predict(x)
num_errors = sum(y_pred != labels)
print('Number of errors: {}'.format(num_errors))

Let’s examine the custom parameters that we specified in the creation of the svm.OneClassSVM object. Note that these parameters are dependent on datasets and usage scenarios; in general, you should always have a good understanding of all tunable parameters offered by a classifier before you use it. To deal with a small proportion of outliers in the data, we set the nu parameter to be roughly equivalent to the outlier ratio. According to the sklearn documentation, this parameter controls the “upper bound on the fraction of training errors and the lower bound of the fraction of support vectors.” In other words, it represents the acceptable range of errors generated by the model that can be caused by stray outliers, allowing the model some flexibility to prevent overfitting the model to outliers in the training set.

The kernel is selected by visually inspecting the dataset’s distribution. Each cluster in the bimodal distribution has Gaussian characteristics, which suggests that the radial basis function (RBF) kernel would be a good fit given that the values of both the Gaussian function and the RBF decrease exponentially as points move radially further away from the center.

The gamma parameter is used to tune the RBF kernel. This parameter defines how much influence any one training sample has on the resulting model. Its default value is 0.5. Smaller values of gamma would result in a “smoother” decision boundary, which might not be able to adequately capture the shape of the dataset. Larger values might result in overfitting. We chose a smaller value of gamma in this case to prevent overfitting to outliers that are close to the decision boundary.

Inspecting the resulting model, we see that the one-class SVM is able to fit this strongly bimodal dataset quite well, generating two mode clusters of inliers, as demonstrated in Figure 3-13. There are 16 misclassifications, so the presence of outliers in the training data did have some effect on the resulting model.

Figure 3-13. Decision boundary for one-class SVM on bimodal synthetic data—trained using both outliers and inliers

Let’s retrain the model on purely the inliers and see if it does any better. Figure 3-14 presents the result.

Figure 3-14. Decision boundary for one-class SVM on bimodal synthetic data—trained using only inliers

Indeed, as can be observed from Figure 3-14, there now are only three classification errors.

One-class SVMs offer a more flexible method for fitting a learned distribution to your dataset than robust covariance estimation. If you are thinking of using one as the engine for your anomaly detection system, however, you need to pay special attention to potential outliers that might slip past detection and cause the gradual degradation of the model’s accuracy.

Isolation forests

Random forest classifiers have a reputation for working well as anomaly detection engines in high-dimensional datasets. Random forests are algorithmic trees, and stream classification on tree data structures is much more efficient compared to models that involve many cluster or distance function computations. The number of feature value comparisons required to classify an incoming data point is the height of the tree (vertical distance between the root node and the terminating leaf node). This makes it very suitable for real-time anomaly detection on time series data.

The sklearn.ensemble.IsolationForest class helps determine the anomaly score of a sample using the Isolation Forest algorithm. This algorithm trains a model by iterating through data points in the training set, randomly selecting a feature and randomly selecting a split value between the maximum and minimum values of that feature (across the entire dataset). The algorithm operates in the context of anomaly detection by computing the number of splits required to isolate a single sample; that is, how many times we need to perform splits on features in the dataset before we end up with a region that contains only the single target sample. The intuition behind this method is that inliers have more feature value similarities, which requires them to go through more splits to be isolated. Outliers, on the other hand, should be easier to isolate with a small number of splits because they will likely have some feature value differences that distinguish them from inliers. By measuring the “path length” of recursive splits from the root of the tree, we have a metric with which we can attribute an anomaly score to data points. Anomalous data points should have shorter path lengths than regular data points. In the sklearn implementation, the threshold for points to be considered anomalous is defined by the contamination ratio. With a contamination ratio of 0.01, the shortest 1% of paths will be considered anomalies.

Let’s see this method in action by applying the Isolation Forest algorithm on the non-Gaussian contaminated dataset we saw in earlier sections (see Figure 3-15):²⁹

from sklearn.ensemble import IsolationForest

rng = np.random.RandomState(99)

classifier = IsolationForest(max_samples=num_samples,
                             contamination=outlier_ratio,
                             random_state=rng)
classifier.fit(x)
y_pred = classifier.predict(x)
num_errors = sum(y_pred != labels)
print('Number of errors: {}'.format(num_errors))

> Number of errors: 8

Figure 3-15. Decision boundary for isolation forest on synthetic non-Gaussian data

Using isolation forests in streaming time series anomaly detection is very similar to using one-class SVMs or robust covariance estimations. The anomaly detector simply maintains a tree of isolation forest splits and updates the model with new incoming points (as long as they are not deemed anomalies) in newly isolated segments of the feature space.

It is important to note that even though testing/classification is efficient, initially training the model is often more resource and time intensive than other methods of anomaly detection discussed earlier. On very low-dimensional data, using isolation forests for anomaly detection might not be suitable, because the small number of features on which we can perform splits can limit the effectiveness of the algorithm.

Local outlier factor

The LOF is an anomaly score that you can generate using the scikit-learn class sklearn.neighbors.LocalOutlierFactor. Similar to the aforementioned k-NN and k-means anomaly detection methods, LOF classifies anomalies using local density around a sample. The local density of a data point refers to the concentration of other points in the immediate surrounding region, where the size of this region can be defined either by a fixed distance threshold or by the closest n neighboring points. LOF measures the isolation of a single data point with respect to its closest n neighbors. Data points with a significantly lower local density than that of their closest n neighbors are considered to be anomalies. Let’s run an example on a similar non-Gaussian, contaminated dataset once again:³¹

from sklearn.neighbors import LocalOutlierFactor

classifier = LocalOutlierFactor(n_neighbors=100)
y_pred = classifier.fit_predict(x)

Z = classifier._decision_function(np.c_[xx.ravel(), yy.ravel()])
Z = Z.reshape(xx.shape)

num_errors = sum(y_pred != labels)
print('Number of errors: {}'.format(num_errors))

> Number of errors: 9

Figure 3-16 presents the result.

Figure 3-16. Decision boundary for local outlier factor on bimodal synthetic distribution

As we can observe from Figure 3-16, LOF works very well even when there is contamination in the training set, and it is not very strongly affected by dimensionality of the data. As long as the dataset maintains the property that outliers have a weaker local density than their neighbors in a majority of the training features, LOF can find clusters of inliers well. Because of the algorithm’s approach, it is able to distinguish between outliers in datasets with varying cluster densities. For instance, a point in a sparse cluster might have a higher distance to its nearest neighbors than another point in a denser cluster (in another area of the same dataset), but because density comparisons are made only with local neighbors, each cluster will have different distance conditions for what constitutes an outlier. Lastly, LOF’s nonparametric nature means that it can easily be generalized across multiple different dimensions as long as the data is numerical and continuous.

Performance and scalability in real-time streaming applications

Many applications of anomaly detection in the context of security require a system that can handle real-time streaming classification requests and deal with shifting trends in the data over time. But unlike with ad hoc machine learning processes, classification accuracy is not the only metric to optimize. Even though they might yield inferior classification results, some algorithms are less time and resource intensive than others and can be the optimal choice for designing systems in resource-critical environments (e.g., for performing machine learning on mobile devices or embedded systems).

Parallelization is the classic computer science answer to performance problems. Parallelizing machine learning algorithms and/or running them in a distributed fashion on MapReduce frameworks such as Apache Spark (Streaming) are good ways to improve performance by orders of magnitude. In designing systems for the real world, keep in mind that some machine learning algorithms cannot easily be parallelized, because internode communication is required (e.g., simple clustering algorithms). Using distributed machine learning libraries such as Apache Spark MLlib can help you to avoid the pain of having to implement and optimize distributed machine learning systems. We further investigate the use of these frameworks in Chapter 7.