mixle.stats.compute.error_receipts module

Optional compensated (Kahan) accumulation, carrying a running numerics-error bound.

Precedent: mixle.inference.precision_plan already tracks a validated summed-LL error band (“~1e-6 relative”) for the reduced-precision fused kernel – a static, offline-verified bound picked once per model/data pair. This module is the dynamic, per-accumulator counterpart: an OPT-IN compensated-summation mode that carries its own running error bound as the computation proceeds, so a partition’s numerics receipt travels with its sufficient statistics through combine() exactly like the statistics themselves.

Design: the receipt an accumulator carries is (abs_total, n) – the running sum of absolute addend magnitudes and the running term count. Both are exactly additive under combine() (no approximation: abs_total_combined = abs_total_a + abs_total_b, n_combined = n_a + n_b, by definition), so they satisfy the roadmap’s “bounds ADD through combine() like the stats” literally – they are themselves additive sufficient statistics. The error BOUND is then the pure function error_bound() of those two additive quantities. This is deliberately the SAFE direction: for the naive bound, recomputing from the merged (n, abs_total) is always >= summing the two children’s bounds separately (since (n_a + n_b - 1) >= (n_a - 1) + (n_b - 1) whenever n_a, n_b >= 1), so composing through the additive receipt never under-reports the merged error.

Bound formulas (Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., secs 4.2 / 8.1):

  • naive summation of n terms: |error| <= (n-1) * eps * sum(|x_i|) (eq. 8.13 style)

  • Kahan compensated summation: |error| <= (2*eps + n*eps**2) * sum(|x_i|) (eq. 8.15 style)

where eps is the float64 machine epsilon (2**-52). Both are asymptotic first/second-order bounds under the idealized no-overflow/no-underflow model; they are verified as non-violated upper bounds (against a high-precision math.fsum reference) in mixle/tests/numerics_error_receipts_test.py.

error_bound(n, abs_total, compensated)[source]

Return a valid upper bound on the float64 summation round-off error.

Parameters:

  • n (int) – Number of terms summed.

  • abs_total (float) – Running sum of the absolute value of each (weighted) addend.

  • compensated (bool) – Whether the sum was accumulated with Kahan compensation.

Returns:

Non-negative upper bound on |computed_sum - true_sum|.

Return type:

float

class CompensatedAccumulator(total=0.0, compensation=0.0, abs_total=0.0, n=0, compensated=True)[source]

Bases: object

A running (optionally Kahan-compensated) sum plus the additive receipt its bound derives from.

total is the running sum (Kahan-corrected when compensated=True, plain float64 accumulation otherwise). abs_total and n are the additive receipt: they compose exactly under combine(), and bound() derives the error bound from them via error_bound().

Parameters:
add(x, weight=1.0)[source]

Fold one (weighted) addend into the running sum.

Parameters:
Return type:

CompensatedAccumulator

combine(other)[source]

Merge another partition’s running sum + receipt into this one, in place.

The other accumulator’s compensation-corrected estimate (other.total - other.compensation) is folded in as a single addend via the same (compensated or plain) addition rule this accumulator uses; abs_total/n – the receipt – add exactly, the same way the underlying sufficient statistics (sum, count, …) do.

Parameters:

other (CompensatedAccumulator)

Return type:

CompensatedAccumulator

bound()[source]

Return the current error-bound receipt (see error_bound()).

Return type:

float

kahan_reduce(values, weights=None)[source]

Reduce a sequence (optionally weighted) to a fresh CompensatedAccumulator.

A convenience one-shot reducer, mainly used by tests and small ad hoc reductions; accumulators that need this incrementally (across update/seq_update calls) should hold their own CompensatedAccumulator instance instead of re-reducing from scratch each time.

Parameters:
Return type:

CompensatedAccumulator

class ConditioningReceipt(eigenvalues, condition_number, near_degenerate, degenerate_ratio_threshold)[source]

Bases: object

A real, computed numerical-conditioning diagnostic for a (multivariate) fit.

Captures the covariance eigenvalue spectrum an estimate() call saw, so a caller can tell a healthy fit from one balanced on a near-degenerate direction without recomputing the eigenspectrum itself.

Parameters:

  • eigenvalues (ndarray)

  • condition_number (float)

  • near_degenerate (bool)

  • degenerate_ratio_threshold (float)

eigenvalues

Eigenvalues of the (raw, pre-regularization) covariance, ascending.

Type:

np.ndarray

condition_number

max_eigenvalue / min_eigenvalue (inf if the smallest eigenvalue is <= 0, i.e. the raw covariance is singular / numerically indefinite).

Type:

float

near_degenerate

True when the smallest-to-largest eigenvalue ratio falls below degenerate_ratio_threshold (or the smallest eigenvalue is non-positive).

Type:

bool

degenerate_ratio_threshold

The ratio threshold used to set near_degenerate.

Type:

float

to_dict()[source]

JSON-friendly view of this receipt.

Return type:

dict[str, Any]

conditioning_receipt(covar, degenerate_ratio=1.0e-6)[source]

Compute a ConditioningReceipt from a (symmetric) covariance matrix.

degenerate_ratio is the smallest/largest eigenvalue ratio below which the covariance is flagged near-degenerate (i.e. it has a direction of near-zero variance relative to the dominant one) – the “near-degenerate variance flag” the roadmap calls for, generalized from a single variance to the eigenvalue spectrum for the multivariate case.

Parameters:
Return type:

ConditioningReceipt