DisparityBounds

where $[L_S, U_S] = \operatorname{ShiftBounds}(\mathbf{x}, \mathbf{y}, \mathrm{misrate}_S)$, $[L_A, U_A] = \operatorname{AvgSpreadBounds}(\mathbf{x}, \mathbf{y}, \mathrm{misrate}_A)$, and $\mathrm{misrate}_S + \mathrm{misrate}_A = \mathrm{misrate}$ (Bonferroni split).

Robust bounds on $\operatorname{Disparity}(\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}$ — probability that true disparity falls outside bounds in the long run (minimum depends on $n$, $m$; see Algorithm)

Output

• Value interval $[L, U]$ bounding $\operatorname{Disparity}(\mathbf{x}, \mathbf{y})$
• Unit dimensionless (spread units)

Notes

• Note Bonferroni split between shift and avg-spread bounds; no independence assumption needed; bounds may be unbounded when pooled spread cannot be certified positive

Properties

• Location invariance $\operatorname{DisparityBounds}(\mathbf{x} + k, \mathbf{y} + k, \mathrm{misrate}) = \operatorname{DisparityBounds}(\mathbf{x}, \mathbf{y}, \mathrm{misrate})$
• Scale invariance $\operatorname{DisparityBounds}(k \cdot \mathbf{x}, k \cdot \mathbf{y}, \mathrm{misrate}) = \operatorname{sign}(k) \cdot \operatorname{DisparityBounds}(\mathbf{x}, \mathbf{y}, \mathrm{misrate})$
• Antisymmetry $\operatorname{DisparityBounds}(\mathbf{x}, \mathbf{y}, \mathrm{misrate}) = -\operatorname{DisparityBounds}(\mathbf{y}, \mathbf{x}, \mathrm{misrate})$ (bounds reversed)
• Monotonicity in misrate smaller $\mathrm{misrate}$ produces wider bounds

Example

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

See also: $\operatorname{Compare2}$ for comparing Disparity against practical thresholds with automatic verdict generation.

Algorithm

The $\operatorname{DisparityBounds}$ estimator constructs bounds on $\operatorname{Disparity}(\mathbf{x}, \mathbf{y})$ by combining $\operatorname{ShiftBounds}$ and $\operatorname{AvgSpreadBounds}$ through a Bonferroni split.

Misrate allocation

The total $\mathrm{misrate}$ budget is split between the shift and avg-spread components. Let $min_S = 2 / \binom{n + m}{n}$ (minimum for $\operatorname{ShiftBounds}$) and $min_A = 2 \cdot \max(2^{1-\lfloor n slash 2 \rfloor}, 2^{1-\lfloor m slash 2 \rfloor})$ (minimum for $\operatorname{AvgSpreadBounds}$). The extra budget beyond the minimums is split equally:

Component bounds

Compute $[L_S, U_S] = \operatorname{ShiftBounds}(\mathbf{x}, \mathbf{y}, \mathrm{misrate}_S)$ and $[L_A, U_A] = \operatorname{AvgSpreadBounds}(\mathbf{x}, \mathbf{y}, \mathrm{misrate}_A)$. By Bonferronis inequality, the probability that both intervals simultaneously contain their respective true values is at least $1 - \mathrm{misrate}_S - \mathrm{misrate}_A = 1 - \mathrm{misrate}$.

Interval division

When $L_A > 0$, the disparity bounds are obtained by dividing the shift interval by the avg-spread interval. Since dividing by a positive interval can flip the ordering depending on the sign of the numerator endpoints, the algorithm computes all four combinations and takes the extremes:

Edge cases

When $L_A = 0$ (the avg-spread interval includes zero), the bounds become partially or fully unbounded depending on the sign of $[L_S, U_S]$:

• $L_S > 0$: $[L_S / U_A, +\infty)$
• $U_S < 0$: $(-\infty, U_S / U_A]$
• $L_S = U_S = 0$: $[0, 0]$
• otherwise: $(-\infty, +\infty)$

When $U_A = 0$ (the avg-spread interval collapses to zero), only the sign of the shift determines the result.

using Pragmastat.Exceptions;
using Pragmastat.Internal;
using Pragmastat.Metrology;

using static Pragmastat.Functions.MinAchievableMisrate;

namespace Pragmastat.Estimators;

