Z-Score Formula: Standardization, Anomaly Detection, and Statistical Process Control in Production Systems

import math
from dataclasses import dataclass
from typing import List, Tuple, Optional


@dataclass
class ZScoreResult:
    """Result of z-score calculation for a single value."""
    value: float
    z_score: float
    is_outlier: bool
    percentile: float
    interpretation: str


class ZScoreCalculator:
    """Production-grade z-score computation with distribution-aware thresholds."""

    def calculate_mean(self, data: List[float]) -> float:
        """Calculate arithmetic mean."""
        if not data:
            raise ValueError("Cannot calculate mean of empty dataset")
        return sum(data) / len(data)

    def calculate_stddev(self, data: List[float], ddof: int = 1) -> float:
        """
        Calculate standard deviation.
        ddof=0 for population, ddof=1 for sample (Bessel's correction).
        """
        if len(data) < 2:
            raise ValueError("Need at least 2 data points for sample stddev")
        mean = self.calculate_mean(data)
        variance = sum((x - mean) ** 2 for x in data) / (len(data) - ddof)
        return math.sqrt(variance)

    def calculate_zscore(self, x: float, mean: float, stddev: float) -> float:
        """
        Calculate z-score for a single value.
        z = (x - mean) / stddev
        """
        if stddev == 0:
            return 0.0  # all values are identical
        return (x - mean) / stddev

    def calculate_zscores(self, data: List[float]) -> List[float]:
        """Calculate z-scores for an entire dataset."""
        mean = self.calculate_mean(data)
        stddev = self.calculate_stddev(data)
        return [self.calculate_zscore(x, mean, stddev) for x in data]

    def detect_outliers(self, data: List[float], threshold: float = 3.0) -> List[ZScoreResult]:
        """
        Detect outliers using z-score threshold.
        Default threshold of 3.0 catches 99.73% of normal data.
        """
        mean = self.calculate_mean(data)
        stddev = self.calculate_stddev(data)

        results = []
        for x in data:
            z = self.calculate_zscore(x, mean, stddev)
            is_outlier = abs(z) > threshold

            results.append(ZScoreResult(
                value=x,
                z_score=round(z, 4),
                is_outlier=is_outlier,
                percentile=round(self._z_to_percentile(z), 4),
                interpretation=self._interpret_zscore(z),
            ))

        return results

    def _z_to_percentile(self, z: float) -> float:
        """
        Convert z-score to percentile using approximation of the normal CDF.
        Uses Abramowitz and Stegun approximation (error < 7.5e-8).
        """
        if z < -8:
            return 0.0
        if z > 8:
            return 100.0

        # Approximation of the standard normal CDF
        sign = 1 if z >= 0 else -1
        z = abs(z)

        t = 1.0 / (1.0 + 0.2316419 * z)
        d = 0.3989422804014327  # 1/sqrt(2*pi)
        p = d * math.exp(-z * z / 2.0) * t * (
            0.319381530 + t * (-0.356563782 + t * (1.781477937 + t * (-1.821255978 + t * 1.330274429)))
        )

        percentile = 1.0 - p
        if sign < 0:
            percentile = 1.0 - percentile

        return percentile * 100.0

    def _interpret_zscore(self, z: float) -> str:
        """Interpret the magnitude of a z-score."""
        abs_z = abs(z)
        if abs_z < 1:
            return "Within 1 standard deviation — common (68% of data)"
        elif abs_z < 2:
            return "Within 2 standard deviations — typical (95% of data)"
        elif abs_z < 3:
            return "Between 2 and 3 standard deviations — unusual (5% of data)"
        else:
            return f"Beyond 3 standard deviations — {abs_z:.2f} sigma outlier (rare, <0.3%)"

import math
import time
from dataclasses import dataclass
from typing import List, Optional, Dict
from collections import deque


@dataclass
class AnomalyEvent:
    """An anomaly detected by the z-score detector."""
    timestamp: float
    value: float
    z_score: float
    threshold: float
    severity: str  # 'warning' or 'critical'
    baseline_mean: float
    baseline_stddev: float


