Anomaly Detection

1. Types of Anomalies

Understanding which type of anomaly you are looking for shapes both the algorithm choice and the feature engineering strategy.

Point Anomaly

A single observation is anomalous with respect to the rest of the data. This is the most common and well-studied case. Example: a single transaction of $50,000 when a customer normally spends $200/month. The transaction is unusual in isolation, regardless of any surrounding context.

Contextual Anomaly

An observation is anomalous only within a specific context, but would be considered normal in a different context. The same value can be normal or anomalous depending on when, where, or from whom it originates. Example: a temperature of 15°C is perfectly normal in spring but is anomalous in the middle of summer. The value itself is not unusual globally — only relative to the seasonal context. In fraud: a $500 restaurant charge is normal on a Friday night but suspicious at 3 AM on a Tuesday in a foreign country.

Collective Anomaly

A collection of observations is anomalous with respect to the entire data set, even though each individual observation may appear normal in isolation. Example: a series of small transactions — each under $10 — individually look like normal micro-purchases, but collectively, 200 such transactions in 10 minutes constitute structured fraud (card testing). In network intrusion, no single packet looks unusual, but the sequence of packets reveals a port scan.

2. Statistical Methods

Statistical methods are interpretable, fast, and appropriate for univariate, roughly normally distributed data. They are the right starting point before reaching for more complex algorithms.

Z-Score

Flag a point as anomalous if its standardized distance from the mean exceeds a threshold (typically 3):

$$z = \frac{x - \mu}{\sigma}, \quad \text{flag if } |z| > 3$$

Works well when data is approximately normal. Critical weakness: the mean and standard deviation themselves are sensitive to outliers. A single extreme value inflates $\sigma$, masking other outliers — the so-called "masking effect."

IQR Method

Non-parametric and robust to distributional assumptions. Compute the interquartile range $\text{IQR} = Q_3 - Q_1$, then flag values outside the fences:

$$\text{lower} = Q_1 - 1.5 \times \text{IQR}, \quad \text{upper} = Q_3 + 1.5 \times \text{IQR}$$

import numpy as np

def iqr_outliers(x: np.ndarray, factor: float = 1.5) -> np.ndarray:
    q1, q3 = np.percentile(x, [25, 75])
    iqr = q3 - q1
    lower = q1 - factor * iqr
    upper = q3 + factor * iqr
    return (x < lower) | (x > upper)  # boolean mask

data = np.array([10, 12, 11, 13, 200, 10, 9, 11])
print(data[iqr_outliers(data)])  # [200]

Robust Z-Score (Median / MAD)

Replace mean and standard deviation with their robust counterparts — the median and Median Absolute Deviation (MAD). The constant 1.4826 makes MAD a consistent estimator of $\sigma$ for normal data:

$$z_{\text{robust}} = \frac{x - \text{median}(x)}{1.4826 \times \text{MAD}(x)}, \quad \text{MAD} = \text{median}(|x - \text{median}(x)|)$$

def robust_zscore(x: np.ndarray) -> np.ndarray:
    med = np.median(x)
    mad = np.median(np.abs(x - med))
    return (x - med) / (1.4826 * mad + 1e-10)  # epsilon prevents div-by-zero

scores = robust_zscore(data)
anomalies = data[np.abs(scores) > 3.5]

Grubbs Test

A formal hypothesis test for detecting a single outlier in a univariate normal sample. The test statistic is:

$$G = \frac{\max_i |x_i - \bar{x}|}{s}$$

Reject the null (no outlier) if $G$ exceeds the critical value for sample size $n$ and significance level $\alpha$. Important limitation: it tests for exactly one outlier — repeated application (iterative Grubbs) is required to detect multiple outliers, but inflates Type I error.

Limitations of Statistical Methods

• Univariate only: All three methods work on a single feature at a time. A transaction with a normal amount AND a normal hour AND a normal merchant can still be anomalous in the combination of all three features.
• Distributional assumptions: Z-score and Grubbs require approximate normality. Skewed or multimodal distributions generate many false positives.
• Curse of dimensionality: In high-dimensional spaces, the concept of distance collapses — all points become approximately equidistant. This makes distance- and density-based outlier scores unreliable beyond ~10–20 dimensions.

3. Isolation Forest

Isolation Forest is the go-to algorithm for general-purpose anomaly detection on tabular data. It is fast, scales well to high dimensions, and requires no assumption about the data distribution.

Algorithm Intuition

The key insight: anomalies are few and different, so they are easier to isolate than normal points. The algorithm builds an ensemble of random isolation trees. At each node, a random feature is selected, and a random split value is drawn uniformly within the feature's range. This partitions the data recursively. Because anomalies occupy sparse regions of the feature space, they are separated from all other points in very few splits. Normal points, clustered in dense regions, require many splits to isolate.

The anomaly score for a point $x$ is based on its average path length $h(x)$ across all trees, normalized by the expected path length for a sample of size $n$:

$$s(x, n) = 2^{-\frac{\mathbb{E}[h(x)]}{c(n)}}, \quad c(n) = 2H(n-1) - \frac{2(n-1)}{n}$$

