Anomaly Detection in Production: From Outlier Statistics to Real-Time ML Systems

Most anomaly detection tutorials stop at toy datasets and z-scores. Production systems don't. You're dealing with imbalanced classes, concept drift, alert fatigue, and the real cost of each false positive—textbook definitions don't cover that.

This article bridges the gap. We start with classic statistical tools—z-scores, IQR, Grubbs' test—then move to unsupervised ML: Isolation Forest, Local Outlier Factor, autoencoders. From there, deployment: streaming detection with windowed statistics, feature engineering for time-series, and threshold management via feedback loops.

The goal isn't just to detect anomalies. It's to build a system that survives production—complete with debugging guides, incident postmortems, and a cheat sheet for when things go wrong.

What is Anomaly Detection? Definitions and Core Concepts

Anomaly detection identifies data points, events, or observations that deviate so significantly from the rest of the dataset that they raise suspicion of being generated by a different mechanism. In production systems, these are not just statistical curiosities—they represent fraud, system failures, sensor faults, or security breaches. The core challenge is that anomalies are rare by definition, often making up less than 1% of the data, and the notion of 'normal' can shift over time (concept drift).

Three operational categories dominate practice. Supervised anomaly detection requires labeled 'normal' vs 'anomaly' data, but is rarely feasible due to class imbalance and labeling cost. Semi-supervised methods use only normal data to build a profile, flagging deviations at inference. Unsupervised methods, the most common in industry, assume anomalies are isolated and few, relying on density, distance, or reconstruction error to separate them from the majority.

A critical distinction often missed: outliers vs novelties. Outliers are present in the training data and can corrupt model fitting. Novelties appear only at inference, representing genuinely new patterns. Production pipelines must handle both, often by combining robust training (e.g., winsorizing) with online scoring. The choice of detection method depends on whether you need to explain why a point is anomalous (interpretability) or simply flag it for review (black-box scoring).

import numpy as np
from sklearn.ensemble import IsolationForest

# Simulate normal data with a few anomalies
np.random.seed(42)
normal = np.random.normal(0, 1, (1000, 2))
anomalies = np.random.uniform(low=-6, high=6, size=(20, 2))
data = np.vstack([normal, anomalies])

# Unsupervised detection
model = IsolationForest(contamination=0.02, random_state=42)
preds = model.fit_predict(data)

print(f"Detected anomalies: {(preds == -1).sum()} out of {len(data)}")
print(f"True anomalies: {len(anomalies)}")

Statistical Methods: Z-Score, IQR, Grubbs' Test, and Their Limits

Statistical outlier tests are the oldest and most interpretable tools. The Z-score method assumes data is Gaussian: a point with |Z| > 3 (or 2.5 in stricter settings) is flagged. For a sample x_i, Z_i = (x_i - μ) / σ. This works well for unimodal, symmetric distributions but fails catastrophically on skewed or multimodal data. In production, Z-scores are often computed on rolling windows (e.g., 30-day mean/std) to adapt to seasonality.

The Interquartile Range (IQR) method is non-parametric: outliers are points below Q1 - 1.5IQR or above Q3 + 1.5IQR. The 1.5 multiplier is a heuristic from Tukey, not a statistical guarantee. For large datasets, this can flag 0.7% of points as outliers under normality, but on heavy-tailed distributions it may flag 5-10%. IQR is robust to extreme values because it uses medians, but it ignores the shape of the distribution beyond quartiles.

Grubbs' test is a formal hypothesis test for a single outlier in a univariate Gaussian sample. It computes G = max|x_i - x̄| / s and compares against a critical value from the Studentized range distribution. Its limits are severe: it assumes exactly one outlier, requires normality, and is sensitive to masking (multiple outliers hiding each other). Generalized Extreme Studentized Deviate (ESD) test handles multiple outliers but still assumes normality. In practice, these tests are used for quality control (e.g., manufacturing) where Gaussian assumptions hold, but rarely for complex production data.