class ZScoreAnomalyDetector:
    """Production z-score anomaly detector with rolling baselines and segmentation."""

    def __init__(self, window_size: int = 1440, warning_threshold: float = 2.5, critical_threshold: float = 3.5):
        """
        window_size: number of data points in rolling baseline (default: 1440 = 24 hours at 1/min)
        warning_threshold: z-score for warning alerts
        critical_threshold: z-score for critical alerts
        """
        self.window_size = window_size
        self.warning_threshold = warning_threshold
        self.critical_threshold = critical_threshold
        self.data_window: deque = deque(maxlen=window_size)
        self.segmented_windows: Dict[int, deque] = {}
        self.segment_size = 60  # 60 data points per segment (1 hour at 1/min)

    def add_value(self, value: float, timestamp: Optional[float] = None) -> Optional[AnomalyEvent]:
        """
        Add a new value and check for anomaly.
        Returns AnomalyEvent if anomaly detected, None otherwise.
        """
        if timestamp is None:
            timestamp = time.time()

        if len(self.data_window) < self.window_size:
            self.data_window.append(value)
            return None  # not enough data for baseline

        # Compute baseline from current window
        mean = self._mean(self.data_window)
        stddev = self._stddev(self.data_window)

        if stddev == 0:
            self.data_window.append(value)
            return None  # no variance in baseline

        z = (value - mean) / stddev

        # Add to window after computing z-score (don't let current anomaly inflate baseline)
        self.data_window.append(value)

        if abs(z) > self.critical_threshold:
            return AnomalyEvent(
                timestamp=timestamp,
                value=value,
                z_score=round(z, 4),
                threshold=self.critical_threshold,
                severity='critical',
                baseline_mean=round(mean, 4),
                baseline_stddev=round(stddev, 4),
            )
        elif abs(z) > self.warning_threshold:
            return AnomalyEvent(
                timestamp=timestamp,
                value=value,
                z_score=round(z, 4),
                threshold=self.warning_threshold,
                severity='warning',
                baseline_mean=round(mean, 4),
                baseline_stddev=round(stddev, 4),
            )

        return None

    def add_value_segmented(self, value: float, hour: int, timestamp: Optional[float] = None) -> Optional[AnomalyEvent]:
        """
        Add value with hour-of-day segmentation.
        Each hour has its own baseline, capturing diurnal patterns.
        """
        if timestamp is None:
            timestamp = time.time()

        if hour not in self.segmented_windows:
            self.segmented_windows[hour] = deque(maxlen=self.segment_size * 7)  # 7 days of this hour

        segment = self.segmented_windows[hour]

        if len(segment) < self.segment_size:
            segment.append(value)
            return None

        mean = self._mean(segment)
        stddev = self._stddev(segment)

        if stddev == 0:
            segment.append(value)
            return None

        z = (value - mean) / stddev
        segment.append(value)

        if abs(z) > self.critical_threshold:
            return AnomalyEvent(
                timestamp=timestamp,
                value=value,
                z_score=round(z, 4),
                threshold=self.critical_threshold,
                severity='critical',
                baseline_mean=round(mean, 4),
                baseline_stddev=round(stddev, 4),
            )
        elif abs(z) > self.warning_threshold:
            return AnomalyEvent(
                timestamp=timestamp,
                value=value,
                z_score=round(z, 4),
                threshold=self.warning_threshold,
                severity='warning',
                baseline_mean=round(mean, 4),
                baseline_stddev=round(stddev, 4),
            )

        return None

    def detect_with_log_transform(self, value: float, timestamp: Optional[float] = None) -> Optional[AnomalyEvent]:
        """
        Detect anomalies using z-scores on log-transformed data.
        Use for right-skewed metrics (latency, revenue, request sizes).
        """
        if value <= 0:
            return None  # log undefined for non-positive values

        log_value = math.log(value)

        # Store log-transformed values in a separate window
        if not hasattr(self, '_log_window'):
            self._log_window: deque = deque(maxlen=self.window_size)

        if len(self._log_window) < self.window_size:
            self._log_window.append(log_value)
            return None

        mean = self._mean(self._log_window)
        stddev = self._stddev(self._log_window)

        if stddev == 0:
            self._log_window.append(log_value)
            return None

        z = (log_value - mean) / stddev
        self._log_window.append(log_value)

        if abs(z) > self.critical_threshold:
            return AnomalyEvent(
                timestamp=timestamp or time.time(),
                value=value,
                z_score=round(z, 4),
                threshold=self.critical_threshold,
                severity='critical',
                baseline_mean=round(math.exp(mean), 4),
                baseline_stddev=round(math.exp(mean + stddev) - math.exp(mean), 4),
            )
        elif abs(z) > self.warning_threshold:
            return AnomalyEvent(
                timestamp=timestamp or time.time(),
                value=value,
                z_score=round(z, 4),
                threshold=self.warning_threshold,
                severity='warning',
                baseline_mean=round(math.exp(mean), 4),
                baseline_stddev=round(math.exp(mean + stddev) - math.exp(mean), 4),
            )

        return None

    def _mean(self, window: deque) -> float:
        return sum(window) / len(window)

    def _stddev(self, window: deque) -> float:
        mean = self._mean(window)
        return math.sqrt(sum((x - mean) ** 2 for x in window) / (len(window) - 1))