/// <summary>
/// Distribution-free bounds for disparity using Bonferroni combination.
/// </summary>
public class DisparityBoundsEstimator : ITwoSampleBoundsEstimator
{
  public static readonly DisparityBoundsEstimator Instance = new();

  /// <summary>
  /// Raw native-array entry point. Returns unitless bounds.
  /// </summary>
  /// <param name="x">First sample values, in the order the disjoint-pair shuffle should run on.</param>
  /// <param name="y">Second sample values, in the order the disjoint-pair shuffle should run on.</param>
  /// <param name="misrate">Misclassification rate.</param>
  /// <param name="assumeSorted">
  /// When true, the inputs are treated as already sorted ascending and the embedded
  /// order-independent sparity and shift-bounds sub-computations skip re-sorting. This flag is
  /// inert (never changes the result) ONLY when the inputs are genuinely sorted: the embedded
  /// shift-bounds consumes the inputs as a sorted view, so on UNSORTED input passing true is
  /// undefined behavior and CAN change the result. The caller is responsible.
  /// </param>
  public Bounds Estimate(double[] x, double[] y, double misrate, bool assumeSorted = false) =>
    EstimateRaw(x, assumeSorted ? x : null, y, assumeSorted ? y : null, misrate, null);

  public Bounds Estimate(double[] x, double[] y, double misrate, string? seed, bool assumeSorted = false) =>
    EstimateRaw(x, assumeSorted ? x : null, y, assumeSorted ? y : null, misrate, seed);

  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/shift-only.
    var rb = EstimateRaw(x.Values, x.SortedValues, y.Values, y.SortedValues, misrate, seed);
    return new Bounds(rb.Lower, rb.Upper, MeasurementUnit.Disparity);
  }

  /// <summary>
  /// Single shared implementation. <paramref name="x"/>/<paramref name="y"/> are always in
  /// ORIGINAL order; <paramref name="sortedX"/>/<paramref name="sortedY"/>, when non-null, are
  /// pre-sorted views used only for the order-independent sparity and shift-bounds sub-computations.
  /// Returns unitless bounds (the Sample path re-attaches the disparity unit).
  /// </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 minShift = TwoSample(n, m);
    double minX = OneSample(n / 2);
    double minY = OneSample(m / 2);
    double minAvg = 2.0 * Math.Max(minX, minY);

    if (misrate < minShift + minAvg)
      throw AssumptionException.Domain(Subject.Misrate);

    double extra = misrate - (minShift + minAvg);
    double alphaShift = minShift + extra / 2.0;
    double alphaAvg = minAvg + extra / 2.0;

    // The Spread > 0 sparity check is performed by AvgSpreadBoundsEstimator below (identical
    // predicate, same Subject.X/Y order). ShiftBounds runs first but cannot throw for these inputs
    // (alphaShift >= the two-sample minimum), so it cannot mask that sparity error.
    // ShiftBounds is order-independent given sorted input; use the sorted views when present.
    Bounds shiftBounds = sortedX != null && sortedY != null
      ? ShiftBoundsEstimator.EstimateRaw(sortedX, sortedY, alphaShift, assumeSorted: true)
      : ShiftBoundsEstimator.EstimateRaw(x, y, alphaShift, assumeSorted: false);
    Bounds avgBounds = AvgSpreadBoundsEstimator.EstimateRaw(x, sortedX, y, sortedY, alphaAvg, seed);

    double la = avgBounds.Lower;
    double ua = avgBounds.Upper;
    double ls = shiftBounds.Lower;
    double us = shiftBounds.Upper;

    if (la > 0.0)
    {
      double r1 = ls / la;
      double r2 = ls / ua;
      double r3 = us / la;
      double r4 = us / ua;
      double lower = Math.Min(Math.Min(r1, r2), Math.Min(r3, r4));
      double upper = Math.Max(Math.Max(r1, r2), Math.Max(r3, r4));
      return new Bounds(lower, upper, MeasurementUnit.Number);
    }

    if (ua <= 0.0)
    {
      if (ls == 0.0 && us == 0.0)
        return new Bounds(0.0, 0.0, MeasurementUnit.Number);
      if (ls >= 0.0)
        return new Bounds(0.0, double.PositiveInfinity, MeasurementUnit.Number);
      if (us <= 0.0)
        return new Bounds(double.NegativeInfinity, 0.0, MeasurementUnit.Number);
      return new Bounds(double.NegativeInfinity, double.PositiveInfinity, MeasurementUnit.Number);
    }

    if (ls > 0.0)
      return new Bounds(ls / ua, double.PositiveInfinity, MeasurementUnit.Number);
    if (us < 0.0)
      return new Bounds(double.NegativeInfinity, us / ua, MeasurementUnit.Number);
    if (ls == 0.0 && us == 0.0)
      return new Bounds(0.0, 0.0, MeasurementUnit.Number);
    if (ls == 0.0 && us > 0.0)
      return new Bounds(0.0, double.PositiveInfinity, MeasurementUnit.Number);
    if (ls < 0.0 && us == 0.0)
      return new Bounds(double.NegativeInfinity, 0.0, MeasurementUnit.Number);

    return new Bounds(double.NegativeInfinity, double.PositiveInfinity, MeasurementUnit.Number);
  }
}