where $H$ is the harmonic number. Scores near 1 indicate anomalies; scores near 0.5 indicate normal points.

Key Properties

• Training complexity: $O(n \log n)$; scoring complexity: $O(\log n)$ per point.
• Does not require distance or density estimation — immune to the curse of dimensionality.
• The contamination parameter sets the decision threshold (fraction of training data expected to be anomalous).
• Typical hyperparameters: n_estimators=100 (increase to 200 for stability), max_samples="auto" (uses 256 samples per tree by default — sufficient for most datasets).

import numpy as np
from sklearn.ensemble import IsolationForest

rng = np.random.default_rng(42)
# 990 normal points + 10 anomalies
X_normal = rng.normal(loc=0, scale=1, size=(990, 5))
X_anomaly = rng.normal(loc=5, scale=0.5, size=(10, 5))
X = np.vstack([X_normal, X_anomaly])

clf = IsolationForest(
    n_estimators=200,
    contamination=0.01,   # 1% expected anomaly rate
    random_state=42
)
clf.fit(X)

scores = clf.decision_function(X)   # more negative = more anomalous
labels = clf.predict(X)             # -1 = anomaly, 1 = normal

print(f"Flagged anomalies: {(labels == -1).sum()}")   # ~10
print(f"Anomaly score range: {scores.min():.3f} to {scores.max():.3f}")

Extended Isolation Forest

Standard Isolation Forest has a known artifact: it assigns anomalous scores to the center of dense clusters due to the axis-aligned split structure. Extended Isolation Forest (EIF) uses random hyperplane splits instead of axis-aligned splits, removing this artifact — worth using when you have multiple tight clusters.

4. Local Outlier Factor (LOF)

LOF is a density-based algorithm that compares the local density of a point to the local densities of its neighbors. It detects anomalies that are in low-density regions relative to their local neighborhood — crucially, it handles datasets with clusters of varying density.

Algorithm Intuition

For each point $p$, LOF computes a local reachability density (lrd) based on the $k$-nearest neighbors. The LOF score is the ratio of the average lrd of $p$'s neighbors to $p$'s own lrd:

$$\text{LOF}_k(p) = \frac{\frac{1}{|N_k(p)|}\sum_{o \in N_k(p)} \text{lrd}_k(o)}{\text{lrd}_k(p)}$$

Interpretation: LOF ≈ 1 means the point has roughly the same density as its neighbors (normal). LOF >> 1 means the point is in a much sparser region than its neighbors (anomalous). Unlike Isolation Forest, LOF produces a relative score — a point in a sparse-but-uniform region will have LOF ≈ 1, while a point in a dense region surrounded by denser neighbors will have LOF > 1.

When to Use LOF vs. Isolation Forest

Use LOF when you have non-uniform density: datasets with multiple clusters of different densities, where a point that is "far from its neighbors" is the relevant definition of anomaly. Isolation Forest is globally biased toward sparse regions of the input space; LOF is locally adaptive.

from sklearn.neighbors import LocalOutlierFactor
import numpy as np

rng = np.random.default_rng(0)
# Two dense clusters + a few outliers
cluster_a = rng.normal(loc=[0, 0], scale=0.5, size=(200, 2))
cluster_b = rng.normal(loc=[8, 8], scale=0.5, size=(200, 2))
outliers  = rng.uniform(low=-2, high=10, size=(10, 2))
X = np.vstack([cluster_a, cluster_b, outliers])

lof = LocalOutlierFactor(
    n_neighbors=20,          # key hyperparameter: larger = smoother scores
    contamination=0.024      # ~10/410
)
labels = lof.fit_predict(X)   # -1 = anomaly, 1 = normal
scores = -lof.negative_outlier_factor_   # higher = more anomalous

print(f"Flagged: {(labels == -1).sum()}")
print(f"Top LOF scores: {sorted(scores)[-5:]}")

Limitations

• Computational cost: $O(n^2)$ for exact k-NN; approximate nearest neighbors (e.g., FAISS, Annoy) are needed for $n > 100$K.
• n_neighbors sensitivity: small values are noisy; large values miss local anomalies. Common default: 20. Tune via precision@k on labeled validation data.
• No incremental updates: adding new data requires recomputing densities — not suitable for streaming anomaly detection.

5. Autoencoders for Anomaly Detection

Autoencoders learn to compress data into a lower-dimensional latent representation and then reconstruct the original input. Trained exclusively on normal data, they reconstruct normal inputs well but fail on anomalous inputs — producing high reconstruction error. This reconstruction error becomes the anomaly score.

Architecture and Training

The encoder maps input $x \in \mathbb{R}^d$ to a latent vector $z \in \mathbb{R}^k$ (with $k \ll d$); the decoder maps $z$ back to $\hat{x} \in \mathbb{R}^d$. The model minimizes reconstruction loss on normal training data:

$$\mathcal{L} = \|x - \hat{x}\|_2^2 \quad \text{(MSE for continuous features)}$$