import math
from dataclasses import dataclass
from typing import List, Tuple, Optional
from enum import Enum


class ScalingMethod(Enum):
    ZSCORE = 'zscore'
    MINMAX = 'minmax'
    ROBUST = 'robust'
    LOG_ZSCORE = 'log_zscore'


@dataclass
class ScalingParams:
    """Parameters needed to apply the same scaling to new data."""
    method: ScalingMethod
    mean: Optional[float] = None
    stddev: Optional[float] = None
    min_val: Optional[float] = None
    max_val: Optional[float] = None
    median: Optional[float] = None
    iqr: Optional[float] = None


class FeatureNormalizer:
    """Production feature normalization with automatic method selection."""

    def zscore_normalize(self, data: List[float]) -> Tuple[List[float], ScalingParams]:
        """
        Standard z-score normalization: (x - mean) / stddev.
        Output has mean=0, stddev=1.
        """
        mean = sum(data) / len(data)
        stddev = math.sqrt(sum((x - mean) ** 2 for x in data) / (len(data) - 1))

        if stddev == 0:
            return [0.0] * len(data), ScalingParams(method=ScalingMethod.ZSCORE, mean=mean, stddev=0)

        normalized = [(x - mean) / stddev for x in data]
        return normalized, ScalingParams(method=ScalingMethod.ZSCORE, mean=mean, stddev=stddev)

    def minmax_normalize(self, data: List[float]) -> Tuple[List[float], ScalingParams]:
        """
        Min-max normalization: (x - min) / (max - min).
        Output is in range [0, 1].
        """
        min_val = min(data)
        max_val = max(data)
        range_val = max_val - min_val

        if range_val == 0:
            return [0.5] * len(data), ScalingParams(method=ScalingMethod.MINMAX, min_val=min_val, max_val=max_val)

        normalized = [(x - min_val) / range_val for x in data]
        return normalized, ScalingParams(method=ScalingMethod.MINMAX, min_val=min_val, max_val=max_val)

    def robust_normalize(self, data: List[float]) -> Tuple[List[float], ScalingParams]:
        """
        Robust scaling: (x - median) / IQR.
        Uses median and interquartile range instead of mean and stddev.
        Resistant to outliers.
        """
        sorted_data = sorted(data)
        n = len(sorted_data)
        median = sorted_data[n // 2] if n % 2 == 1 else (sorted_data[n // 2 - 1] + sorted_data[n // 2]) / 2

        q1_idx = n // 4
        q3_idx = 3 * n // 4
        q1 = sorted_data[q1_idx]
        q3 = sorted_data[q3_idx]
        iqr = q3 - q1

        if iqr == 0:
            return [0.0] * len(data), ScalingParams(method=ScalingMethod.ROBUST, median=median, iqr=0)

        normalized = [(x - median) / iqr for x in data]
        return normalized, ScalingParams(method=ScalingMethod.ROBUST, median=median, iqr=iqr)

    def log_zscore_normalize(self, data: List[float]) -> Tuple[List[float], ScalingParams]:
        """
        Log-transform followed by z-score normalization.
        Use for right-skewed features (latency, revenue, counts).
        """
        log_data = [math.log(x) if x > 0 else math.log(1e-10) for x in data]
        mean = sum(log_data) / len(log_data)
        stddev = math.sqrt(sum((x - mean) ** 2 for x in log_data) / (len(log_data) - 1))

        if stddev == 0:
            return [0.0] * len(data), ScalingParams(method=ScalingMethod.LOG_ZSCORE, mean=mean, stddev=0)

        normalized = [(x - mean) / stddev for x in log_data]
        return normalized, ScalingParams(method=ScalingMethod.LOG_ZSCORE, mean=mean, stddev=stddev)

    def auto_select_method(self, data: List[float]) -> ScalingMethod:
        """
        Automatically select the best normalization method based on data characteristics.
        """
        n = len(data)
        sorted_data = sorted(data)
        mean = sum(data) / n
        median = sorted_data[n // 2]
        stddev = math.sqrt(sum((x - mean) ** 2 for x in data) / (n - 1))

        # Check for zeros or negatives
        if any(x <= 0 for x in data):
            # Cannot use log-transform
            # Check for outliers
            q1 = sorted_data[n // 4]
            q3 = sorted_data[3 * n // 4]
            iqr = q3 - q1
            outlier_count = sum(1 for x in data if x < q1 - 1.5 * iqr or x > q3 + 1.5 * iqr)
            outlier_pct = outlier_count / n

            if outlier_pct > 0.05:
                return ScalingMethod.ROBUST
            return ScalingMethod.ZSCORE

        # Check skewness
        skewness = sum(((x - mean) / stddev) ** 3 for x in data) / n if stddev > 0 else 0

        if abs(skewness) > 1:
            return ScalingMethod.LOG_ZSCORE

        # Check for outliers
        q1 = sorted_data[n // 4]
        q3 = sorted_data[3 * n // 4]
        iqr = q3 - q1
        outlier_count = sum(1 for x in data if x < q1 - 1.5 * iqr or x > q3 + 1.5 * iqr)
        outlier_pct = outlier_count / n

        if outlier_pct > 0.05:
            return ScalingMethod.ROBUST

        return ScalingMethod.ZSCORE

    def normalize_with_params(self, data: List[float], params: ScalingParams) -> List[float]:
        """Apply pre-computed scaling parameters to new data (e.g., test set)."""
        if params.method == ScalingMethod.ZSCORE:
            if params.stddev == 0:
                return [0.0] * len(data)
            return [(x - params.mean) / params.stddev for x in data]
        elif params.method == ScalingMethod.MINMAX:
            range_val = params.max_val - params.min_val
            if range_val == 0:
                return [0.5] * len(data)
            return [(x - params.min_val) / range_val for x in data]
        elif params.method == ScalingMethod.ROBUST:
            if params.iqr == 0:
                return [0.0] * len(data)
            return [(x - params.median) / params.iqr for x in data]
        elif params.method == ScalingMethod.LOG_ZSCORE:
            log_data = [math.log(x) if x > 0 else math.log(1e-10) for x in data]
            if params.stddev == 0:
                return [0.0] * len(log_data)
            return [(x - params.mean) / params.stddev for x in log_data]
        return data

import math
from dataclasses import dataclass
from typing import List, Optional, Tuple


@dataclass
class ControlChartResult:
    """Result of control chart analysis."""
    value: float
    z_score: float
    within_control_limits: bool
    center_line: float
    ucl: float  # upper control limit (+3 sigma)
    lcl: float  # lower control limit (-3 sigma)


@dataclass
class ProcessCapability:
    """Process capability indices."""
    cp: float
    cpk: float
    process_mean: float
    process_stddev: float
    usl: float
    lsl: float
    capability_rating: str
    recommendation: str


class ProcessCapabilityAnalyzer:
    """Statistical process control analysis using z-scores."""

    def __init__(self, usl: float, lsl: float):
        """
        USL: Upper Specification Limit (e.g., max acceptable latency)
        LSL: Lower Specification Limit (e.g., min acceptable throughput)
        """
        self.usl = usl
        self.lsl = lsl

    def calculate_cp(self, stddev: float) -> float:
        """
        Cp = (USL - LSL) / (6 * sigma)
        Measures process spread relative to specification spread.
        Cp > 1: process spread fits within specifications.
        """
        if stddev == 0:
            return float('inf')
        return (self.usl - self.lsl) / (6 * stddev)

    def calculate_cpk(self, mean: float, stddev: float) -> float:
        """
        Cpk = min((USL - mean) / (3 * sigma), (mean - LSL) / (3 * sigma))
        Measures both spread and centering.
        Cpk < 1: process cannot consistently meet specifications.
        """
        if stddev == 0:
            return float('inf')

        cpu = (self.usl - mean) / (3 * stddev)  # upper capability
        cpl = (mean - self.lsl) / (3 * stddev)  # lower capability
        return min(cpu, cpl)

    def analyze(self, data: List[float]) -> ProcessCapability:
        """
        Full process capability analysis.
        """
        n = len(data)
        mean = sum(data) / n
        stddev = math.sqrt(sum((x - mean) ** 2 for x in data) / (n - 1))

        cp = self.calculate_cp(stddev)
        cpk = self.calculate_cpk(mean, stddev)

        if cpk >= 2.0:
            rating = 'Excellent'
            recommendation = 'Process is highly capable. Monitor for drift but no action needed.'
        elif cpk >= 1.33:
            rating = 'Capable'
            recommendation = 'Process meets specifications with margin. Continue monitoring.'
        elif cpk >= 1.0:
            rating = 'Marginal'
            recommendation = 'Process barely meets specifications. Investigate sources of variation.'
        else:
            rating = 'Incapable'
            recommendation = 'Process cannot consistently meet specifications. Reduce variation or adjust specifications.'

        return ProcessCapability(
            cp=round(cp, 4),
            cpk=round(cpk, 4),
            process_mean=round(mean, 4),
            process_stddev=round(stddev, 4),
            usl=self.usl,
            lsl=self.lsl,
            capability_rating=rating,
            recommendation=recommendation,
        )

    def control_chart(self, data: List[float]) -> List[ControlChartResult]:
        """
        Generate control chart analysis for a dataset.
        Flags points outside +/- 3 sigma control limits.
        """
        n = len(data)
        mean = sum(data) / n
        stddev = math.sqrt(sum((x - mean) ** 2 for x in data) / (n - 1))

        ucl = mean + 3 * stddev
        lcl = mean - 3 * stddev

        results = []
        for x in data:
            z = (x - mean) / stddev if stddev > 0 else 0
            results.append(ControlChartResult(
                value=x,
                z_score=round(z, 4),
                within_control_limits=lcl <= x <= ucl,
                center_line=round(mean, 4),
                ucl=round(ucl, 4),
                lcl=round(lcl, 4),
            ))

        return results

    def detect_runs(self, data: List[float], run_length: int = 7) -> List[dict]:
        """
        Detect runs (consecutive points on one side of the mean).
        A run of 7+ points suggests the process has shifted.
        """
        mean = sum(data) / len(data)
        runs = []
        current_run_start = 0
        current_side = 'above' if data[0] > mean else 'below'

        for i in range(1, len(data)):
            side = 'above' if data[i] > mean else 'below'
            if side != current_side:
                run_length_actual = i - current_run_start
                if run_length_actual >= run_length:
                    runs.append({
                        'start_index': current_run_start,
                        'end_index': i - 1,
                        'length': run_length_actual,
                        'side': current_side,
                        'interpretation': (
                            f'Run of {run_length_actual} points {current_side} mean '
                            f'suggests process has shifted. Investigate root cause.'
                        ),
                    })
                current_run_start = i
                current_side = side

        # Check final run
        run_length_actual = len(data) - current_run_start
        if run_length_actual >= run_length:
            runs.append({
                'start_index': current_run_start,
                'end_index': len(data) - 1,
                'length': run_length_actual,
                'side': current_side,
                'interpretation': (
                    f'Run of {run_length_actual} points {current_side} mean '
                    f'suggests process has shifted. Investigate root cause.'
                ),
            })

        return runs

import math
from dataclasses import dataclass
from typing import List, Tuple, Optional


@dataclass
class RobustBaseline:
    """Robust baseline using median and MAD instead of mean and stddev."""
    median: float
    mad: float  # median absolute deviation
    modified_z_threshold: float


class RobustAnomalyDetector:
    """Anomaly detection using robust statistics (median, MAD) instead of mean/stddev."""

    def compute_median(self, data: List[float]) -> float:
        """Compute median of a dataset."""
        sorted_data = sorted(data)
        n = len(sorted_data)
        if n % 2 == 1:
            return sorted_data[n // 2]
        return (sorted_data[n // 2 - 1] + sorted_data[n // 2]) / 2

    def compute_mad(self, data: List[float]) -> float:
        """
        Compute Median Absolute Deviation (MAD).
        MAD = median(|x_i - median(x)|)
        Robust alternative to standard deviation.
        """
        median = self.compute_median(data)
        abs_deviations = [abs(x - median) for x in data]
        return self.compute_median(abs_deviations)

    def compute_modified_zscore(self, x: float, median: float, mad: float) -> float:
        """
        Modified z-score using MAD.
        z_mad = 0.6745 * (x - median) / MAD
        The 0.6745 constant scales MAD to be comparable to stddev for normal data.
        """
        if mad == 0:
            return 0.0
        return 0.6745 * (x - median) / mad

    def detect_outliers_robust(self, data: List[float], threshold: float = 3.5) -> List[dict]:
        """
        Detect outliers using modified z-scores with MAD.
        Threshold of 3.5 is standard for modified z-scores (Iglewicz and Hoaglin).
        """
        median = self.compute_median(data)
        mad = self.compute_mad(data)

        results = []
        for x in data:
            modified_z = self.compute_modified_zscore(x, median, mad)
            results.append({
                'value': x,
                'modified_z_score': round(modified_z, 4),
                'is_outlier': abs(modified_z) > threshold,
                'method': 'MAD-based (robust to outliers)',
            })

        return results

    def iqr_outliers(self, data: List[float], multiplier: float = 1.5) -> List[dict]:
        """
        Detect outliers using Interquartile Range (IQR) method.
        Outliers: x < Q1 - 1.5*IQR or x > Q3 + 1.5*IQR
        """
        sorted_data = sorted(data)
        n = len(sorted_data)
        q1 = sorted_data[n // 4]
        q3 = sorted_data[3 * n // 4]
        iqr = q3 - q1

        lower_fence = q1 - multiplier * iqr
        upper_fence = q3 + multiplier * iqr

        results = []
        for x in data:
            is_outlier = x < lower_fence or x > upper_fence
            results.append({
                'value': x,
                'is_outlier': is_outlier,
                'lower_fence': round(lower_fence, 4),
                'upper_fence': round(upper_fence, 4),
                'method': 'IQR-based',
            })

        return results

    def compare_methods(self, data: List[float]) -> dict:
        """
        Compare z-score, modified z-score (MAD), and IQR outlier detection.
        Shows where methods agree and disagree.
        """
        mean = sum(data) / len(data)
        stddev = math.sqrt(sum((x - mean) ** 2 for x in data) / (len(data) - 1))
        median = self.compute_median(data)
        mad = self.compute_mad(data)

        sorted_data = sorted(data)
        n = len(sorted_data)
        q1 = sorted_data[n // 4]
        q3 = sorted_data[3 * n // 4]
        iqr = q3 - q1

        results = []
        for x in data:
            z = (x - mean) / stddev if stddev > 0 else 0
            modified_z = self.compute_modified_zscore(x, median, mad)
            iqr_outlier = x < q1 - 1.5 * iqr or x > q3 + 1.5 * iqr

            results.append({
                'value': x,
                'z_score': round(z, 4),
                'modified_z_score': round(modified_z, 4),
                'z_outlier': abs(z) > 3,
                'mad_outlier': abs(modified_z) > 3.5,
                'iqr_outlier': iqr_outlier,
                'agreement': 'all' if (abs(z) > 3) == (abs(modified_z) > 3.5) == iqr_outlier else 'disagree',
            })

        return {
            'baseline_stats': {
                'mean': round(mean, 4),
                'stddev': round(stddev, 4),
                'median': round(median, 4),
                'mad': round(mad, 4),
                'q1': round(q1, 4),
                'q3': round(q3, 4),
                'iqr': round(iqr, 4),
            },
            'results': results,
        }