Anomaly Detection

Anomaly Detection

Anomaly detection identifies data points, events, or observations that deviate significantly from the expected pattern, potentially indicating issues or important events.

Types of Anomalies

• Point Anomalies: Individual data points that are far from normal
• Contextual Anomalies: Anomalies that are context-dependent
• Collective Anomalies: A collection of data points that are anomalous together

Anomaly Detection Architecture

Statistical Methods

Z-Score Method

import numpy as np
from typing import Dict, List, Tuple

class ZScoreAnomalyDetector:
    def __init__(self, threshold=3.0, window_size=None):
        self.threshold = threshold
        self.window_size = window_size
        self.mean = None
        self.std = None
    
    def fit(self, data: np.ndarray):
        """Fit the detector to normal data"""
        self.mean = np.mean(data)
        self.std = np.std(data)
        
        if self.std == 0:
            self.std = 1e-10  # Avoid division by zero
    
    def detect(self, data: np.ndarray) -> Dict:
        """Detect anomalies in data"""
        if self.mean is None or self.std is None:
            raise ValueError("Detector not fitted. Call fit() first.")
        
        z_scores = np.abs((data - self.mean) / self.std)
        anomalies = z_scores > self.threshold
        
        return {
            "anomalies": anomalies,
            "z_scores": z_scores,
            "anomaly_indices": np.where(anomalies)[0],
            "anomaly_count": np.sum(anomalies),
            "anomaly_rate": np.mean(anomalies)
        }
    
    def detect_streaming(self, value: float) -> bool:
        """Detect anomaly in streaming data"""
        if self.mean is None or self.std is None:
            return False
        
        z_score = abs((value - self.mean) / self.std)
        return z_score > self.threshold

class RollingZScoreDetector:
    def __init__(self, window_size: int, threshold: float = 3.0):
        self.window_size = window_size
        self.threshold = threshold
        self.data_buffer = []
    
    def detect(self, value: float) -> Dict:
        """Detect anomaly using rolling statistics"""
        self.data_buffer.append(value)
        
        if len(self.data_buffer) > self.window_size:
            self.data_buffer.pop(0)
        
        if len(self.data_buffer) < 10:  # Need minimum data
            return {"is_anomaly": False, "z_score": 0.0}
        
        mean = np.mean(self.data_buffer)
        std = np.std(self.data_buffer)
        
        if std == 0:
            return {"is_anomaly": False, "z_score": 0.0}
        
        z_score = abs((value - mean) / std)
        is_anomaly = z_score > self.threshold
        
        return {
            "is_anomaly": is_anomaly,
            "z_score": z_score,
            "mean": mean,
            "std": std
        }

IQR Method

import numpy as np
from typing import Dict

class IQRAnomalyDetector:
    def __init__(self, multiplier=1.5):
        self.multiplier = multiplier
        self.q1 = None
        self.q3 = None
        self.iqr = None
    
    def fit(self, data: np.ndarray):
        """Fit the detector to normal data"""
        self.q1 = np.percentile(data, 25)
        self.q3 = np.percentile(data, 75)
        self.iqr = self.q3 - self.q1
        
        if self.iqr == 0:
            self.iqr = 1e-10  # Avoid division by zero
    
    def detect(self, data: np.ndarray) -> Dict:
        """Detect anomalies using IQR"""
        if self.q1 is None or self.q3 is None:
            raise ValueError("Detector not fitted. Call fit() first.")
        
        lower_bound = self.q1 - self.multiplier * self.iqr
        upper_bound = self.q3 + self.multiplier * self.iqr
        
        anomalies = (data < lower_bound) | (data > upper_bound)
        
        return {
            "anomalies": anomalies,
            "lower_bound": lower_bound,
            "upper_bound": upper_bound,
            "anomaly_indices": np.where(anomalies)[0],
            "anomaly_count": np.sum(anomalies),
            "anomaly_rate": np.mean(anomalies)
        }
    
    def detect_single(self, value: float) -> bool:
        """Detect if a single value is anomalous"""
        if self.q1 is None or self.q3 is None:
            return False
        
        lower_bound = self.q1 - self.multiplier * self.iqr
        upper_bound = self.q3 + self.multiplier * self.iqr
        
        return value < lower_bound or value > upper_bound

ML-Based Methods

Isolation Forest

import numpy as np
from typing import Dict, List
import random

class IsolationTree:
    def __init__(self, max_depth=10):
        self.max_depth = max_depth
        self.tree = None
    
    def fit(self, data: np.ndarray, depth=0):
        """Build isolation tree"""
        if depth >= self.max_depth or len(data) <= 1:
            return {"type": "leaf", "size": len(data)}
        
        # Random feature and split value
        n_features = data.shape[1]
        feature_idx = random.randint(0, n_features - 1)
        
        feature_min = np.min(data[:, feature_idx])
        feature_max = np.max(data[:, feature_idx])
        
        if feature_min == feature_max:
            return {"type": "leaf", "size": len(data)}
        
        split_value = random.uniform(feature_min, feature_max)
        
        left_mask = data[:, feature_idx] < split_value
        right_mask = ~left_mask
        
        left_data = data[left_mask]
        right_data = data[right_mask]
        
        return {
            "type": "internal",
            "feature_idx": feature_idx,
            "split_value": split_value,
            "left": self.fit(left_data, depth + 1),
            "right": self.fit(right_data, depth + 1)
        }
    
    def path_length(self, point: np.ndarray, node=None, depth=0) -> float:
        """Calculate path length for a point"""
        if node is None:
            node = self.tree
        
        if node["type"] == "leaf":
            return depth
        
        if point[node["feature_idx"]] < node["split_value"]:
            return self.path_length(point, node["left"], depth + 1)
        else:
            return self.path_length(point, node["right"], depth + 1)

