IF&PAonline laboratory
RUEN
Try it

Implementations in five languages

Each implementation is a single file with no dependencies. It contains the preferential median and compare_estimators, which computes nine other estimators on the same data: copy it and test the method on your data against common practice right away. The version implemented is IF&PA-L: the Borda rule, no self-refinement (versions of the method).

Verified. All five implementations give the same answers on 430 data sets: examples from the original paper, degenerate and non-converging cases, data on scales from 10−34 to 109. The preferential median was checked against the authors' reference program. Besides it, the code implements nine common estimators — from the mean and median to DerSimonian–Laird and the mixture model mode (table with the original sources); they were checked across languages and against independent computations by the formulas of the original sources. The data sets with answers can be downloaded to check your own implementation (see “Reference values”).

Python 3.8+, standard library only. File: ifpa.py.

from ifpa import preferential_median, compare_estimators

x = [10.012, 10.018, 10.009, 10.015, 10.021, 10.011, 10.064, 10.016]
u = [0.010, 0.012, 0.015, 0.008, 0.014, 0.010, 0.004, 0.020]

r = preferential_median(x, u)
print(f"PM = {r['pm']:.4f} ± {r['u_pm']:.4f}")      # PM = 10.0125 ± 0.0040

for e in compare_estimators(x, u):                   # PM and nine other estimators
    print(e["name"], e["value"], e["u"])
Full code of ifpa.py
"""
IF&PA — interval fusion with preference aggregation.
Preferential median (PM) of a set of intervals [x_k - u_k, x_k + u_k].

Pure Python, no dependencies. Python 3.8+.

References:
  Muravyov S.V., Khudonogova L.I., Emelyanova E.Yu. Interval data fusion with
  preference aggregation. Measurement, 2018, 116, 621-630. doi:10.1016/j.measurement.2017.08.045
  Muravyov S.V., Khudonogova L.I., Shabramov V.S., Ignatyev V.D., Andreyev D.I.
  Kemeny and Borda rules in constructing central tendency estimators based on the
  preferential median. Proc. AISP'25, 2025. doi:10.1109/AISP68263.2025.11396113
  (the Borda rule gives the same preferential median as the Kemeny rule).
  Muravyov S.V., Khudonogova L.I., Shabramov V.S., Ignatyev V.D. Comparative studies
  of preferential and weighted medians accuracy in robust estimation of central
  tendency. Int. J. Data Sci. Anal., 2026, 22(1), 262. doi:10.1007/s41060-026-01222-6
  (standard uncertainty of the PM due to the RAV partition: 0.41 h).

Version of the method: IF&PA-L — Borda rule, no self-refinement.
Uncertainty: u_pm returned here is the combined standard uncertainty
  sqrt(u_grid^2 + u_sampling^2), written u(x_pm) on the portal's Method page;
  u_grid = 0.41 h is what the papers call u_pm; u_sampling = 1.2533 * 1.4826 * MAD / sqrt(m)
  is an addition of this implementation, not part of the published method.

Usage:
  >>> from ifpa import preferential_median
  >>> r = preferential_median([4.375, 4.9], [0.175, 0.2], n=3)
  >>> round(r["pm"], 3), r["n"]
  (4.65, 4)
"""
import math

PM_UNCERTAINTY_FACTOR = 0.41   # u_grid = 0.41 h  (h / sqrt(6), triangular distribution)
MEDIAN_SE_FACTOR = 1.2533      # sqrt(pi/2): standard error of the median
MAD_TO_SD = 1.4826             # MAD -> standard deviation for a normal sample
N_MIN_DEFAULT = 11
N_MAX_DEFAULT = 201
FLOAT_EPS = 2.220446049250313e-16   # machine epsilon of double


