AvgSpreadBounds

where $[L_x, U_x] = \operatorname{SpreadBounds}(\mathbf{x}, \mathrm{misrate} / 2)$, $[L_y, U_y] = \operatorname{SpreadBounds}(\mathbf{y}, \mathrm{misrate} / 2)$, $w_x = n / (n + m)$, $w_y = m / (n + m)$, and $[L_A, U_A] = [w_x L_x + w_y L_y, w_x U_x + w_y U_y]$.

Robust bounds on $\operatorname{AvgSpread}(\mathbf{x}, \mathbf{y})$ with specified coverage.

Input

• $\mathbf{x} = (x_1, x_2, \ldots, x_n)$ — first sample of measurements, where $n \geq 2$, requires sparity(x)
• $\mathbf{y} = (y_1, y_2, \ldots, y_m)$ — second sample of measurements, where $m \geq 2$, requires sparity(y)
• $\mathrm{misrate} \geq 2 \cdot \max(2^{1-\lfloor n/2 \rfloor}, 2^{1-\lfloor m/2 \rfloor})$ — probability that true avg spread falls outside bounds in the long run

Output

• Value — interval $[L, U]$ bounding $\operatorname{AvgSpread}(\mathbf{x}, \mathbf{y})$
• Unit — same unit as $\mathbf{x}$, $\mathbf{y}$

Notes

• Note — Bonferroni combination of two $\operatorname{SpreadBounds}$ calls with equal split $\mathrm{misrate} / 2$; no independence assumption needed; randomized pairing and cutoff, conservative with ties

Properties

• Symmetry $\operatorname{AvgSpreadBounds}(\mathbf{x}, \mathbf{y}, \mathrm{misrate}) = \operatorname{AvgSpreadBounds}(\mathbf{y}, \mathbf{x}, \mathrm{misrate})$ (equal split)
• Shift invariance adding constants to $\mathbf{x}$ and/or $\mathbf{y}$ does not change bounds
• Scale equivariance $\operatorname{AvgSpreadBounds}(k \cdot \mathbf{x}, k \cdot \mathbf{y}, \mathrm{misrate}) = \lvert k \rvert \cdot \operatorname{AvgSpreadBounds}(\mathbf{x}, \mathbf{y}, \mathrm{misrate})$
• Non-negativity bounds are non-negative
• Monotonicity in misrate smaller $\mathrm{misrate}$ produces wider bounds

Example

• AvgSpreadBounds([1..30], [21..50], 10^(-3)) returns bounds containing AvgSpread

$\operatorname{AvgSpreadBounds}$ provides distribution-free uncertainty bounds for the pooled spread: the weighted average of the two sample spreads. The algorithm computes separate $\operatorname{SpreadBounds}$ for each sample using an equal Bonferroni split and then combines them linearly with weights $n/(n+m)$ and $m/(n+m)$. This guarantees that the probability of missing the true $\operatorname{AvgSpread}$ is at most $\mathrm{misrate}$ without requiring independence between samples.

Minimum misrate — because $\mathrm{misrate} / 2$ must satisfy the per-sample minimum, the overall misrate must be large enough for both samples:

Algorithm

The $\operatorname{AvgSpreadBounds}$ estimator constructs bounds on the pooled spread by combining two independent $\operatorname{SpreadBounds}$ calls through a Bonferroni split.

The algorithm proceeds as follows:

• Equal Bonferroni split — Use $\mathrm{misrate} / 2$ for each per-sample bounds call. Each call uses half the total error budget.

• Per-sample bounds — Compute $[L_x, U_x] = \operatorname{SpreadBounds}(\mathbf{x}, \mathrm{misrate} / 2)$ and $[L_y, U_y] = \operatorname{SpreadBounds}(\mathbf{y}, \mathrm{misrate} / 2)$ (see SpreadBounds).

• Weighted linear combination — With weights $w_x = n / (n + m)$ and $w_y = m / (n + m)$, return:

By Bonferroni’s inequality, the probability that both per-sample bounds simultaneously cover their respective true spreads is at least $1 - \mathrm{misrate}$. Since $\operatorname{AvgSpread}$ is a weighted average of the individual spreads, the linear combination of the bounds covers the true $\operatorname{AvgSpread}$ whenever both individual bounds hold.

using Pragmastat.Algorithms;
using Pragmastat.Exceptions;
using Pragmastat.Functions;
using Pragmastat.Internal;
using Pragmastat.Metrology;

namespace Pragmastat.Estimators;

/// <summary>
/// Distribution-free bounds for AvgSpread via Bonferroni combination of SpreadBounds.
/// Uses equal split alpha = misrate / 2.
/// </summary>
internal class AvgSpreadBoundsEstimator : ITwoSampleBoundsEstimator
{
  internal static readonly AvgSpreadBoundsEstimator Instance = new();