class IsolationForest:
    def __init__(self, n_trees=100, max_depth=10, contamination=0.1):
        self.n_trees = n_trees
        self.max_depth = max_depth
        self.contamination = contamination
        self.trees = []
        self.threshold = None
    
    def fit(self, data: np.ndarray):
        """Fit isolation forest"""
        self.trees = []
        
        for _ in range(self.n_trees):
            # Sample subset
            sample_size = min(256, len(data))
            sample_indices = random.sample(range(len(data)), sample_size)
            sample_data = data[sample_indices]
            
            # Build tree
            tree = IsolationTree(self.max_depth)
            tree.tree = tree.fit(sample_data)
            self.trees.append(tree)
        
        # Calculate threshold
        scores = self._calculate_scores(data)
        self.threshold = np.percentile(scores, 100 * self.contamination)
    
    def _calculate_scores(self, data: np.ndarray) -> np.ndarray:
        """Calculate anomaly scores"""
        scores = []
        
        for point in data:
            avg_path_length = np.mean([tree.path_length(point) for tree in self.trees])
            scores.append(self._anomaly_score(avg_path_length))
        
        return np.array(scores)
    
    def _anomaly_score(self, path_length: float) -> float:
        """Calculate anomaly score from path length"""
        c = self._average_path_length(256)  # Average path length for 256 samples
        return 2 ** (-path_length / c)
    
    def _average_path_length(self, n: int) -> float:
        """Calculate average path length for unsuccessful search in BST"""
        if n <= 1:
            return 0
        return 2.0 * (np.log(n - 1) + 0.5772156649) - (2.0 * (n - 1) / n)
    
    def detect(self, data: np.ndarray) -> Dict:
        """Detect anomalies"""
        scores = self._calculate_scores(data)
        anomalies = scores > self.threshold
        
        return {
            "anomalies": anomalies,
            "scores": scores,
            "threshold": self.threshold,
            "anomaly_indices": np.where(anomalies)[0],
            "anomaly_count": np.sum(anomalies),
            "anomaly_rate": np.mean(anomalies)
        }

Autoencoder-based Detection

import numpy as np
from typing import Dict, Tuple

class AutoencoderDetector:
    def __init__(self, input_dim: int, encoding_dim: int = 32):
        self.input_dim = input_dim
        self.encoding_dim = encoding_dim
        self.encoder_weights = None
        self.decoder_weights = None
        self.threshold = None
    
    def _sigmoid(self, x):
        return 1 / (1 + np.exp(-np.clip(x, -500, 500)))
    
    def _sigmoid_derivative(self, x):
        return x * (1 - x)
    
    def fit(self, data: np.ndarray, epochs=100, learning_rate=0.01):
        """Train autoencoder"""
        n_samples = data.shape[0]
        
        # Initialize weights
        self.encoder_weights = np.random.randn(self.input_dim, self.encoding_dim) * 0.1
        self.decoder_weights = np.random.randn(self.encoding_dim, self.input_dim) * 0.1
        
        for epoch in range(epochs):
            # Forward pass
            encoded = self._sigmoid(np.dot(data, self.encoder_weights))
            decoded = self._sigmoid(np.dot(encoded, self.decoder_weights))
            
            # Calculate loss
            loss = np.mean((data - decoded) ** 2)
            
            # Backward pass
            output_error = data - decoded
            output_delta = output_error * self._sigmoid_derivative(decoded)
            
            hidden_error = np.dot(output_delta, self.decoder_weights.T)
            hidden_delta = hidden_error * self._sigmoid_derivative(encoded)
            
            # Update weights
            self.decoder_weights += learning_rate * np.dot(encoded.T, output_delta) / n_samples
            self.encoder_weights += learning_rate * np.dot(data.T, hidden_delta) / n_samples
            
            if epoch % 10 == 0:
                print(f"Epoch {epoch}, Loss: {loss:.6f}")
        
        # Calculate threshold
        reconstructions = self._reconstruct(data)
        mse = np.mean((data - reconstructions) ** 2, axis=1)
        self.threshold = np.percentile(mse, 95)  # 95th percentile
    
    def _reconstruct(self, data: np.ndarray) -> np.ndarray:
        """Reconstruct data"""
        encoded = self._sigmoid(np.dot(data, self.encoder_weights))
        decoded = self._sigmoid(np.dot(encoded, self.decoder_weights))
        return decoded
    
    def detect(self, data: np.ndarray) -> Dict:
        """Detect anomalies using reconstruction error"""
        reconstructions = self._reconstruct(data)
        mse = np.mean((data - reconstructions) ** 2, axis=1)
        anomalies = mse > self.threshold
        
        return {
            "anomalies": anomalies,
            "reconstruction_errors": mse,
            "threshold": self.threshold,
            "anomaly_indices": np.where(anomalies)[0],
            "anomaly_count": np.sum(anomalies),
            "anomaly_rate": np.mean(anomalies)
        }

