Add statistical algorithms with finite-domain and exact-arithmetic tests

README.md

@@ -33,7 +33,7 @@ python -m doctest algorithm_lab/*.py
 Optional installation: `python -m pip install -e .`. For linting and formatting,
 install `.[dev]` and run `ruff check .` and `ruff format --check .`.
 
-## Catalog (148 algorithms)
+## Catalog (156 algorithms)
 
 | Module | Topic | What it demonstrates |
 | --- | --- | --- |
@@ -185,6 +185,14 @@ install `.[dev]` and run `ruff check .` and `ruff format --check .`.
 | [interval_intersection](algorithm_lab/interval_intersection.py) | Intervals | Intersect two sorted disjoint lists of closed integer intervals |
 | [stable_matching](algorithm_lab/stable_matching.py) | Graphs | Gale-Shapley matching with complete strict preferences |
 | [zero_one_bfs](algorithm_lab/zero_one_bfs.py) | Graphs | Shortest distances in directed graphs with zero-or-one edge costs |
+| [percentile](algorithm_lab/percentile.py) | Statistics | Linearly interpolated sample quantiles |
+| [weighted_median](algorithm_lab/weighted_median.py) | Statistics | Lower weighted medians with exact integer weights |
+| [entropy](algorithm_lab/entropy.py) | Statistics | Shannon entropy in bits from nonnegative weights |
+| [jensen_shannon_divergence](algorithm_lab/jensen_shannon_divergence.py) | Statistics | Symmetric Jensen-Shannon divergence in bits |
+| [softmax](algorithm_lab/softmax.py) | Statistics | Stable softmax by subtracting the maximum logit |
+| [isotonic_regression](algorithm_lab/isotonic_regression.py) | Statistics | Exact unweighted isotonic regression with pooled adjacent violators |
+| [kendall_tau](algorithm_lab/kendall_tau.py) | Statistics | Tie-adjusted Kendall rank correlation tau-b |
+| [online_covariance](algorithm_lab/online_covariance.py) | Statistics | One-pass covariance for paired numeric observations |
 
 ## Validation
 

algorithm_lab/_numeric.py

@@ -0,0 +1,25 @@
+"""Finite floating-point validation and scale-safe probability normalization."""
+
+import math
+
+
+def finite(value):
+    if isinstance(value, bool) or not isinstance(value, (int, float)):
+        raise ValueError("values must be finite integers or floats")
+    try:
+        result = float(value)
+    except OverflowError as exc:
+        raise ValueError("value exceeds the floating-point range") from exc
+    if not math.isfinite(result):
+        raise ValueError("values must be finite")
+    return result
+
+
+def probabilities(weights):
+    values = [finite(value) for value in weights]
+    if not values or any(value < 0 for value in values) or max(values) == 0:
+        raise ValueError("weights must be nonnegative with positive total mass")
+    scale = max(values)
+    scaled = [value / scale for value in values]
+    total = math.fsum(scaled)
+    return [value / total for value in scaled]

algorithm_lab/entropy.py

@@ -0,0 +1,18 @@
+"""Shannon entropy in bits from nonnegative weights.
+
+Accept finite int/float weights, excluding bool, with positive total mass.
+Normalize by the largest weight before summing to avoid overflow. Zero-mass
+categories contribute zero. Return entropy in bits. O(n) time and space.
+Extremely small probabilities may underflow to zero in floating point.
+
+>>> entropy([1, 1, 1, 1])
+2.0
+"""
+
+from math import fsum, log2
+
+from algorithm_lab._numeric import probabilities
+
+
+def entropy(weights):
+    return fsum(-p * log2(p) for p in probabilities(weights) if p > 0)

algorithm_lab/isotonic_regression.py

@@ -0,0 +1,25 @@
+"""Exact unweighted isotonic regression with pooled adjacent violators.
+
+Accept int/Fraction values, excluding bool and float. Return nondecreasing
+Fraction fitted values minimizing the sum of squared residuals. Adjacent
+blocks merge only for a strict decrease. Empty input returns []. O(n)
+rational operations and O(n) space; no observations are reordered.
+
+>>> isotonic_regression([3, 1, 2, 5])
+[Fraction(2, 1), Fraction(2, 1), Fraction(2, 1), Fraction(5, 1)]
+"""
+
+from fractions import Fraction
+
+
+def isotonic_regression(values):
+    blocks = []
+    for value in values:
+        if isinstance(value, bool) or not isinstance(value, (int, Fraction)):
+            raise ValueError("values must be integers or Fractions")
+        blocks.append((Fraction(value), 1))
+        while len(blocks) > 1 and blocks[-2][0] * blocks[-1][1] > blocks[-1][0] * blocks[-2][1]:
+            right, count_right = blocks.pop()
+            left, count_left = blocks.pop()
+            blocks.append((left + right, count_left + count_right))
+    return [total / count for total, count in blocks for _ in range(count)]

algorithm_lab/jensen_shannon_divergence.py

@@ -0,0 +1,28 @@
+"""Symmetric Jensen-Shannon divergence in bits.
+
+Inputs are equal-length vectors of finite nonnegative int/float weights,
+excluding bool, each with positive mass. Normalize each independently and
+compare against their equally weighted mixture. Return divergence in [0,1]
+up to floating-point rounding; this is not the square-root distance.
+O(n) time and space. Tiny probabilities may underflow to zero.
+
+>>> jensen_shannon_divergence([1, 0], [0, 1])
+1.0
+"""
+
+from math import fsum, log2
+
+from algorithm_lab._numeric import probabilities
+
+
+def jensen_shannon_divergence(left, right):
+    left, right = probabilities(left), probabilities(right)
+    if len(left) != len(right):
+        raise ValueError("distributions must have equal length")
+    terms = []
+    for p, q in zip(left, right, strict=True):
+        if p:
+            terms.append(p * (log2(p) - log2(p + q) + 1) / 2)
+        if q:
+            terms.append(q * (log2(q) - log2(p + q) + 1) / 2)
+    return max(0.0, min(1.0, fsum(terms)))

algorithm_lab/kendall_tau.py

@@ -0,0 +1,31 @@
+"""Tie-adjusted Kendall rank correlation tau-b.
+
+Inputs are equal-length finite int/float observations, excluding bool, with
+at least two pairs. Return tau-b, accounting separately for ties on each
+axis. A constant axis raises ValueError because the denominator is zero.
+O(n^2) comparisons and O(n) copied storage. Ordering uses validated floats.
+
+>>> kendall_tau([1, 2, 3], [3, 2, 1])
+-1.0
+"""
+
+from math import sqrt
+
+from algorithm_lab._numeric import finite
+
+
+def kendall_tau(left, right):
+    left, right = [finite(x) for x in left], [finite(x) for x in right]
+    if len(left) != len(right) or len(left) < 2:
+        raise ValueError("at least two paired observations are required")
+    concordance, untied_left, untied_right = 0, 0, 0
+    for i in range(len(left)):
+        for j in range(i):
+            a = (left[i] > left[j]) - (left[i] < left[j])
+            b = (right[i] > right[j]) - (right[i] < right[j])
+            concordance += a * b
+            untied_left += a != 0
+            untied_right += b != 0
+    if not untied_left or not untied_right:
+        raise ValueError("correlation is undefined for a constant axis")
+    return concordance / sqrt(untied_left * untied_right)

algorithm_lab/online_covariance.py

@@ -0,0 +1,34 @@
+"""One-pass covariance for paired numeric observations.
+
+observations yields (x,y) pairs of finite int/float values, excluding bool.
+sample is a bool selecting n-1 rather than n as divisor. Require one pair
+for population or two for sample covariance. Return a float; unrepresentable
+intermediate state raises ValueError. O(n) time and O(1) state using a
+paired Welford recurrence; observations need not fit in memory.
+
+>>> online_covariance([(1, 2), (2, 4), (3, 6)], sample=True)
+2.0
+"""
+
+from math import isfinite
+
+from algorithm_lab._numeric import finite
+
+
+def online_covariance(observations, *, sample=False):
+    if not isinstance(sample, bool):
+        raise ValueError("sample must be a boolean")
+    count, mean_x, mean_y, moment = 0, 0.0, 0.0, 0.0
+    for x, y in observations:
+        x, y = finite(x), finite(y)
+        count += 1
+        delta_x, delta_y = x - mean_x, y - mean_y
+        mean_x += delta_x / count
+        mean_y += delta_y / count
+        moment += delta_x * (y - mean_y)
+        if not all(isfinite(v) for v in (mean_x, mean_y, moment)):
+            raise ValueError("covariance exceeds floating-point range")
+    denominator = count - int(sample)
+    if denominator <= 0:
+        raise ValueError("not enough observations")
+    return moment / denominator

algorithm_lab/percentile.py

@@ -0,0 +1,27 @@
+"""Linearly interpolated sample quantiles.
+
+Values are nonempty finite int/float observations excluding bool. q is a
+finite number in [0,1], excluding bool. Return a float at position
+q*(n-1) in sorted data (the common type-7 convention). O(n log n) time and
+O(n) space. Input order is irrelevant and input is unchanged.
+
+>>> percentile([0, 10, 20, 30], 0.25)
+7.5
+"""
+
+from math import floor
+
+from algorithm_lab._numeric import finite
+
+
+def percentile(values, q):
+    values = sorted(finite(value) for value in values)
+    q = finite(q)
+    if not values or not 0 <= q <= 1:
+        raise ValueError("nonempty values and q in [0,1] are required")
+    position = q * (len(values) - 1)
+    index = floor(position)
+    if index == len(values) - 1:
+        return values[index]
+    fraction = position - index
+    return (1 - fraction) * values[index] + fraction * values[index + 1]

algorithm_lab/softmax.py

@@ -0,0 +1,24 @@
+"""Stable softmax by subtracting the maximum logit.
+
+Accept a nonempty iterable of finite int/float logits, excluding bool.
+Return nonnegative probabilities summing to approximately 1. Subtracting
+the maximum avoids exponential overflow; very small tails can underflow
+to zero. O(n) time and space, leaving input unchanged.
+
+>>> softmax([1000, 1000])
+[0.5, 0.5]
+"""
+
+from math import exp, fsum
+
+from algorithm_lab._numeric import finite
+
+
+def softmax(logits):
+    logits = [finite(value) for value in logits]
+    if not logits:
+        raise ValueError("logits must be nonempty")
+    maximum = max(logits)
+    values = [exp(value - maximum) for value in logits]
+    total = fsum(values)
+    return [value / total for value in values]

algorithm_lab/weighted_median.py

@@ -0,0 +1,31 @@
+"""Lower weighted medians with exact integer weights.
+
+Values and weights are equally sized integer iterables, excluding bool.
+Weights are nonnegative with positive total; zero weights have no influence.
+Return the smallest observed value whose cumulative weight reaches half the
+total. This minimizes weighted absolute deviation. O(n log n) time, O(n) space.
+
+>>> weighted_median([10, 1, 5], [1, 2, 1])
+1
+"""
+
+from algorithm_lab._validation import integer
+
+
+def weighted_median(values, weights):
+    values, weights = list(values), list(weights)
+    if len(values) != len(weights):
+        raise ValueError("values and weights must have equal length")
+    for value in values:
+        integer(value)
+    for weight in weights:
+        integer(weight, "weight", minimum=0)
+    total = sum(weights)
+    if not total:
+        raise ValueError("total weight must be positive")
+    cumulative = 0
+    for value, weight in sorted(zip(values, weights, strict=True)):
+        cumulative += weight
+        if cumulative * 2 >= total:
+            return value
+    raise AssertionError("positive total must produce a median")

tests/test_entropy.py

@@ -0,0 +1,19 @@
+import math
+import unittest
+
+from algorithm_lab.entropy import entropy
+
+
+class Tests(unittest.TestCase):
+    def test_known_distributions_and_scaling(self):
+        for n in range(1, 30):
+            self.assertAlmostEqual(entropy([1] * n), math.log2(n))
+        self.assertEqual(entropy([0, 8, 0]), 0)
+        self.assertAlmostEqual(entropy([1, 3]), 0.8112781244591328)
+        self.assertAlmostEqual(entropy([1e308, 1e308]), 1)
+        self.assertAlmostEqual(entropy([1, 3]), entropy([1e100, 3e100]))
+
+    def test_invalid(self):
+        for weights in [[], [0, 0], [-1, 2], [math.nan], [math.inf], [True], [10**1000]]:
+            with self.assertRaises(ValueError):
+                entropy(weights)

tests/test_isotonic_regression.py

@@ -0,0 +1,32 @@
+import unittest
+from fractions import Fraction
+from itertools import product
+
+from algorithm_lab.isotonic_regression import isotonic_regression
+
+
+class Tests(unittest.TestCase):
+    def test_all_contiguous_block_partitions(self):
+        self.assertEqual(isotonic_regression([]), [])
+        for n in range(1, 6):
+            for values in product(range(3), repeat=n):
+                costs = []
+                for mask in range(1 << (n - 1)):
+                    ends = [0] + [i for i in range(1, n) if mask >> (i - 1) & 1] + [n]
+                    candidate = []
+                    for a, b in zip(ends, ends[1:], strict=False):
+                        candidate += [Fraction(sum(values[a:b]), b - a)] * (b - a)
+                    if candidate == sorted(candidate):
+                        costs.append(
+                            sum((a - b) ** 2 for a, b in zip(values, candidate, strict=True))
+                        )
+                actual = isotonic_regression(values)
+                self.assertEqual(actual, sorted(actual))
+                self.assertEqual(
+                    sum((a - b) ** 2 for a, b in zip(values, actual, strict=True)), min(costs)
+                )
+
+    def test_invalid(self):
+        for values in [[True], [1.5]]:
+            with self.assertRaises(ValueError):
+                isotonic_regression(values)

tests/test_jensen_shannon_divergence.py

@@ -0,0 +1,27 @@
+import math
+import unittest
+
+from algorithm_lab.jensen_shannon_divergence import jensen_shannon_divergence
+
+
+class Tests(unittest.TestCase):
+    def test_known_values_and_symmetry(self):
+        self.assertEqual(jensen_shannon_divergence([1, 0], [0, 1]), 1)
+        self.assertEqual(jensen_shannon_divergence([1, 3], [2, 6]), 0)
+        expected = 0.31127812445913283
+        self.assertAlmostEqual(jensen_shannon_divergence([1, 0], [1, 1]), expected)
+        for a, b in [([1, 2, 0], [2, 1, 4]), ([1e308, 0], [1e308, 1e308])]:
+            self.assertAlmostEqual(jensen_shannon_divergence(a, b), jensen_shannon_divergence(b, a))
+        self.assertTrue(math.isfinite(jensen_shannon_divergence([5e-324, 1], [0, 1])))
+
+    def test_invalid(self):
+        for a, b in [
+            ([], []),
+            ([1], [1, 2]),
+            ([0], [1]),
+            ([-1], [1]),
+            ([math.nan], [1]),
+            ([True], [1]),
+        ]:
+            with self.assertRaises(ValueError):
+                jensen_shannon_divergence(a, b)

tests/test_kendall_tau.py

@@ -0,0 +1,27 @@
+import math
+import unittest
+from itertools import permutations
+
+from algorithm_lab.kendall_tau import kendall_tau
+
+
+class Tests(unittest.TestCase):
+    def test_permutation_inversion_oracle(self):
+        for n in range(2, 7):
+            pairs = n * (n - 1) // 2
+            for p in permutations(range(n)):
+                inversions = sum(p[i] > p[j] for i in range(n) for j in range(i + 1, n))
+                self.assertAlmostEqual(kendall_tau(range(n), p), (pairs - 2 * inversions) / pairs)
+        self.assertAlmostEqual(kendall_tau([1, 1, 2], [1, 2, 3]), 2 / math.sqrt(6))
+
+    def test_invalid(self):
+        for a, b in [
+            ([], []),
+            ([1], [1]),
+            ([1, 1], [1, 2]),
+            ([1, 2], [1]),
+            ([True, 2], [1, 2]),
+            ([math.nan, 2], [1, 2]),
+        ]:
+            with self.assertRaises(ValueError):
+                kendall_tau(a, b)