import numpy as np
from scipy import stats

data = np.random.normal(0, 1, 1000)
data[0] = 10  # inject outlier

# Z-score
z_scores = np.abs(stats.zscore(data))
z_outliers = np.where(z_scores > 3)[0]
print(f"Z-score outliers: {len(z_outliers)}")

# IQR
Q1, Q3 = np.percentile(data, [25, 75])
iqr = Q3 - Q1
lower = Q1 - 1.5 * iqr
upper = Q3 + 1.5 * iqr
iqr_outliers = np.where((data < lower) | (data > upper))[0]
print(f"IQR outliers: {len(iqr_outliers)}")

# Grubbs' test (two-sided)
from scipy.stats import t as t_dist
alpha = 0.05
n = len(data)
G = np.max(np.abs(data - np.mean(data))) / np.std(data, ddof=1)
t_crit = t_dist.ppf(1 - alpha/(2*n), n-2)
G_crit = (n-1) / np.sqrt(n) * np.sqrt(t_crit**2 / (n-2 + t_crit**2))
print(f"Grubbs G={G:.3f}, critical={G_crit:.3f}, outlier={G > G_crit}")

Unsupervised Machine Learning: Isolation Forest, LOF, and One-Class SVM

Isolation Forest (iForest) is the go-to algorithm for high-dimensional anomaly detection. It works by randomly partitioning the feature space with decision trees; anomalies are isolated in fewer splits because they are few and different. The anomaly score is based on the average path length: s(x) = 2^(-E(h(x))/c(n)), where E(h(x)) is the expected path length and c(n) is the average path length for a random tree. IForest handles up to hundreds of features, is linear in time and memory, and requires no distance computation. In production, train with 100 trees and subsample 256 points per tree—this is the default and works well.

Local Outlier Factor (LOF) measures local density deviation. For each point, it computes the ratio of its local density to its neighbors' densities. A point with LOF >> 1 is an outlier. LOF is excellent for detecting local anomalies (e.g., a point that is normal globally but anomalous in its neighborhood) but is O(n^2) in naive implementations. Use KD-trees or ball trees for up to 10,000 points; beyond that, approximate nearest neighbors (e.g., ANNOY) are necessary. LOF's sensitivity to the 'n_neighbors' parameter is a common pitfall—set it to 20-30 for most datasets.

One-Class SVM (OC-SVM) learns a decision boundary that encloses most of the data in a high-dimensional feature space using a kernel trick (typically RBF). It solves: minimize 0.5||w||^2 + (1/νn)Σξ_i - ρ, subject to w·Φ(x_i) ≥ ρ - ξ_i, ξ_i ≥ 0. The parameter ν (nu) controls the fraction of outliers expected (upper bound) and the fraction of support vectors (lower bound). OC-SVM works well for moderate-sized datasets (up to 100k points) but scales poorly with sample size (O(n^2) to O(n^3)). It also requires careful tuning of kernel bandwidth (gamma) and nu. In practice, iForest often outperforms OC-SVM on large, noisy datasets.

import numpy as np
from sklearn.ensemble import IsolationForest
from sklearn.neighbors import LocalOutlierFactor
from sklearn.svm import OneClassSVM

# Generate data: 2 clusters + anomalies
np.random.seed(42)
X = np.vstack([
    np.random.normal([0, 0], 0.5, (500, 2)),
    np.random.normal([3, 3], 0.5, (500, 2)),
    np.random.uniform(low=-2, high=5, size=(20, 2))
])

# Isolation Forest
iForest = IsolationForest(contamination=0.02, random_state=42)
y_pred_if = iForest.fit_predict(X)
print(f"iForest anomalies: {(y_pred_if == -1).sum()}")

# LOF
lof = LocalOutlierFactor(contamination=0.02, n_neighbors=20)
y_pred_lof = lof.fit_predict(X)
print(f"LOF anomalies: {(y_pred_lof == -1).sum()}")