Mathematical Foundation

Z-Score Formula

Where:

• ( X ) is the data point
• ( \mu ) is the mean
• ( \sigma ) is the standard deviation

IQR Formula

Anomaly boundaries:

Isolation Forest Anomaly Score

Where:

• ( E(h(x)) ) is the average path length
• ( c(n) = 2H(n-1) - \frac{2(n-1)}{n} ) is the average path length of unsuccessful search in BST
• ( H(i) ) is the harmonic number

Ensemble Methods

import numpy as np
from typing import Dict, List

class AnomalyEnsemble:
    def __init__(self):
        self.detectors = []
        self.weights = []
    
    def add_detector(self, detector, weight=1.0):
        """Add detector to ensemble"""
        self.detectors.append(detector)
        self.weights.append(weight)
    
    def detect(self, data: np.ndarray) -> Dict:
        """Detect anomalies using ensemble"""
        all_scores = []
        
        for detector, weight in zip(self.detectors, self.weights):
            result = detector.detect(data)
            scores = result.get("scores", result.get("reconstruction_errors", None))
            
            if scores is not None:
                # Normalize scores
                normalized_scores = (scores - np.min(scores)) / (np.max(scores) - np.min(scores) + 1e-10)
                all_scores.append(normalized_scores * weight)
        
        if not all_scores:
            return {"anomalies": np.array([]), "scores": np.array([])}
        
        # Combine scores
        combined_scores = np.mean(all_scores, axis=0)
        
        # Determine threshold
        threshold = np.percentile(combined_scores, 95)
        anomalies = combined_scores > threshold
        
        return {
            "anomalies": anomalies,
            "scores": combined_scores,
            "threshold": threshold,
            "anomaly_indices": np.where(anomalies)[0],
            "anomaly_count": np.sum(anomalies),
            "anomaly_rate": np.mean(anomalies)
        }

Best Practices

1. Feature Engineering

class AnomalyFeatureEngine:
    def __init__(self):
        self.feature_extractors = []
    
    def add_extractor(self, extractor):
        """Add feature extractor"""
        self.feature_extractors.append(extractor)
    
    def extract_features(self, data: np.ndarray) -> np.ndarray:
        """Extract features for anomaly detection"""
        features = []
        
        for extractor in self.feature_extractors:
            feature = extractor(data)
            features.append(feature)
        
        return np.column_stack(features)

# Example feature extractors
def statistical_features(data):
    """Extract statistical features"""
    return [
        np.mean(data),
        np.std(data),
        np.min(data),
        np.max(data),
        np.median(data)
    ]

def trend_features(data):
    """Extract trend features"""
    x = np.arange(len(data))
    slope = np.polyfit(x, data, 1)[0]
    return [slope]

2. Threshold Optimization

class ThresholdOptimizer:
    def __init__(self):
        self.thresholds = {}
    
    def optimize_threshold(self, data: np.ndarray, labels: np.ndarray = None):
        """Optimize anomaly threshold"""
        if labels is not None:
            # Supervised optimization
            return self._supervised_optimization(data, labels)
        else:
            # Unsupervised optimization
            return self._unsupervised_optimization(data)
    
    def _supervised_optimization(self, data: np.ndarray, labels: np.ndarray):
        """Optimize threshold with labeled data"""
        best_threshold = 0
        best_f1 = 0
        
        for threshold in np.percentile(data, np.arange(80, 99, 1)):
            predictions = data > threshold
            
            # Calculate F1 score
            tp = np.sum(predictions & labels)
            fp = np.sum(predictions & ~labels)
            fn = np.sum(~predictions & labels)
            
            precision = tp / (tp + fp) if (tp + fp) > 0 else 0
            recall = tp / (tp + fn) if (tp + fn) > 0 else 0
            f1 = 2 * precision * recall / (precision + recall) if (precision + recall) > 0 else 0
            
            if f1 > best_f1:
                best_f1 = f1
                best_threshold = threshold
        
        return best_threshold
    
    def _unsupervised_optimization(self, data: np.ndarray):
        """Optimize threshold without labels"""
        # Use elbow method or other heuristics
        thresholds = np.percentile(data, np.arange(80, 99, 1))
        scores = []
        
        for threshold in thresholds:
            anomaly_rate = np.mean(data > threshold)
            scores.append(anomaly_rate)
        
        # Find elbow point
        diffs = np.diff(scores)
        elbow_idx = np.argmax(np.abs(diffs))
        
        return thresholds[elbow_idx]

Summary

Anomaly detection is essential for identifying unusual patterns in ML systems. By implementing statistical and ML-based methods, organizations can detect issues early and maintain system reliability through proactive monitoring and response.