  public Bounds Estimate(Sample x, Sample y, Probability misrate)
  {
    return Estimate(x, y, misrate, null);
  }

  public Bounds Estimate(Sample x, Sample y, Probability misrate, string? seed)
  {
    Assertion.NonWeighted("x", x);
    Assertion.NonWeighted("y", y);
    Assertion.CompatibleUnits(x, y);
    (x, y) = Assertion.ConvertToFiner(x, y);
    // The shuffle runs on the original order; the cached sorted views are sparity-only.
    var rb = EstimateRaw(x.Values, x.SortedValues, y.Values, y.SortedValues, misrate, seed);
    return new Bounds(rb.Lower, rb.Upper, x.Unit);
  }

  /// <summary>
  /// Single shared implementation. <paramref name="x"/>/<paramref name="y"/> are always in
  /// ORIGINAL order (the disjoint-pair shuffle is order-dependent). <paramref name="sortedX"/>/
  /// <paramref name="sortedY"/>, when non-null, are pre-sorted views used only to speed up the
  /// order-independent sparity check. Returns unitless bounds.
  /// </summary>
  internal static Bounds EstimateRaw(
    IReadOnlyList<double> x, IReadOnlyList<double>? sortedX,
    IReadOnlyList<double> y, IReadOnlyList<double>? sortedY,
    double misrate, string? seed)
  {
    Assertion.Validity(x, Subject.X);
    Assertion.Validity(y, Subject.Y);

    if (double.IsNaN(misrate) || misrate < 0 || misrate > 1)
      throw AssumptionException.Domain(Subject.Misrate);

    int n = x.Count;
    int m = y.Count;
    if (n < 2)
      throw AssumptionException.Domain(Subject.X);
    if (m < 2)
      throw AssumptionException.Domain(Subject.Y);

    double alpha = misrate / 2.0;
    double minX = MinAchievableMisrate.OneSample(n / 2);
    double minY = MinAchievableMisrate.OneSample(m / 2);
    if (alpha < minX || alpha < minY)
      throw AssumptionException.Domain(Subject.Misrate);

    if (SpreadBoundsEstimator.SpreadForSparity(x, sortedX) <= 0)
      throw AssumptionException.Sparity(Subject.X);
    if (SpreadBoundsEstimator.SpreadForSparity(y, sortedY) <= 0)
      throw AssumptionException.Sparity(Subject.Y);

    // Validity/domain/sparity already checked above; reuse the inner shuffle directly to avoid
    // re-running SpreadImpl per sample. The shuffle operates on the ORIGINAL order.
    Bounds boundsX = SpreadBoundsEstimator.EstimateInner(x, alpha, seed);
    Bounds boundsY = SpreadBoundsEstimator.EstimateInner(y, alpha, seed);

    double weightX = (double)n / (n + m);
    double weightY = (double)m / (n + m);

    double lower = weightX * boundsX.Lower + weightY * boundsY.Lower;
    double upper = weightX * boundsX.Upper + weightY * boundsY.Upper;
    return new Bounds(lower, upper, MeasurementUnit.Number);
  }
}

Tests

Use an equal Bonferroni split. Compute $[L_x, U_x] = \operatorname{SpreadBounds}(\mathbf{x}, \mathrm{misrate} / 2)$ and $[L_y, U_y] = \operatorname{SpreadBounds}(\mathbf{y}, \mathrm{misrate} / 2)$ using disjoint-pair sign-test inversion (see $\operatorname{SpreadBounds}$). Let $w_x = n / (n + m)$ and $w_y = m / (n + m)$. Return $[L_A, U_A] = [w_x L_x + w_y L_y, w_x U_x + w_y U_y]$.

The $\operatorname{AvgSpreadBounds}$ test suite contains 40 test cases (3 demo + 5 natural + 6 property + 6 edge + 2 distro + 5 misrate + 6 unsorted + 7 error). Since $\operatorname{AvgSpreadBounds}$ returns bounds rather than a point estimate, tests validate that bounds are well-formed and satisfy equivariance properties. Each test case output is a JSON object with lower and upper fields representing the interval bounds. Because $\operatorname{SpreadBounds}$ is randomized, tests fix a seed to make outputs deterministic. Both $\operatorname{SpreadBounds}$ calls use the same seed (two identical RNG streams).

Minimum misrate constraint — the equal split requires

and

,

so

.

Demo examples ($n = m = 30$, $n = m = 20$) — from manual introduction:

• demo-1: $\mathbf{x} = (1, \ldots, 30)$, $\mathbf{y} = (21, \ldots, 50)$, baseline fixture misrate
• demo-2: $\mathbf{x} = (1, \ldots, 30)$, $\mathbf{y} = (21, \ldots, 50)$, stricter fixture misrate, wider bounds
• demo-3: $\mathbf{x} = (1, \ldots, 20)$, $\mathbf{y} = (5, \ldots, 24)$, looser fixture misrate