# One-Class SVM
ocsvm = OneClassSVM(nu=0.02, kernel='rbf', gamma='scale')
y_pred_oc = ocsvm.fit_predict(X)
print(f"OC-SVM anomalies: {(y_pred_oc == -1).sum()}")

Deep Learning Approaches: Autoencoders, VAEs, and Time-Series Models

Autoencoders (AEs) learn a compressed representation of normal data by minimizing reconstruction error. An AE consists of an encoder f: X → Z and a decoder g: Z → X' such that the reconstruction loss L = ||X - X'||^2 is minimized. At inference, points with high reconstruction error are anomalies. The threshold is typically set as the 95th or 99th percentile of reconstruction errors on a validation set. AEs work well on high-dimensional data (images, sequences) but require large amounts of normal data (10k+ samples) and careful regularization to avoid overfitting to anomalies.

Variational Autoencoders (VAEs) extend AEs by learning a probabilistic latent space. The loss is the Evidence Lower Bound (ELBO): L = -E_{z~q}[log p(x|z)] + KL(q(z|x) || p(z)). The KL divergence term regularizes the latent space, making VAEs more robust to overfitting than vanilla AEs. For anomaly detection, the reconstruction probability (not error) is used: p(x|z) averaged over multiple samples from the latent distribution. VAEs are superior for detecting subtle anomalies that AEs might reconstruct well due to memorization.

For time-series anomaly detection, recurrent architectures (LSTM-AEs) or temporal convolutional networks (TCNs) are standard. An LSTM-AE encodes a sequence into a fixed-size vector and decodes it back. The reconstruction error at each timestep is aggregated (e.g., max or mean) to score the entire sequence. More advanced methods use attention or transformers (e.g., Anomaly Transformer) to capture long-range dependencies. In production, sliding windows of 100-500 timesteps are common, and the model is retrained weekly on the last 30 days of normal data. A critical detail: always normalize time-series per-channel with rolling statistics to handle non-stationarity.

import numpy as np
import torch
import torch.nn as nn
import torch.optim as optim

# Simple Autoencoder for anomaly detection
class Autoencoder(nn.Module):
    def __init__(self, input_dim=10, latent_dim=3):
        super().__init__()
        self.encoder = nn.Sequential(
            nn.Linear(input_dim, 8),
            nn.ReLU(),
            nn.Linear(8, latent_dim)
        )
        self.decoder = nn.Sequential(
            nn.Linear(latent_dim, 8),
            nn.ReLU(),
            nn.Linear(8, input_dim)
        )
    def forward(self, x):
        return self.decoder(self.encoder(x))

# Generate normal data
np.random.seed(42)
X_normal = np.random.normal(0, 1, (1000, 10)).astype(np.float32)
X_anomaly = np.random.uniform(-5, 5, (50, 10)).astype(np.float32)

# Train
model = Autoencoder()
optimizer = optim.Adam(model.parameters(), lr=0.001)
criterion = nn.MSELoss()
for epoch in range(50):
    model.train()
    optimizer.zero_grad()
    output = model(torch.from_numpy(X_normal))
    loss = criterion(output, torch.from_numpy(X_normal))
    loss.backward()
    optimizer.step()

# Score
model.eval()
with torch.no_grad():
    recon_normal = model(torch.from_numpy(X_normal))
    recon_anomaly = model(torch.from_numpy(X_anomaly))
    errors_normal = ((X_normal - recon_normal.numpy())**2).mean(axis=1)
    errors_anomaly = ((X_anomaly - recon_anomaly.numpy())**2).mean(axis=1)

threshold = np.percentile(errors_normal, 95)
print(f"Threshold (95th percentile): {threshold:.4f}")
print(f"Anomalies detected: {(errors_anomaly > threshold).sum()} out of {len(errors_anomaly)}")

Feature Engineering for Anomaly Detection: Aggregations, Windows, and Ratios

