""" 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}")