These cases illustrate how tighter misrates produce wider bounds.

Natural sequences ($\mathrm{misrate}$ varies across achievable fixture levels) — 5 tests:

• natural-10-10: $\mathbf{x} = (1, \ldots, 10)$, $\mathbf{y} = (1, \ldots, 10)$
• natural-10-15: $\mathbf{x} = (1, \ldots, 10)$, $\mathbf{y} = (1, \ldots, 15)$
• natural-15-10: $\mathbf{x} = (1, \ldots, 15)$, $\mathbf{y} = (1, \ldots, 10)$
• natural-15-15: $\mathbf{x} = (1, \ldots, 15)$, $\mathbf{y} = (1, \ldots, 15)$
• natural-20-20: $\mathbf{x} = (1, \ldots, 20)$, $\mathbf{y} = (1, \ldots, 20)$

Property validation ($n = m = 10$) — 6 tests:

• property-identity: $\mathbf{x} = (0, 2, \ldots, 18)$, $\mathbf{y} = (0, 2, \ldots, 18)$, expected output: $[2, 16]$
• property-location-shift: $\mathbf{x} = (10, 12, \ldots, 28)$, $\mathbf{y} = (12, 14, \ldots, 30)$, expected output: $[2, 16]$ (shift invariance)
• property-scale-2x: $\mathbf{x} = (0, 4, \ldots, 36)$, $\mathbf{y} = (4, 8, \ldots, 40)$, expected output: $[4, 32]$ (= 2× identity bounds, scale equivariance)
• property-scale-neg: $\mathbf{x} = (0, -2, \ldots, -18)$, $\mathbf{y} = (-2, -4, \ldots, -20)$, expected output: $[2, 16]$ (= identity bounds, $\lvert k \rvert$ scaling)
• property-symmetry: $\mathbf{x} = (1, \ldots, 10)$, $\mathbf{y} = (6, \ldots, 15)$
• property-symmetry-swapped: $\mathbf{x} = (6, \ldots, 15)$, $\mathbf{y} = (1, \ldots, 10)$, same output as property-symmetry (swap symmetry with equal Bonferroni split)

Edge cases — boundary conditions and extreme scenarios (6 tests):

• edge-small: $n = m = 4$, minimum non-trivial fixture misrate
• edge-negative: $\mathbf{x} = (-10, \ldots, -1)$, $\mathbf{y} = (-20, \ldots, -11)$ (negative values)
• edge-mixed-signs: mixed positive/negative values
• edge-wide-range: powers of 10 from $1$ to $10^9$ (extreme value range)
• edge-asymmetric-8-30: $n = 8$, $m = 30$ (unbalanced sizes)
• edge-duplicates-mixed: $\mathbf{x} = (1, 1, 1, 2, 3, 4)$, $\mathbf{y} = (2, 2, 2, 3, 4, 5)$ (partial ties)

Distribution tests (reference fixture misrates) — 2 tests:

• additive-20-20: $n = m = 20$, $\underline{\operatorname{Additive}}(10, 1)$
• uniform-20-20: $n = m = 20$, $\underline{\operatorname{Uniform}}(0, 1)$

Misrate variation ($\mathbf{x} = (1, \ldots, 20)$, $\mathbf{y} = (5, \ldots, 24)$) — 5 tests spanning progressively stricter fixture misrates:

These tests validate monotonicity: smaller misrates produce wider bounds.

Unsorted tests — verify independent sorting of $\mathbf{x}$ and $\mathbf{y}$ (6 tests):

• unsorted-reverse-x: X reversed, Y sorted
• unsorted-reverse-y: X sorted, Y reversed
• unsorted-reverse-both: both reversed
• unsorted-shuffle-x: X shuffled, Y sorted
• unsorted-shuffle-y: X sorted, Y shuffled
• unsorted-wide-range: wide value range, both unsorted

These tests validate that $\operatorname{AvgSpreadBounds}$ produces identical results regardless of input order.

Error cases — inputs that violate assumptions (7 tests):

• error-empty-x: $\mathbf{x} = ()$ (empty X array) — validity error
• error-empty-y: $\mathbf{y} = ()$ (empty Y array) — validity error
• error-single-element-x: $|\mathbf{x}| = 1$ (too few elements for pairing) — domain error
• error-single-element-y: $|\mathbf{y}| = 1$ (too few elements for pairing) — domain error
• error-constant-x: constant $\mathbf{x}$ violates sparity ($\operatorname{Spread} = 0$)
• error-constant-y: constant $\mathbf{y}$ violates sparity ($\operatorname{Spread} = 0$)
• error-misrate-below-min: misrate below minimum achievable — domain error