The bottleneck forces the model to learn the most informative structure of normal data. Anomalies, lying in regions of the input space not seen during training, map to poor reconstructions.

Threshold Selection

The reconstruction error threshold is a hyperparameter. Common approach: compute reconstruction errors on a held-out set of normal data and set the threshold at the 95th or 99th percentile. Observations exceeding this threshold in production are flagged as anomalies. The choice of percentile trades off precision vs. recall.

import numpy as np
import torch
import torch.nn as nn
from torch.utils.data import DataLoader, TensorDataset

# --- Model definition ---
class Autoencoder(nn.Module):
    def __init__(self, input_dim: int, latent_dim: int):
        super().__init__()
        self.encoder = nn.Sequential(
            nn.Linear(input_dim, 64),
            nn.ReLU(),
            nn.Linear(64, latent_dim),
        )
        self.decoder = nn.Sequential(
            nn.Linear(latent_dim, 64),
            nn.ReLU(),
            nn.Linear(64, input_dim),
        )

    def forward(self, x):
        return self.decoder(self.encoder(x))

# --- Training on normal data only ---
def train_autoencoder(X_normal: np.ndarray, latent_dim: int = 8,
                      epochs: int = 50, lr: float = 1e-3) -> Autoencoder:
    model = Autoencoder(X_normal.shape[1], latent_dim)
    opt = torch.optim.Adam(model.parameters(), lr=lr)
    criterion = nn.MSELoss()

    X_t = torch.tensor(X_normal, dtype=torch.float32)
    loader = DataLoader(TensorDataset(X_t), batch_size=256, shuffle=True)

    model.train()
    for _ in range(epochs):
        for (batch,) in loader:
            loss = criterion(model(batch), batch)
            opt.zero_grad(); loss.backward(); opt.step()
    return model

# --- Scoring and threshold ---
def reconstruction_scores(model: Autoencoder, X: np.ndarray) -> np.ndarray:
    model.eval()
    with torch.no_grad():
        X_t = torch.tensor(X, dtype=torch.float32)
        recon = model(X_t).numpy()
    return np.mean((X - recon) ** 2, axis=1)  # per-sample MSE

# Set threshold at 95th percentile of normal validation errors
# threshold = np.percentile(reconstruction_scores(model, X_val_normal), 95)
# anomalies  = reconstruction_scores(model, X_test) > threshold

Variants

• Variational Autoencoder (VAE): Learns a probabilistic latent space. Anomaly score can be the reconstruction probability (ELBO), which is better calibrated than raw MSE.
• LSTM Autoencoder: Replaces dense layers with LSTM layers for sequential data (time series, log sequences). Effective for collective anomaly detection.
• Sparse Autoencoder: Adds L1 penalty on latent activations. Forces the model to learn sparse, interpretable features.

6. Evaluation

Anomaly detection evaluation is more challenging than standard classification because of severe class imbalance and, often, the complete absence of labels.

When Labels Exist

Why AUC-PR over AUC-ROC: With 0.1% anomalies, a model that flags everything as normal achieves AUC-ROC ≈ 0.5 but precision ≈ 0. Conversely, a weak model with many false positives can achieve high AUC-ROC because the true-negative rate stays near 1.0 across thresholds. The PR curve collapses this blind spot — a random baseline scores AUC-PR equal to the prevalence (0.001), making gains visible.

When No Labels Exist

• Contamination rate calibration: Use domain knowledge to set the expected anomaly rate and verify that the flagged volume matches this prior. If fraud is known to be ~0.5%, the model should flag roughly 0.5% of transactions.
• Domain expert audit: Sample the top-k flagged cases and have analysts label them. This provides an empirical precision estimate and reveals systematic false positive patterns.
• Score stability across seeds: Run the model with multiple random seeds. Consistent flagging of the same observations provides confidence that the signal is real, not noise.
• Temporal consistency: In a streaming setting, anomaly rate should be roughly stable over time. Sudden spikes in flagging rate indicate model or data drift, not necessarily more anomalies.

Business Cost Matrix

Translate model metrics into business value using a cost matrix. For fraud detection:

$$\text{Expected Cost} = C_{FN} \times \text{FN Rate} \times P(\text{anomaly}) + C_{FP} \times \text{FP Rate} \times P(\text{normal})$$

where $C_{FN}$ is the cost of a missed fraud (dollar loss) and $C_{FP}$ is the cost of a false alarm (analyst investigation time, customer friction). Optimizing for F1 implicitly assumes $C_{FN} = C_{FP}$, which is almost never true in practice. Setting the decision threshold to minimize expected cost rather than maximize F1 gives the operationally correct result.

7. Choosing the Right Method

Decision Heuristic

• Start with Isolation Forest as the baseline for any tabular anomaly detection task.
• If you suspect non-uniform density or have meaningful local cluster structure, add LOF as a comparison.
• Use statistical methods for single-feature monitoring (e.g., data pipeline QA, feature drift alerting).
• Use an autoencoder when you have high-dimensional or sequential inputs, or when you need the model to learn complex normal patterns from raw features.
• For contextual anomalies, segment the data by context first, then apply any of the above within segments.