Feature engineering is the single highest-leverage activity in production anomaly detection. Raw metrics—CPU load, transaction amount, packet count—are rarely anomalous in isolation. The signal lives in derived features: rolling statistics, rate-of-change, and ratios that capture behavioral context. For univariate time series, a common pattern is to compute a sliding window mean μ_t and standard deviation σ_t over the last N observations, then define the anomaly score as the z-score z_t = (x_t - μ_t) / σ_t. This is simple, interpretable, and works for many monitoring use cases. The window size N must be tuned: too small and you react to noise; too large and you miss fast-evolving anomalies. A good starting point is N = 100 for 1-second metrics, or N = 20 for daily aggregates.

For multivariate data, aggregations become more powerful. In fraud detection, features like 'number of transactions in last hour per merchant category' or 'ratio of transaction amount to average amount for that user over 7 days' capture behavioral baselines. Ratios are particularly robust because they normalize for scale—a $10,000 transaction is normal for a high-net-worth user but anomalous for a student. The ratio r = amount / avg_amount_7d directly encodes this. When building these features, beware of lookahead bias: never use future data to compute a feature for the current timestamp. Use expanding or rolling windows that only include past observations. In streaming systems, this means maintaining stateful aggregators (e.g., with Flink's sliding windows or Redis sorted sets).

Windowed features also enable detection of slow-drift anomalies. Instead of comparing a point to a fixed threshold, compare it to a moving baseline. For example, a 10% increase in error rate over 5 minutes might be normal, but a 10% increase over 5 days is suspicious. Use multiple windows—short (1 min), medium (1 hour), long (24 hours)—and compute the ratio of short-term to long-term averages. This gives a 'burstiness' score. A common production pattern is to store these features in a feature store (e.g., Feast, Tecton) so they are consistent across training and serving. Missing values in windows should be handled explicitly: forward-fill for short gaps, or mark as NaN and skip scoring if too many gaps exist.

Finally, consider domain-specific transformations. For network intrusion detection, features like 'number of distinct destination IPs in last 60 seconds' or 'ratio of SYN to ACK packets' are far more informative than raw byte counts. For system health, use 'ratio of garbage collection time to total runtime' or 'p99 latency minus p50 latency' to capture tail behavior. Always validate feature importance using a simple model (e.g., Isolation Forest feature importances or SHAP values) before deploying. A feature that never fires is dead weight; a feature that fires too often is noise. The goal is a sparse set of high-signal features that make anomalies pop out like a sore thumb.

import pandas as pd
import numpy as np

def compute_rolling_features(df: pd.DataFrame, value_col: str, windows: list) -> pd.DataFrame:
    """Compute z-scores and ratio features for anomaly detection."""
    df = df.sort_values('timestamp').reset_index(drop=True)
    for w in windows:
        roll = df[value_col].rolling(window=w, min_periods=max(10, w//10))
        mu = roll.mean()
        sigma = roll.std(ddof=0)
        df[f'zscore_{w}'] = (df[value_col] - mu) / sigma.replace(0, np.nan)
        df[f'ratio_{w}_to_1h'] = df[value_col] / df[value_col].rolling(window=3600, min_periods=100).mean()
    # Burstiness: ratio of 1-min mean to 1-hour mean
    df['burstiness'] = (df[value_col].rolling(60).mean() / df[value_col].rolling(3600).mean()).fillna(1.0)
    return df

# Example usage
np.random.seed(42)
times = pd.date_range('2024-01-01', periods=10000, freq='1s')
df = pd.DataFrame({'timestamp': times, 'cpu': np.random.normal(50, 5, 10000)})
df.loc[5000:5050, 'cpu'] += 30  # inject anomaly
df = compute_rolling_features(df, 'cpu', windows=[60, 300])
print(df[['timestamp', 'cpu', 'zscore_60', 'burstiness']].tail(10))

Production Deployment: Streaming Detection, Threshold Tuning, and Alert Management

Deploying anomaly detection in production is fundamentally different from offline experimentation. The model must process data in real time, adapt to changing baselines, and not drown operators in false alarms. The first architectural decision is whether to use batch scoring (e.g., hourly Spark jobs) or streaming (e.g., Flink, Kafka Streams, or a simple Python service with Redis). For sub-second latency requirements, streaming is mandatory. A common pattern is to deploy a lightweight scoring service that consumes from a Kafka topic, computes features using stateful aggregators, and emits anomaly scores to an output topic. The model itself should be a serialized artifact (e.g., ONNX, pickle with versioning) that can be hot-swapped without downtime.

Threshold tuning is the most painful part of production anomaly detection. Static thresholds (e.g., z-score > 3) fail when data distributions shift seasonally. Instead, use dynamic thresholds based on historical percentiles. For example, set the threshold at the 99.5th percentile of the anomaly score over the last 7 days, recomputed daily. This adapts to gradual drift. For streaming, maintain an online quantile estimator (e.g., t-digest or Greenwald-Khanna) that updates with each new score. The threshold should be configurable per metric or per entity (e.g., per user, per server). A common mistake is to use a single global threshold; this guarantees either too many false positives for noisy metrics or missed anomalies for quiet ones. Use a holdout validation set to tune the threshold via precision-recall tradeoff, targeting a specific false positive rate (e.g., 1 alert per 1000 observations).

Alert management is where production systems live or die. Every alert must have a clear severity level (info, warning, critical) and a runbook. Implement deduplication: if the same anomaly persists for 5 consecutive windows, fire only one alert, not five. Use escalation policies: if an alert is not acknowledged within 10 minutes, page the on-call engineer. Integrate with incident management tools (PagerDuty, Opsgenie) and include a link to a dashboard showing the relevant metric and features. Avoid alert fatigue by setting a maximum alert rate (e.g., no more than 10 alerts per hour per service). If the system exceeds this, automatically suppress lower-severity alerts and escalate to the ML team for retuning.

Finally, build a feedback loop. Every alert that an operator dismisses as a false positive should be logged with a reason. Use this labeled data to periodically retrain or fine-tune the model. A simple approach: collect all false positives over a week, compute their feature vectors, and train a small classifier (e.g., logistic regression) to predict whether an anomaly is a false positive. This can be used as a post-filter to reduce noise. The feedback loop is the only way to improve over time; without it, the system degrades as data drifts.

import json
import time
from collections import deque
import numpy as np
from kafka import KafkaConsumer, KafkaProducer

class StreamingAnomalyDetector:
    def __init__(self, window_size=100, threshold_percentile=99.5):
        self.window = deque(maxlen=window_size)
        self.scores = deque(maxlen=10000)
        self.threshold = None
        self.threshold_percentile = threshold_percentile

    def update_threshold(self):
        if len(self.scores) > 1000:
            self.threshold = np.percentile(list(self.scores), self.threshold_percentile)

    def score(self, value):
        if len(self.window) < 10:
            self.window.append(value)
            return 0.0
        mu = np.mean(self.window)
        sigma = np.std(self.window) + 1e-9
        z = (value - mu) / sigma
        self.window.append(value)
        self.scores.append(abs(z))
        self.update_threshold()
        return abs(z)

consumer = KafkaConsumer('metrics', bootstrap_servers='localhost:9092')
producer = KafkaProducer(bootstrap_servers='localhost:9092')
detector = StreamingAnomalyDetector()

for msg in consumer:
    data = json.loads(msg.value)
    score = detector.score(data['cpu'])
    if detector.threshold and score > detector.threshold:
        alert = {'metric': 'cpu', 'value': data['cpu'], 'score': score, 'timestamp': time.time()}
        producer.send('alerts', json.dumps(alert).encode())
        print(f"ALERT: {alert}")

Monitoring and Debugging: Drift Detection, False Positive Analysis, and Incident Response

Once an anomaly detection system is live, the work shifts to monitoring the monitor. The model itself can drift: data distributions change, new patterns emerge, and the threshold that worked last month may now generate 100 false positives per hour. The first line of defense is drift detection on the input features. Use a two-sample statistical test (e.g., Kolmogorov-Smirnov test or Population Stability Index) to compare the current day's feature distribution to a reference distribution (e.g., the last 30 days). If the PSI exceeds 0.2, trigger an alert to the ML team. For streaming, implement a lightweight drift detector using a sliding window of the last 10,000 observations and compare to a baseline window. This catches gradual shifts before they cause alert storms.

False positive analysis is a continuous process. Every alert that an operator dismisses should be logged with a reason (e.g., 'scheduled maintenance', 'known traffic spike', 'data quality issue'). Aggregate these reasons weekly to identify systemic issues. For example, if 30% of false positives are due to scheduled deployments, add a feature flag to suppress alerts during maintenance windows. If 20% are due to missing data (e.g., NaN values), improve the data pipeline. Use a confusion matrix over the last 7 days to compute precision and recall. If precision drops below 0.5, the threshold is too loose; if recall drops below 0.8, it's too tight. Automate threshold tuning by running a grid search over the last 7 days of labeled data and selecting the threshold that maximizes F1 score.

Incident response for anomaly detection systems follows a standard playbook. When a high-severity alert fires, the on-call engineer should first verify the data: is the metric real or a sensor glitch? Check the raw data source and the feature computation pipeline. Next, check if the anomaly is part of a known pattern (e.g., a weekly batch job). If not, escalate to the appropriate team (e.g., infrastructure, security, fraud). After the incident, conduct a post-mortem: what caused the anomaly? Was it a true positive or false positive? How can the model be improved to catch it earlier or avoid it? Document the findings in a runbook. Over time, the runbook becomes a knowledge base that reduces mean time to resolution (MTTR).

Finally, monitor the monitor's performance metrics: alert rate, false positive rate, detection latency, and model staleness. Set up a dashboard that shows these metrics over time. If the alert rate suddenly drops to zero, the model may have broken (e.g., a feature is always NaN). If latency spikes, the streaming pipeline may be backlogged. Use these metrics to trigger alerts on the monitoring system itself—meta-alerts. Without this, you risk a silent failure where the anomaly detector stops detecting anything, and you only find out after a major incident.

import numpy as np
from scipy.stats import ks_2samp

def compute_psi(expected: np.ndarray, actual: np.ndarray, bins: int = 10) -> float:
    """Population Stability Index (PSI) for drift detection."""
    expected_hist, edges = np.histogram(expected, bins=bins, range=(0, 100))
    actual_hist, _ = np.histogram(actual, bins=edges)
    expected_pct = expected_hist / len(expected) + 1e-9
    actual_pct = actual_hist / len(actual) + 1e-9
    psi = np.sum((actual_pct - expected_pct) * np.log(actual_pct / expected_pct))
    return psi

def detect_drift(reference: np.ndarray, current: np.ndarray, threshold: float = 0.2) -> bool:
    """Return True if drift detected via PSI or KS test."""
    psi = compute_psi(reference, current)
    ks_stat, ks_p = ks_2samp(reference, current)
    if psi > threshold or ks_p < 0.05:
        return True
    return False

# Example: compare last 7 days to previous 30 days
np.random.seed(42)
ref = np.random.normal(50, 5, 30000)  # 30 days of 1-second data
curr = np.random.normal(55, 5, 7000)  # drift: mean shifted to 55
print(f"Drift detected: {detect_drift(ref, curr)}")  # True
print(f"PSI: {compute_psi(ref, curr):.4f}")

Case Studies: Fraud Detection, Intrusion Detection, and System Health Monitoring

Fraud detection is the canonical application of anomaly detection. In credit card transactions, the goal is to flag transactions that deviate from a user's typical spending pattern. A production system at a major bank might process 10,000 transactions per second, scoring each with an ensemble of models: a lightweight rule-based filter (e.g., transaction amount > 3x the user's 30-day average) followed by a gradient-boosted tree (e.g., XGBoost) trained on features like 'time since last transaction', 'merchant category frequency', and 'geographic distance from last transaction'. The key insight is that fraudsters often test small amounts first, so features like 'number of transactions in last hour' are highly predictive. The system must balance false positives (annoying customers) against false negatives (financial loss). A typical target is a 1% false positive rate with 80% recall. The model is retrained daily on the last 30 days of labeled data (confirmed fraud + chargebacks).

Intrusion detection in network security uses anomaly detection to identify malicious traffic. A classic approach is to model 'normal' network behavior using statistical profiles of packet headers and flow records. For example, the KDD Cup 1999 dataset (still used as a benchmark) includes features like 'duration', 'protocol type', 'service', 'flag', and 'src_bytes'. Modern systems use deep packet inspection and flow-level features aggregated over sliding windows. A production IDS at a large cloud provider might process 1 million flows per second, using a streaming Isolation Forest model that updates its isolation trees every hour. The model flags flows that are isolated in shallow trees (short path length). The challenge is the extreme class imbalance: malicious traffic is often less than 0.001% of all traffic. To handle this, the system uses a two-stage approach: a cheap pre-filter (e.g., rule-based) that passes 1% of traffic to the expensive model, reducing the scoring load by 99%.

System health monitoring uses anomaly detection to catch infrastructure failures before they cause outages. For example, a cloud platform monitors hundreds of metrics per server: CPU, memory, disk I/O, network latency, error rates. The goal is to detect anomalies that precede a service degradation. A common pattern is to use a multivariate model like a Variational Autoencoder (VAE) trained on normal operation data. The reconstruction error serves as the anomaly score. In production, this model runs on a 1-minute sliding window, scoring each server's metric vector. When the reconstruction error exceeds a threshold (e.g., 99.9th percentile of the last 24 hours), an alert fires. A real-world deployment at a large SaaS company reduced mean time to detection (MTTD) from 15 minutes to 2 minutes, and reduced false positives by 40% compared to static thresholds. The key was feature engineering: adding ratios like 'disk_io / cpu_usage' to capture correlated anomalies.

Across all three domains, the common success factors are: (1) rich feature engineering with domain knowledge, (2) adaptive thresholds that handle seasonality and drift, (3) a feedback loop to continuously improve the model, and (4) a robust alerting system that avoids fatigue. The failures all share the same root cause: treating anomaly detection as a one-time model-building exercise rather than an ongoing operational discipline. The most successful teams treat their anomaly detector as a living system that requires daily attention, much like a production database.

import pandas as pd
import numpy as np
from sklearn.ensemble import IsolationForest
from sklearn.metrics import precision_recall_curve

# Simulate fraud detection data
np.random.seed(42)
n = 100000
normal = pd.DataFrame({
    'amount': np.random.exponential(50, n),
    'time_since_last': np.random.exponential(3600, n),
    'merchant_freq': np.random.poisson(5, n),
    'is_fraud': 0
})
# Inject fraud: high amount, short time since last, rare merchant
fraud = pd.DataFrame({
    'amount': np.random.exponential(500, 100),
    'time_since_last': np.random.exponential(60, 100),
    'merchant_freq': np.random.poisson(0.5, 100),
    'is_fraud': 1
})
df = pd.concat([normal, fraud], ignore_index=True).sample(frac=1)

# Train Isolation Forest
model = IsolationForest(contamination=0.001, random_state=42)
model.fit(df[['amount', 'time_since_last', 'merchant_freq']])
df['score'] = -model.decision_function(df[['amount', 'time_since_last', 'merchant_freq']])

# Evaluate
precision, recall, thresholds = precision_recall_curve(df['is_fraud'], df['score'])
# Find threshold for 80% recall
target_recall = 0.8
idx = np.argmin(np.abs(recall - target_recall))
print(f"Threshold for {target_recall*100:.0f}% recall: {thresholds[idx]:.4f}")
print(f"Precision at that threshold: {precision[idx]:.4f}")