Machine Learning: Detecting Anomalies with Python
Anomaly detection is a great tool in helping engineers identify unusual patterns and deviations in system performance, which could indicate potential failures or inefficiencies.
By proactively detecting anomalies, engineers can build and maintain reliable systems, optimize resource management, and ensure operational efficiency.
This is the first in a twopart series exploring machine learning for anomaly detection. The accompanying post is titled: Machine Learning with a Vibration Sensor.
In this post, we’ll go through the steps of building an anomaly detection model in Python using triaxial acceleration, orientation, and rotation data.
Topics covered include:
 Preprocessing and Feature Extraction
 The Raw Data
 Feature Representation
 Import Python Libraries
 Creating the Training Set
 Normalizing the Data
 Training
 Building the Model
 Model Architecture
 Adam Optimization
 Training the Model
 Testing/Evaluation
 Creating the Test Set
 Evaluating the Model
 Classifying Anomalies
 Summary
Preprocessing and Feature Extraction
The Raw Data
`lib.process.file` import decimal import idelib import endaq import pandas as pd def metadata(doc: idelib.dataset.Dataset) > pd.DataFrame: """ Retrieve metadata about a Dataset object :param doc: an idelib.dataset.Dataset to extract metadata from :return: metadata to extract from Dataset object such as filename, sensor type, recording length, etc. """ filename = doc.name t0, t1 = (t * decimal.Decimal('1e6') for t in (doc.sessions[0].firstTime, doc.sessions[0].lastTime)) recording_length = t1  t0 serial_number_id = str(doc.recorderInfo['RecorderSerial']) sensor_type = doc.recorderInfo['PartNumber'] data = { "Name": filename, "Sensor Type": sensor_type, "Recording Length (seconds)": recording_length, "Serial Number": serial_number_id } return pd.DataFrame([data]) def sensor_info(doc) > pd.DataFrame: """ Retrieve sensor info about a Dataset object :param doc: an idelib.dataset.Dataset to extract sensor information from :return: a table containing information of each sensor present in the Dataset object """ channel_table = endaq.ide.get_channel_table(doc).data return channel_table.drop(columns=['channel']) def get_sensor_data(doc, sensor: tuple[str, str]) > pd.DataFrame: """ :param doc: an idelib.dataset.Dataset to extract sensor data from :param sensor: the sensor to select, as a tuple consisting of (name, type). See process.sensor_info to retrieve options. :return: a table containing the collected data from the sensor """ channel_table = endaq.ide.get_channel_table(doc).data sensor_names = zip(channel_table['name'].values, channel_table['type'].values) found = False # check if sensor is present in recording for s in sensor_names: if s[0] == sensor[0] and s[1] == sensor[1]: found = True break if not found: raise ValueError(f"The sensor {sensor} not present in {doc.name}") sensor_data = channel_table.loc[(channel_table['name'] == sensor[0]) & (channel_table['type'] == sensor[1])].squeeze() return endaq.ide.to_pandas(sensor_data.channel, time_mode="datetime", tz="utc") `lib.transform.ts` import numpy as np import pandas as pd def resample_ts(data: pd.DataFrame, rate: float): """ Resample a dataframe to `rate` :param data: a pandas dataframe with timestamps as indices and observed values at each timestamp as columns :param rate: the rate to resample the dataframe to, in seconds :return: the resampled dataframe """ rate_str = f'{rate}s' # resample the dataframe to specified rate resampled_df = data.resample(rate_str).mean() return resampled_df def split_series(ts: np.ndarray, segment_length: int, pad=None) > np.ndarray: """ Split a time series into multiple time series of length `segment_length`. If last part of `ts` is shorter than `segment_length`, rest of series is padded with value `pad` (default is last value in series) :param ts: a np.ndarray of shape (n_timesteps, n_features) :param segment_length: the length of the shortened time series :return: a np.ndarray of shape (n_series, n_timesteps, n_features) """ start = 0 n = len(ts) data = [] while start < n: # check if what's left is less than segment_length if n  start < segment_length: pad_value = ts[1] if pad is None else pad last_ts = np.concatenate((ts[start:], [pad_value for _ in range(segment_length  (n  start))])) data.append(last_ts) break split_ts = ts[start: start + segment_length] data.append(split_ts) start += segment_length return np.array(data)
The above contains several functions used for preprocessing.
The raw data comes in the form of IDE files. For this test setup, we are using a simple structure (fan) as a representative example. We train on “good” fan data, then damage fan to see how good the model is at finding anomalies.
The data from the intact fan we will consider as normal. The files normal_1.ide
, normal_2.ide
, and normal_3.ide
are 10minute recordings from each of the fan’s speed settings: lowest (1), middle (2), highest (3).
I have written a more indepth blog regarding the test: Machine Learning with a Vibration Sensor.
Feature Representation
The model will be trained on multivariate time series containing triaxle acceleration data, orientation, and rotation. Since the sampling rates of the accelerometer (20.0 kHz) and gyroscope (0.2 kHz) are different for this sensor, we’ll need to resample the collected data from each sensor to the same frequency. For this example, we chose to resample the data to 100 Hz. The exact frequency chosen here is arbitrary, but reduces the dimensionality of the data without losing information that could significantly impact the performance of the model.
The input to the autoencoder is a matrix of shape (n_samples, n_time_steps, n_features)
, where each sample represents a time series consisting of n_timesteps=10
and n_features=10
(acceleration and rotation have 3 axes each, and orientation is measured in quaternions with 4 axes). As a result, each sample will represent a time series that spans 100 milliseconds.
Import Python Libraries
In a Juypter notebook, we first import several Python libraries we will use:
from typing import List, Tuple import endaq from sklearn.model_selection import train_test_split from sklearn.preprocessing import StandardScaler import tensorflow as tf from tensorflow.keras import layers import numpy as np import pandas as pd import matplotlib.pyplot as plt # local modules from lib.process import file as process from lib.ts import transform import ml.autoencoder as ae
Then, we’ll define a utility function to help extract data from multiple sensors:
def multi_ts(file, sensors: List[Tuple[str, str]], rate: float) > pd.DataFrame: """ A helper to extract multivariate time series of 1/rate Hz from an IDE file :param file: the file to extract data from :param sensors: the set of sensors to extract data from :param rate: the rate, in seconds, to resample the data to :return: a dataframe containing sensor data resampled to 1/rate """ if len(sensors) == 0: raise ValueError("Please provide at least one sensor to extract data from") doc = endaq.ide.get_doc(file) dataframes = [] for sensor in sensors: try: sensor_data = process.get_sensor_data(doc, sensor) # resample to 1/rate Hz data_resampled = transform.resample_ts(sensor_data, rate) dataframes.append(data_resampled) except ValueError: # ignore if sensor not present in recording continue # concatenate dataframes along columns return pd.concat(dataframes, axis=1)
Creating the Training Set
We then construct the training dataset:
# set parameters rate = 0.01 segment_length = 10 sensors_to_get = [ ('X (40g)', 'Acceleration'), ('Y (40g)', 'Acceleration'), ('Z (40g)', 'Acceleration'), ('X', 'Rotation'), ('Y', 'Rotation'), ('Z', 'Rotation'), ('X', 'Quaternion'), ('Y', 'Quaternion'), ('Z', 'Quaternion'), ('W', 'Quaternion'), ('Internal Pressure', 'Pressure') # get multivariate data from each of the normal files normal_1 = multi_ts("../data/recordings/normal_1.ide", sensors=sensors_to_get, rate=rate) normal_2 = multi_ts("../data/recordings/normal_2.ide", sensors=sensors_to_get, rate=rate) normal_3 = multi_ts("../data/recordings/normal_3.ide", sensors=sensors_to_get, rate=rate) # some columns could include NaNs/missing values, so fill missing data with column mean normal_1 = normal_1.fillna(normal_1.mean()) normal_2 = normal_2.fillna(normal_2.mean()) normal_3 = normal_3.fillna(normal_3.mean()) # now split each dataset into multiple time series X1 = transform.split_series(normal_1.to_numpy(), segment_length) X2 = transform.split_series(normal_2.to_numpy(), segment_length) X3 = transform.split_series(normal_3.to_numpy(), segment_length) # series for training data # get multivariate data from each of the normal files normal_1 = multi_ts("../data/recordings/normal_1.ide", sensors=sensors_to_get, rate=rate) normal_2 = multi_ts("../data/recordings/normal_2.ide", sensors=sensors_to_get, rate=rate) normal_3 = multi_ts("../data/recordings/normal_3.ide", sensors=sensors_to_get, rate=rate) # some columns could include NaNs/missing values, so fill missing data with column mean normal_1 = normal_1.fillna(normal_1.mean()) normal_2 = normal_2.fillna(normal_2.mean()) normal_3 = normal_3.fillna(normal_3.mean()) # now split each dataset into multiple time series X1 = transform.split_series(normal_1.to_numpy(), segment_length) X2 = transform.split_series(normal_2.to_numpy(), segment_length) X3 = transform.split_series(normal_3.to_numpy(), segment_length) # combine sets into one X = np.concatenate([X1, X2, X3], axis=0) # create train and validation set X_train, X_val = train_test_split(X, test_size=0.2, random_state=42)
It is a common practice in machine learning to split a dataset into a train and validation set. Here, the validation set is a subset (20% of the entire size) of the normal dataset that the model will not encounter during training. We use the validation set to evaluate whether the trained model accurately identifies unseen data points that come from the same distribution of the training set as normal.
Normalizing the Data
It is often good practice to normalize the data before training the model, so that all feature values are on a common scale. Normalization also reduces the magnitude of the features, which benefits the algorithm when it calculates gradients during training in order to minimize the loss function.
def normalize(data): """ Normalize the input data :param data: the multivariate time series to normalize, a matrix of shape (n_samples, n_timesteps, n_features) :return: the multivariate time series with features normalized """ n_samples, n_timesteps, n_features = data.shape # reshape to 2D for normalization data_reshaped = data.reshape(1, n_features) # create the scaler and transform data scaler = StandardScaler().fit(data_reshaped) data_normalized = scaler.transform(data_reshaped) # reshape data back to original shape return data_normalized.reshape(n_samples, n_timesteps, n_features) # normalize training and validation data X_train = normalize(X_train) X_val = normalize(X_val)
`ml.autoencoder` import numpy as np import tensorflow as tf from tensorflow.keras import layers, models import numpy as np import tensorflow as tf from tensorflow.keras import layers, models class Autoencoder(models.Model): def __init__(self, input_shape, latent_dim, encoder, decoder): """ A class for building autoencoders :param input_shape: the shape of the input :param latent_dim: the dimensions of the latent space :param encoder: the encoder model :param decoder: the decoder model """ super(Autoencoder, self).__init__() self.latent_dim = latent_dim self.input_shape = input_shape self.encoder = encoder self.decoder = decoder self.threshold = None def call(self, x): """ A method that passes the input through the autoencoder and returns the reconstructed input :param x: the input :return: the reconstructed input """ encoded = self.encoder(x) decoded = self.decoder(encoded) return decoded def reconstruction_errors(self, x, axis=0): """ A method that calculates the reconstruction error for each sample in x :param x: the data to calculate reconstruction errors for :param x: the axis to calculate reconstruction error along :return: an Ndimensional array containing reconstruction errors for the samples """ x_reconstructed = self.call(x) return np.mean((x  x_reconstructed) ** 2, axis=axis) def set_threshold(self, threshold): """ A method to set the threshold for which the Autoencoder classifies a data point as an anomaly :param threshold: the reconstruction error threshold :return: the new threshold """ self.threshold = threshold return self.threshold def classify(self, x, axis=0): """ A method to classify points in x as normal or anomalies :param x: the data to classify :param axis: the axis to calculation reconstruction error along :return: returns Ndimensional array of 1s and 0s where 1 is an inlier and 0 is an outlier """ if self.threshold is None: raise ValueError("Threshold is none. Set threshold first before classifying.") re_errors = self.reconstruction_errors(x, axis=axis) return (re_errors <= self.threshold).astype(int)
We then define the architecture for the model:
# define architecture for model encoder = tf.keras.Sequential([ layers.Input(shape=(X_train.shape[1], X_train.shape[2])), layers.LSTM(128, activation="softsign", return_sequences=True), layers.LSTM(32, activation="softsign", return_sequences=False) ]) decoder = tf.keras.Sequential([ layers.RepeatVector(X_train.shape[1]), layers.LSTM(32, activation="softsign", return_sequences=True), layers.LSTM(128, activation="softsign", return_sequences=True), layers.TimeDistributed(layers.Dense(X_train.shape[2])) ])
Note that the architecture will depend on the application the model is being used for. When developing a model, we will often test different architectures in an effort to curate a more accurate model.
Model Architecture
Here, the encoder has 2 LSTM layers with the softsign
activation function. The second LSTM layer is the latent layer, which has 32 nodes. The encoder compresses a sequence into a reduced context vector. The RepeatVector layer prepares the context vector to be the initial input into the decoder.
The decoder also has 2 LSTM layers activated with softsign
and then a final TimeDistributed Dense layer, which is used to ensure the output sequence is of the same shape as the input sequence.
Adam Optimization
After defining the architecture, we can now instantiate the model:
# instaniate model and configuration model = ae.Autoencoder(latent_dim=32, input_shape=(X_train.shape[1], X_train.shape[2]), encoder=encoder, decoder=decoder) model.compile(optimizer='adam', loss='mse')
The compile
method configures the model to use mean squared error as the objective function, and to use Adam optimization, an iterative algorithm that is often used to minimize the objective function when training neural networks.
Training the Model
Now, we train the model:
history = model.fit(X_train, X_train, epochs=300, validation_data=(X_val, X_val))
where the epochs parameter specifies the number of iterations to run the optimization algorithm for. On each training step, training and validation loss will be calculated and stored in the history
object. Once the model has been trained, we can plot the losses as shown below:
# plot training and validation loss plt.plot(history.history['loss'], label="Train Loss") plt.plot(history.history['val_loss'], label="Validation Loss") plt.title(f"Loss vs Epoch") plt.grid(True) plt.legend() plt.xlabel("Epoch") plt.ylabel("Loss (MSE)") plt.show()
Overall, for any good model, we should expect training and validation loss to decrease, signaling that the model is getting better at reconstructing the input (i.e. the model is “learning”).
Testing/Evaluation
Now, let’s see how the model performs when evaluated on data from the modified fans. First, we process the data from the files containing the test data.
Creating the Test Set
# get data from fan with blade broken broken_data = multi_ts("../data/recordings/broken.ide", sensors=sensors_to_get, rate=rate) # impute columns broken_data = broken_data.fillna(broken_data.mean()) # split into multiple time series X_broken = transform.split_series(broken_data.to_numpy(), segment_length) X_broken = normalize(X_broken) # get data from fan with mass attached mass_1 = multi_ts("../data/recordings/mass_1.ide", sensors=sensors_to_get, rate=rate) mass_2 = multi_ts("../data/recordings/mass_2.ide", sensors=sensors_to_get, rate=rate) mass_3 = multi_ts("../data/recordings/mass_3.ide", sensors=sensors_to_get, rate=rate) # some columns could include NaNs/missing values, so fill missing data with column mean mass_1 = mass_1.fillna(mass_1.mean()) mass_2 = mass_2.fillna(mass_2.mean()) mass_3 = mass_3.fillna(mass_3.mean()) # now split each dataset into multiple time series X_mass_1 = transform.split_series(mass_1.to_numpy(), segment_length) X_mass_2 = transform.split_series(mass_2.to_numpy(), segment_length) X_mass_3 = transform.split_series(mass_3.to_numpy(), segment_length) # combine sets into one and normalize X_test_1 = np.concatenate([X_mass_1, X_mass_2, X_mass_3], axis=0) X_test_1 = normalize(X_test_1) # get data from fan in closet closet_data = multi_ts("../data/recordings/closet.ide", sensors_to_get, rate=rate) # impute columns closet_data = closet_data.fillna(closet_data.mean()) # split into time series X_closet = transform.split_series(closet_data.to_numpy(), segment_length) X_closet = normalize(X_closet) # get data fan with brush obstructing blades brush_data = multi_ts("../data/recordings/brush.ide", sensors_to_get, rate=rate) # impute columns brush_data = brush_data.fillna(brush_data.mean()) # split into time series X_brush = transform.split_series(brush_data.to_numpy(), segment_length) X_brush = normalize(X_brush)
Defined below is a function that produces some summary metrics for how the model performs on data X
such as the average, standard deviation, of several quantiles of reconstruction errors that occur across samples. We should expect the values that come from the distribution of reconstruction errors for the modified fans to be considerably higher than the distribution of errors for the training and validation set.
Evaluating the Model
def summarize_metrics(model, X, decimal_places=5): """ Produce a few summary statistics from the trained model :param model: the trained model :param X: the data to calculate summary metrics for :param decimal_places: the number of decimal places to round output to :return: a dict containing summary statistics from model evaluated on data """ errs = model.reconstruction_errors(X, axis=(1, 2)) qs = [0.05, 0.25, 0.5, 0.75, 0.95] quantiles = np.quantile(errs, q=qs) mse = np.mean((X  model.call(X)) ** 2) std = np.std(errs) data = { "MSE": f"{mse:.{decimal_places}f}", "STD RE": f"{std:.{decimal_places}f}", "Min RE": f"{np.min(errs):.{decimal_places}f}", "Max RE": f"{np.max(errs):.{decimal_places}f}" } for i, q in enumerate(quantiles): data[f"{qs[i]}q"] = q return data
We can then produce a dataframe that displays these summary metrics for the train, validation, and test data:
s_training = summarize_metrics(model, X_train) s_val = summarize_metrics(model, X_val) s_test_1 = summarize_metrics(model, X_test_1) s_test_2 = summarize_metrics(model, X_broken) s_test_3 = summarize_metrics(model, X_closet) s_test_4 = summarize_metrics(model, X_brush) lst = [s_training, s_val, s_test_1, s_test_2, s_test_3, s_test_4] indices = ["Train", "Validation", "Mass Attached", "Balde Broken", "Inside Closet", "Brush"] summary_data = pd.DataFrame(lst, index=indices)
The above plot shows the quantiles of reconstruction errors for the different datasets. As displayed above, the distribution of reconstruction errors for the training and validation set are similar, while those in the test sets for the modified fans are considerably different. Specifically, the reconstruction errors on the test sets are significantly larger than those from training and validation. The model therefore shows desired behavior, as it is reconstructing the time series from the normal fan better than how it reconstructs data from the test sets.
Classifying Anomalies
With the trained model, we can treat the reconstruction error on a particular sample as an anomaly score that we can use to categorize data. For instance, here we could set a threshold such that when the reconstruction error falls above the threshold, the data point (i.e. a series spanning 10 timesteps or 100 milliseconds) as an anomaly and those that fall below as normal. We could also create a more sophisticated system where we set multiple different thresholds, such that anomalies are identified in different levels.
If we consider the case of a simple threshold, we want to set a threshold such that the model classifies unseen normal data as normal frequently (let’s say above 95%) and anomalous data (e.g. that comes from the broken fan) as normal rarely (below 5%). The threshold can be chosen in different ways such as using trialanderror or following mathematical properties of statistical distributions. For instance, if the reconstruction errors from the normal data followed the normal distribution, then taking the threshold to be 2 standard deviations above the average reconstruction error should result in 97.5% of data points from the intact fan being classified as normal.
Suppose here we take the threshold to be the 0.97quantile (or 97th percentile) of reconstruction errors from the training set. So, we should expect the number of points classified as normal from a dataset with normal data to be about 97%. We’ll then use that threshold to classify data points as 1 (normal) or 0 (anomaly).
# set threshold train_errors = model.reconstruction_errors(X_train, axis=(1, 2)) threshold = np.quantile(train_errors, q=0.97) model.set_threshold(threshold) # classify points and get proportion for each train_labels = model.classify(X_train, axis=(1, 2)) val_labels = model.classify(X_val, axis=(1,2)) # mass attached test_1_labels = model.classify(X_test_1, axis=(1, 2)) # blade broken test_2_labels = model.classify(X_broken, axis=(1, 2)) # inside closet test_3_labels = model.classify(X_closet, axis=(1, 2)) # brush test_4_labels = model.classify(X_brush, axis=(1, 2)) # create dataframe showing proportion classified as Normal indices = ["Train", "Validation", "Mass Attached", "Blade Broken", "Inside Closet", "Brush"] lst = [train_labels, val_labels, test_1_labels, test_2_labels, test_3_labels, test_4_labels] data = [{"Proportion Classified as Normal": (labels == 1).mean()} for labels in lst] table = pd.DataFrame(data, index=indices) table

Proportion Classified as Normal 

Train 
0.969933 
Validation 
0.911690 
Mass Attached 
0.000000 
Blade Broken 
0.000000 
Inside Closet 
0.000000 
Brush 
0.000000 
As shown above, the model performs relatively well at classifying normal data that it did not see during training (i.e. validation set) with 91.2% accuracy. On the other hand, the number of points it classifies as normal from the anomalous settings are all approximately 0.
Summary
By establishing a robust foundation for anomaly detection using Python, we have demonstrated the initial steps necessary for creating a machinelearning model tailored to identify irregularities in complex systems. In the 2nd post of this series (Machine Learning with a Vibration Sensor), we will delve deeper into the practical applications of these models, focusing on the development of a vibration sensorbased anomaly detection model to monitor equipment health and predict failures.