def _median(values):
    v = sorted(values)
    p = len(v)
    if p == 0:
        raise ValueError("empty sample")
    if p % 2 == 1:
        return v[(p - 1) // 2]
    return (v[p // 2 - 1] + v[p // 2]) / 2.0


def _sampling_uncertainty(centres):
    """Standard error of the median of interval centres, robust scale (MAD)."""
    m = len(centres)
    if m < 2:
        return 0.0
    med = _median(centres)
    scale = MAD_TO_SD * _median([abs(c - med) for c in centres])
    return MEDIAN_SE_FACTOR * scale / math.sqrt(m)


def _coverage(grid, lo, hi, eps):
    """c_i = number of k with lo_k - eps <= a_i <= hi_k + eps, for an ascending grid.
    Two pointers over the sorted bounds: O(n + m log m) instead of O(n m)."""
    left = sorted(l - eps for l in lo)
    right = sorted(r + eps for r in hi)
    m, i, j, out = len(left), 0, 0, []
    for a in grid:
        while i < m and left[i] <= a:
            i += 1
        while j < m and right[j] < a:
            j += 1
        out.append(i - j)   # opened up to a minus closed before a
    return out


def _widest_block(idx, grid, centre):
    """Non-adjacent winners after n_max: the widest block of consecutive indices;
    on a tie — the block whose median is closest to the median of centres."""
    blocks, cur = [], [idx[0]]
    for i in idx[1:]:
        if i == cur[-1] + 1:
            cur.append(i)
        else:
            blocks.append(cur)
            cur = [i]
    blocks.append(cur)
    best, best_key = None, None
    for b in blocks:
        key = (len(b), -abs(_median([grid[i] for i in b]) - centre))
        if best_key is None or key > best_key:
            best, best_key = b, key
    return best


def preferential_median(x, u, n=None, n_min=N_MIN_DEFAULT, n_max=N_MAX_DEFAULT):
    """
    Preferential median of intervals [x_k - u_k, x_k + u_k].

    x : interval centres
    u : half-widths (standard uncertainties); a single number means equal u_k
    n : number of discrete values of the range of actual values (RAV);
        None -> n = max(ceil((a_n - a_1) / min 2u_k), n_min)

    Returns a dict: pm, u_pm (total), u_grid (0.41 h), u_sampling, n, h,
    winners (0-based indices), converged, grid, coverage.
    """
    x = [float(v) for v in x]
    if isinstance(u, (int, float)):
        u = [float(u)] * len(x)
    u = [float(v) for v in u]
    if len(x) != len(u):
        raise ValueError("x and u have different lengths")
    if not x:
        raise ValueError("empty set of intervals")
    if any(v < 0 for v in u):
        raise ValueError("uncertainties must be non-negative")
    if not all(math.isfinite(v) for v in x + u):
        raise ValueError("NaN or infinity in data")

    lo = [a - b for a, b in zip(x, u)]
    hi = [a + b for a, b in zip(x, u)]
    a1, an = min(lo), max(hi)
    centres = [(a + b) / 2.0 for a, b in zip(lo, hi)]

    if an - a1 <= 0:  # all intervals collapse to one point
        return dict(pm=a1, u_pm=0.0, u_grid=0.0, u_sampling=0.0, n=1, h=0.0,
                    winners=[0], converged=True, grid=[a1], coverage=[len(lo)])

    if n is None:
        width_min = min(b - a for a, b in zip(lo, hi))  # the narrowest interval, 2 min u_k
        n = math.ceil((an - a1) / width_min) if width_min > 0 else n_min
        n = max(n, n_min)
    n = max(int(n), 2)

    while True:
        h = (an - a1) / (n - 1)
        grid = [a1 + h * i for i in range(n)]
        # tolerance in the units of the data: a share of the grid step and a few ulps of the values
        eps = max(1e-9 * h, 8 * FLOAT_EPS * max(abs(a1), abs(an)))
        # coverage c_i: how many intervals contain a_i. The Borda score is
        # z_i = sum_k (n - |A_k| - 1) + n c_i, so Borda winners are the values
        # with the largest coverage.
        coverage = _coverage(grid, lo, hi, eps)
        best = max(coverage)
        idx = [i for i, c in enumerate(coverage) if c == best]
        adjacent = all(j - i == 1 for i, j in zip(idx, idx[1:]))
        if adjacent or n >= n_max:
            if not adjacent:
                idx = _widest_block(idx, grid, _median(centres))
            u_grid = PM_UNCERTAINTY_FACTOR * h
            u_samp = _sampling_uncertainty(centres)
            return dict(pm=_median([grid[i] for i in idx]),
                        u_pm=math.hypot(u_grid, u_samp),
                        u_grid=u_grid, u_sampling=u_samp, n=n, h=h,
                        winners=idx, converged=adjacent,
                        grid=grid, coverage=coverage)
        n += 1


# ---------------------------------------------------------------------------
# Alternatives for comparison: common practice for combining x_k +/- u_k
# ---------------------------------------------------------------------------
CHI2_Z95 = 1.6448536269514722   # standard normal 0.95 quantile


def chi2_quantile_95(df):
    """0.95 quantile of chi-square, Wilson-Hilferty approximation (error < 1 %)."""
    a = 2.0 / (9.0 * df)
    return df * (1.0 - a + CHI2_Z95 * math.sqrt(a)) ** 3


def weighted_median(x, w):
    """Weighted median: the smallest x with cumulative weight >= half of the total;
    exactly half -> midpoint with the next value."""
    pairs = sorted(zip(x, w), key=lambda p: p[0])
    total = sum(w)
    cum = 0.0
    for j, (v, wj) in enumerate(pairs):
        cum += wj
        if abs(cum - total / 2) <= 1e-12 * total and j + 1 < len(pairs):
            return (v + pairs[j + 1][0]) / 2.0
        if cum >= total / 2:
            return v
    return pairs[-1][0]


def hodges_lehmann(x):
    """Hodges-Lehmann estimate: median of all pairwise (Walsh) averages, i <= j."""
    m = len(x)
    return _median([(x[i] + x[j]) / 2.0 for i in range(m) for j in range(i, m)])


def algorithm_a(x, tol=1e-10, max_iter=1000):
    """Robust mean and SD by ISO 13528 Algorithm A (Huber-type, cut-off 1.5 s*).
    Returns (x*, s*, u) with u = 1.25 s* / sqrt(m)."""
    m = len(x)
    xs = _median(x)
    ss = 1.483 * _median([abs(v - xs) for v in x])
    if ss == 0 or m < 2:
        return xs, ss, 0.0
    for _ in range(max_iter):
        d = 1.5 * ss
        clipped = [min(max(v, xs - d), xs + d) for v in x]
        new_x = sum(clipped) / m
        new_s = 1.134 * math.sqrt(sum((v - new_x) ** 2 for v in clipped) / (m - 1))
        done = abs(new_x - xs) <= tol * max(abs(xs), 1e-300) and abs(new_s - ss) <= tol * ss
        xs, ss = new_x, new_s
        if done or ss == 0:
            break
    return xs, ss, 1.25 * ss / math.sqrt(m)


def dersimonian_laird(x, u):
    """Random-effects consensus (DerSimonian-Laird), as in NIST Consensus Builder.
    Returns (value, standard uncertainty, tau) where tau is the dark uncertainty."""
    m = len(x)
    w = [1.0 / (v * v) for v in u]
    sw = sum(w)
    mean_w = sum(a * b for a, b in zip(w, x)) / sw
    q = sum(a * (b - mean_w) ** 2 for a, b in zip(w, x))
    denom = sw - sum(a * a for a in w) / sw
    tau2 = max(0.0, (q - (m - 1)) / denom) if denom > 0 else 0.0
    ws = [1.0 / (v * v + tau2) for v in u]
    value = sum(a * b for a, b in zip(ws, x)) / sum(ws)
    return value, 1.0 / math.sqrt(sum(ws)), math.sqrt(tau2)


def largest_consistent_subset(x, u):
    """Weighted mean of the largest consistent subset (after Cox, 2007), sequential
    form: while the chi-square test fails at 0.05, exclude the result with the
    largest |x_k - y| / sqrt(u_k^2 - u(y)^2). Returns (value, u, number excluded)."""
    keep = list(range(len(x)))
    while True:
        w = [1.0 / (u[k] ** 2) for k in keep]
        sw = sum(w)
        y = sum(wk * x[k] for wk, k in zip(w, keep)) / sw
        chi2 = sum(wk * (x[k] - y) ** 2 for wk, k in zip(w, keep))
        if len(keep) <= 2 or chi2 <= chi2_quantile_95(len(keep) - 1):
            return y, 1.0 / math.sqrt(sw), len(x) - len(keep)
        uy2 = 1.0 / sw
        worst = max(keep, key=lambda k: abs(x[k] - y) / math.sqrt(max(u[k] ** 2 - uy2, 1e-300)))
        keep.remove(worst)


def mixture_mode(x, u, points=2001, refine=100):
    """Mode of the mixture of normal densities N(x_k, u_k^2) with equal weights
    (Ciarlini, Cox, Pavese, Regoliosi, 2004; Duewer, 2004): grid search over [min(x - 3u), max(x + 3u)], then
    golden-section refinement around the best grid point."""
    def f(t):
        return sum(math.exp(-0.5 * ((t - a) / b) ** 2) / b for a, b in zip(x, u))
    lo = min(a - 3 * b for a, b in zip(x, u))
    hi = max(a + 3 * b for a, b in zip(x, u))
    step = (hi - lo) / (points - 1)
    best_i, best_f = 0, -1.0
    for i in range(points):
        fi = f(lo + step * i)
        if fi > best_f:
            best_i, best_f = i, fi
    a, b = lo + step * max(best_i - 1, 0), lo + step * min(best_i + 1, points - 1)
    g = (math.sqrt(5.0) - 1.0) / 2.0
    c, d = b - g * (b - a), a + g * (b - a)
    for _ in range(refine):
        if f(c) >= f(d):
            b = d
        else:
            a = c
        c, d = b - g * (b - a), a + g * (b - a)
    return (a + b) / 2.0


# Names are shared by all implementations; the site translates them.
ESTIMATORS = (
    "Preferential median (IF&PA)", "Arithmetic mean", "Median", "Hodges-Lehmann",
    "ISO 13528 Algorithm A", "Weighted mean (1/u^2)", "Weighted median (1/u^2)",
    "DerSimonian-Laird", "Largest consistent subset (Cox)", "Mixture model mode",
)


def compare_estimators(x, u, n=None):
    """
    The preferential median next to common alternatives on the same data.

    Returns a list of dicts {name, value, u}; u is None where the estimator
    has no standard uncertainty formula. Estimators that use u_k (weighted,
    DerSimonian-Laird, consistent subset, mixture mode) need all u_k > 0.
    """
    x = [float(v) for v in x]
    if isinstance(u, (int, float)):
        u = [float(u)] * len(x)
    u = [float(v) for v in u]
    m = len(x)
    r = preferential_median(x, u, n=n)
    out = [dict(name=ESTIMATORS[0], value=r["pm"], u=r["u_pm"])]

    mean = sum(x) / m
    sd = math.sqrt(sum((v - mean) ** 2 for v in x) / (m - 1)) if m > 1 else float("nan")
    out.append(dict(name=ESTIMATORS[1], value=mean, u=sd / math.sqrt(m)))
    out.append(dict(name=ESTIMATORS[2], value=_median(x), u=_sampling_uncertainty(x)))
    out.append(dict(name=ESTIMATORS[3], value=hodges_lehmann(x), u=None))
    xa, _sa, ua = algorithm_a(x)
    out.append(dict(name=ESTIMATORS[4], value=xa, u=ua))

    if all(v > 0 for v in u):
        w = [1.0 / (v * v) for v in u]
        out.append(dict(name=ESTIMATORS[5], value=sum(a * b for a, b in zip(w, x)) / sum(w),
                        u=1.0 / math.sqrt(sum(w))))
        out.append(dict(name=ESTIMATORS[6], value=weighted_median(x, w), u=None))
        if m >= 2:
            value, unc, _tau = dersimonian_laird(x, u)
            out.append(dict(name=ESTIMATORS[7], value=value, u=unc))
            value, unc, _excl = largest_consistent_subset(x, u)
            out.append(dict(name=ESTIMATORS[8], value=value, u=unc))
        out.append(dict(name=ESTIMATORS[9], value=mixture_mode(x, u), u=None))
    return out


if __name__ == "__main__":
    r = preferential_median([4.375, 4.9], [0.175, 0.2], n=3)
    print(f"PM = {r['pm']:.3f} +/- {r['u_pm']:.3f} (n = {r['n']})")   # PM = 4.650

    # interlaboratory comparison where Lab-7 understated its uncertainty
    x = [10.012, 10.018, 10.009, 10.015, 10.021, 10.011, 10.064, 10.016]
    u = [0.010, 0.012, 0.015, 0.008, 0.014, 0.010, 0.004, 0.020]
    for e in compare_estimators(x, u):
        unc = "" if e["u"] is None else f" +/- {e['u']:.4f}"
        print(f"{e['name']:30s} {e['value']:.4f}{unc}")

Notation in the code

Names in the code relate to the notation on the Method page and in the papers as follows:

In the codeOn the Method pageMeaning
pmxpmpreferential median
u_pm (uPm in JavaScript)u(xpm)combined standard uncertainty √(u_grid² + u_sampling²)
u_gridupm = 0.41hcomponent due to the RAV partition; this is what the papers call upm
u_samplingussampling component 1.2533·1.4826·MAD/√m — an addition of the site's implementations, not in the papers
n, hn, hfinal number of RAV values and the partition norm
winnersAbestindices of the Borda winners
coveragecicoverage — how many intervals contain each RAV value
converged—false if the winners never became adjacent and the widest block was taken (n reached n_max)

Reference values

Check your copy or your own implementation on these data sets:

Data setParametersPMufinal n
Five instruments, bounds [22.009, 24.991], [25.417, 25.843], [24.565, 24.991], [25.843, 26.269], [24.991, 24.991]n = 1124.9910.55911
Non-adjacent winners: [4.2, 4.55], [4.7, 5.1]n = 34.650.3664
Interlaboratory comparison (example 1)auto10.01250.004011

All ten estimators of compare_estimators on the comparison example — a check of the other estimators:

EstimatorValueu
Preferential median (IF&PA)10.01250.00401
Arithmetic mean10.020750.00633
Median10.01550.00263
Hodges–Lehmann estimate10.01575—
Robust mean, Algorithm A10.015870.00268
Weighted mean, weights 1/u²10.040780.00292
Weighted median, weights 1/u²10.064—
DerSimonian–Laird10.021630.01105
Largest consistent subset10.014220.00428
Mixture model mode10.01408—

All 430 check data sets with answers — ifpa_test_vectors.json (486 KB): examples from the paper, degenerate, non-converging and random sets, and the same data in other units (×10−34, ×109, a shift by 109): the answer must not depend on the units. Each has inputs x, u, n (null — automatic) and answers: pm, u_pm, u_grid, u_sampling, the final n_final, converged and the values of all ten estimators. Tolerance — relative 10−9, 10−6 for the mixture model mode (it is found numerically).

What compare_estimators returns

EstimatorUncertaintySource
Preferential median (IF&PA)√[(0.41h)² + (1.2533·MAD′/√m)²]method — Muravyov et al., Measurement, 2018; Borda rule — AISP'25, 2025; 0.41h — IJDSA, 2026; sampling component — a site addition
Arithmetic means/√muncertainty — GUM, JCGM 100:2008, clause 4.2 · open text
Median1.2533·MAD′/√masymptotic standard error of the median for a normal distribution
Hodges–Lehmann estimate—Hodges, Lehmann, Ann. Math. Statist., 1963 · open text
Robust mean, Algorithm A1.25·s*/√mISO 13528:2022, Annex C (paid standard)
Weighted mean, weights 1/u²1/√ΣwCox, Metrologia, 2002 (procedure A)
Weighted median, weights 1/u²—Gurwitz, BIT, 1990
DerSimonian–Laird1/√Σw*DerSimonian, Laird, Control. Clin. Trials, 1986; Koepke et al., Metrologia, 2017; NIST Consensus Builder
Largest consistent subset1/√Σw over the subsetCox, Metrologia, 2007 · open text; sequential exclusion, χ² by Wilson, Hilferty, PNAS, 1931 · open text
Mixture model mode—Ciarlini, Cox, Pavese, Regoliosi, Metrologia, 2004 · open text; Duewer, CCQM document CCQM/04-15, 2004

MAD′ = 1.4826 × median absolute deviation from the median. Estimators that use the uncertainties are returned only if all uk > 0.

This code does not reproduce the configurator's numbers directly. compare_estimators computes the PM from the intervals x ± u with the RAV partition of IF&PA-L. By default the configurator takes the intervals x ± 2u (expanded uncertainty) and a fine partition — a norm of at most 0.05 of the median u. To get the PM as in the configurator, use pm_estimate from benchmark.py (or pmEstimate from benchmark.js):

from benchmark import pm_estimate
pm, u = pm_estimate(x, u, {"pm_k": 2, "pm_grid": "fine"})

Comparison on simulated data

benchmark.py reproduces the applicability configurator: same conditions, same random generator, same numbers. Needs ifpa.py in the same folder.

python benchmark.py --list                     # ready situations
python benchmark.py --preset overconfident     # outliers with understated uncertainty
python benchmark.py --m 12 --contamination 0.2 --shift 8 --outlier-u 0.3
python benchmark.py --preset bounded --kernel triangular   # another error distribution
python benchmark.py --preset clean --pm-k 1 --pm-grid auto # PM with x ± 1u, IF&PA-L partition

Distributions (--kernel): normal, uniform, triangular, laplace, t3, skewed. PM settings: --pm-k is the interval x ± k·u (default 2), --pm-grid is fine (norm h ≤ 0.05 of the median u, default), medium (0.2) or auto (as in IF&PA-L).

For JavaScript — benchmark.js (needs ifpa.js next to it); in Node.js:

const B = require("./benchmark.js");

const res = B.run(Object.assign({}, B.PRESETS.overconfident, { replicates: 500, seed: 1 }));
console.table(res.rows.map((r) => ({ name: r.name, rmse: r.rmse, mae: r.mae, p90: r.p90 })));

// PM as in the configurator: intervals x ± 2u, fine RAV partition
const pm = B.pmEstimate(x, u, { pmK: 2, pmGrid: "fine" });   // { value, u }

In a browser the files are included with <script> tags (ifpa.js first); the object is window.IFPABench.

Terms of use

The code is provided for research and teaching. If the method or the code helped your work, please cite the paper describing the method. Open licence terms will be announced separately.