IF&PAинтернет-лаборатория
RUEN
Попробовать

Реализации на пяти языках

Каждая реализация — один файл без внешних зависимостей. В нём преференциальная медиана и функция compare_estimators, которая сразу считает девять других оценок на тех же данных: скопировали — и проверили метод на своих данных против привычных практик. Реализована версия IF&PA-L: правило Борда, без самоуточнения (версии метода).

Проверено. Все пять реализаций дают одинаковые ответы на 430 наборах данных: примеры из статьи-первоисточника, вырожденные и несходящиеся случаи, данные в масштабах от 10−34 до 109. Преференциальная медиана сверена с эталонной программой авторов. Кроме неё, в коде реализованы девять распространённых оценок — от среднего и медианы до DerSimonian–Laird и моды смеси (таблица с первоисточниками); они сверены между языками и с независимыми расчётами по формулам первоисточников. Наборы с ответами можно скачать и проверить ими свою реализацию (см. «Контрольные значения»).

Python 3.8+, только стандартная библиотека. Файл: 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 и девять других оценок
    print(e["name"], e["value"], e["u"])
Весь код 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}")

Обозначения в коде

Имена в коде и обозначения на странице «Метод» и в статьях соотносятся так:

В кодеНа странице «Метод»Смысл
pmxpmпреференциальная медиана
u_pm (uPm в JavaScript)u(xpm)суммарная стандартная неопределённость √(u_grid² + u_sampling²)
u_gridupm = 0,41hсоставляющая от разбиения ДАЗ; в статьях именно она называется upm
u_samplingusвыборочная составляющая 1,2533·1,4826·MAD/√m — дополнение реализаций сайта, в статьях её нет
n, hn, hитоговое число значений ДАЗ и норма разбиения
winnersAbestномера победителей Борда
coverageciсводные индикаторы — сколько интервалов накрывает каждое значение ДАЗ
converged—false, если победители так и не пошли подряд и взят самый широкий блок (n достигло n_max)

Контрольные значения

Проверьте свою копию или собственную реализацию на этих наборах:

НаборПараметрыПМuитоговое n
Пять средств измерений, границы [22,009; 24,991], [25,417; 25,843], [24,565; 24,991], [25,843; 26,269], [24,991; 24,991]n = 1124,9910,55911
Несмежные победители: [4,2; 4,55], [4,7; 5,1]n = 34,650,3664
Межлабораторное сличение (пример 1)авто10,01250,004011

Все десять оценок compare_estimators на примере сличения — проверка остальных оценок:

ОценкаЗначениеu
Преференциальная медиана (IF&PA)10,01250,00401
Среднее арифметическое10,020750,00633
Медиана10,01550,00263
Оценка Ходжеса–Лемана10,01575—
Робастное среднее, алгоритм A10,015870,00268
Взвешенное среднее, веса 1/u²10,040780,00292
Взвешенная медиана, веса 1/u²10,064—
DerSimonian–Laird10,021630,01105
Наибольшее согласованное подмножество10,014220,00428
Мода смеси распределений10,01408—

Все 430 наборов сверки с ответами — ifpa_test_vectors.json (486 КБ): примеры статьи, вырожденные, несходящиеся и случайные наборы, а также те же данные в других единицах (×10−34, ×109, сдвиг на 109): ответ не должен зависеть от единиц измерения. Для каждого — входы x, u, n (null — автоматически) и ответы: pm, u_pm, u_grid, u_sampling, итоговое n_final, converged и значения всех десяти оценок. Допуск — относительный 10−9, для моды смеси 10−6 (она ищется численно).

Что возвращает compare_estimators

ОценкаНеопределённостьИсточник
Преференциальная медиана (IF&PA)√[(0,41h)² + (1,2533·MAD′/√m)²]метод — Muravyov et al., Measurement, 2018; правило Борда — AISP'25, 2025; 0,41h — IJDSA, 2026; выборочная составляющая — дополнение сайта
Среднее арифметическоеs/√mнеопределённость — GUM, JCGM 100:2008, п. 4.2 · открытый текст
Медиана1,2533·MAD′/√mасимптотическая стандартная ошибка медианы для нормального распределения
Оценка Ходжеса–Лемана—Hodges, Lehmann, Ann. Math. Statist., 1963 · открытый текст
Робастное среднее, алгоритм A1,25·s*/√mISO 13528:2022, прил. C (стандарт, платный)
Взвешенное среднее, веса 1/u²1/√ΣwCox, Metrologia, 2002 (процедура A)
Взвешенная медиана, веса 1/u²—Gurwitz, BIT, 1990
DerSimonian–Laird1/√Σw*DerSimonian, Laird, Control. Clin. Trials, 1986; Koepke et al., Metrologia, 2017; NIST Consensus Builder
Наибольшее согласованное подмножество1/√Σw по подмножествуCox, Metrologia, 2007 · открытый текст; последовательное исключение, χ² по Wilson, Hilferty, PNAS, 1931 · открытый текст
Мода смеси распределений—Ciarlini, Cox, Pavese, Regoliosi, Metrologia, 2004 · открытый текст; Duewer, документ CCQM/04-15, 2004

MAD′ = 1,4826 · медиана абсолютных отклонений от медианы. Оценки, которые используют неопределённости, возвращаются, только если все uk > 0.

Числа конфигуратора этим кодом напрямую не повторяются. compare_estimators считает ПМ по интервалам x ± u с разбиением ДАЗ, как в IF&PA-L. Конфигуратор по умолчанию берёт интервалы x ± 2u (расширенная неопределённость) и мелкое разбиение — норма не больше 0,05 медианы u. Чтобы получить ПМ, как в конфигураторе, используйте pm_estimate из benchmark.py (или pmEstimate из benchmark.js):

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

Сравнение на имитированных данных

benchmark.py повторяет конфигуратор применимости: те же условия, тот же генератор случайных чисел, те же числа. Нужен ifpa.py в той же папке.

python benchmark.py --list                     # готовые ситуации
python benchmark.py --preset overconfident     # выбросы с заниженной неопределённостью
python benchmark.py --m 12 --contamination 0.2 --shift 8 --outlier-u 0.3
python benchmark.py --preset bounded --kernel triangular   # другое распределение погрешности
python benchmark.py --preset clean --pm-k 1 --pm-grid auto # ПМ по x ± 1u, разбиение IF&PA-L

Распределения (--kernel): normal, uniform, triangular, laplace, t3, skewed. Настройки ПМ: --pm-k — интервал x ± k·u (по умолчанию 2), --pm-grid — fine (норма h ≤ 0,05 медианы u, по умолчанию), medium (0,2) или auto (как в IF&PA-L).

Для JavaScript — benchmark.js (рядом нужен ifpa.js); в 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 })));

// ПМ, как в конфигураторе: интервалы x ± 2u, мелкое разбиение ДАЗ
const pm = B.pmEstimate(x, u, { pmK: 2, pmGrid: "fine" });   // { value, u }

В браузере файлы подключаются тегами <script> (сначала ifpa.js), объект — window.IFPABench.

Условия использования

Код предоставляется для исследовательских и учебных целей. Если метод или код помогли в вашей работе, пожалуйста, сошлитесь на статью о методе. Условия открытой лицензии будут объявлены дополнительно.