Notes

Width Convergence

The table below shows how $\text{Width} = U - L$ narrows as $N$ grows, for $\mathbf{x} = \mathbf{y} = (1, 1 + 1/(N-1), \ldots, 2)$ ($N$ evenly spaced points on $[1, 2]$) and $\mathrm{misrate} = 10^{-3}$. Dashes indicate $N$ too small to achieve the target misrate.

Tests

The $\operatorname{DisparityBounds}$ test suite contains 39 test cases (3 demo + 5 natural + 6 property + 5 edge + 5 misrate + 2 distro + 6 unsorted + 7 error). Since $\operatorname{DisparityBounds}$ returns bounds rather than a point estimate, tests validate that the bounds contain $\operatorname{Disparity}(\mathbf{x}, \mathbf{y})$ and satisfy equivariance properties. Each test case output is a JSON object with lower and upper fields representing the interval bounds. Because the denominator ($\operatorname{AvgSpreadBounds}$) uses randomized $\operatorname{SpreadBounds}$, tests fix a seed to keep outputs deterministic.

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 (reference fixture misrates) 5 tests:

• natural-10-10: $\mathbf{x} = (1, \ldots, 10)$, $\mathbf{y} = (1, \ldots, 10)$, bounds containing $0$
• 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)$, bounds containing $0$
• natural-20-20: $\mathbf{x} = (1, \ldots, 20)$, $\mathbf{y} = (1, \ldots, 20)$, bounds containing $0$

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

• property-identity: $\mathbf{x} = (0, 2, \ldots, 18)$, $\mathbf{y} = (0, 2, \ldots, 18)$, expected output: $[-1.5, 1.5]$
• property-location-shift: $\mathbf{x} = (10, 12, \ldots, 28)$, $\mathbf{y} = (12, 14, \ldots, 30)$, expected output: $[-2, 1]$
• property-scale-2x: $\mathbf{x} = (0, 4, \ldots, 36)$, $\mathbf{y} = (4, 8, \ldots, 40)$ (= 2× location-shift), expected output: $[-2, 1]$ (scale invariance of disparity)
• property-scale-neg: $\mathbf{x} = (0, -2, \ldots, -18)$, $\mathbf{y} = (-2, -4, \ldots, -20)$ (negated), expected output: $[-1, 2]$ (anti-symmetry under sign flip)
• property-symmetry: $\mathbf{x} = (1, \ldots, 10)$, $\mathbf{y} = (6, \ldots, 15)$, observed bounds
• property-symmetry-swapped: $\mathbf{x}$ and $\mathbf{y}$ swapped, bounds negated (anti-symmetry)

Edge cases boundary conditions (5 tests):

• edge-small: $n = m = 6$ (small samples)
• edge-negative: negative values for both samples
• edge-mixed-signs: mixed positive/negative values
• edge-wide-range: extreme value range
• edge-asymmetric-10-20: $n = 10$, $m = 20$ (unbalanced sizes)

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.

Distribution tests ($\mathrm{misrate}$ varies) 2 tests:

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

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{DisparityBounds}$ produces identical results regardless of input order.

Error cases inputs that violate assumptions (7 tests):

• error-empty-x: $\mathbf{x} = ()$ (empty X array)
• error-empty-y: $\mathbf{y} = ()$ (empty Y array)
• error-single-element-x: $|\mathbf{x}| = 1$ (too few elements for pairing)
• error-single-element-y: $|\mathbf{y}| = 1$ (too few elements for pairing)
• 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