function T = compare_estimators(x, u, n)
% COMPARE_ESTIMATORS  The preferential median next to common alternatives.
%   T = compare_estimators(x, u) returns a struct array with fields name,
%   value, u (NaN where no standard uncertainty formula exists):
%     preferential median (IF&PA), arithmetic mean, median, Hodges-Lehmann,
%     ISO 13528 Algorithm A, and — when all u(k) > 0 — weighted mean (1/u^2),
%     weighted median (1/u^2), DerSimonian-Laird (as in NIST Consensus Builder),
%     largest consistent subset (Cox, 2007) and mixture model mode.
%
%   Needs preferential_median.m in the same folder. MATLAB and GNU Octave,
%   no toolboxes.
%
%   The first row is the preferential median (IF&PA-L); its u is the combined
%   uncertainty — see preferential_median.m for the references and notation.
%
%   Example:
%     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];
%     T = compare_estimators(x, u);
%     for k = 1:numel(T), fprintf('%-32s %.4f\n', T(k).name, T(k).value); end

if nargin < 3, n = []; end
x = double(x(:)');
u = double(u(:)');
if isscalar(u), u = repmat(u, size(x)); end
m = numel(x);

r = preferential_median(x, u, n);
T = struct('name', 'Preferential median (IF&PA)', 'value', r.pm, 'u', r.u_pm);

if m > 1, u_mean = std(x) / sqrt(m); else, u_mean = NaN; end
T(end+1) = struct('name', 'Arithmetic mean', 'value', mean(x), 'u', u_mean);

med = median(x);
if m < 2, u_med = 0; else, u_med = 1.2533 * 1.4826 * median(abs(x - med)) / sqrt(m); end
T(end+1) = struct('name', 'Median', 'value', med, 'u', u_med);

% Hodges-Lehmann: median of all pairwise (Walsh) averages, i <= j
[I, J] = ndgrid(1:m, 1:m);
walsh = (x(I(I <= J)) + x(J(I <= J))) / 2;
T(end+1) = struct('name', 'Hodges-Lehmann', 'value', median(walsh), 'u', NaN);

[xa, ua] = algorithm_a(x);
T(end+1) = struct('name', 'ISO 13528 Algorithm A', 'value', xa, 'u', ua);

if all(u > 0)
    w = 1 ./ u.^2;
    sw = sum(w);
    mean_w = sum(w .* x) / sw;
    T(end+1) = struct('name', 'Weighted mean (1/u^2)', 'value', mean_w, 'u', 1 / sqrt(sw));
    T(end+1) = struct('name', 'Weighted median (1/u^2)', 'value', weighted_median(x, w), 'u', NaN);
    if m >= 2
        % random-effects consensus (DerSimonian-Laird)
        q = sum(w .* (x - mean_w).^2);
        denom = sw - sum(w.^2) / sw;
        if denom > 0, tau2 = max(0, (q - (m - 1)) / denom); else, tau2 = 0; end
        wr = 1 ./ (u.^2 + tau2);
        T(end+1) = struct('name', 'DerSimonian-Laird', 'value', sum(wr .* x) / sum(wr), ...
                          'u', 1 / sqrt(sum(wr)));
        [yl, ul] = largest_consistent_subset(x, u);
        T(end+1) = struct('name', 'Largest consistent subset (Cox)', 'value', yl, 'u', ul);
    end
    T(end+1) = struct('name', 'Mixture model mode', 'value', mixture_mode(x, u), 'u', NaN);
end
end

function v = weighted_median(x, w)
% smallest x with cumulative weight >= half; exactly half -> midpoint
[xs, o] = sort(x);
ws = w(o);
total = sum(ws);
cum = cumsum(ws);
j = find(cum >= total / 2 | abs(cum - total / 2) <= 1e-12 * total, 1);
if abs(cum(j) - total / 2) <= 1e-12 * total && j < numel(xs)
    v = (xs(j) + xs(j + 1)) / 2;
else
    v = xs(j);
end
end

function [xs, u] = algorithm_a(x)
% robust mean by ISO 13528 Algorithm A; u = 1.25 s* / sqrt(m)
m = numel(x);
xs = median(x);
ss = 1.483 * median(abs(x - xs));
if ss == 0 || m < 2, u = 0; return; end
for it = 1:1000
    d = 1.5 * ss;
    c = min(max(x, xs - d), xs + d);
    nx = mean(c);
    ns = 1.134 * sqrt(sum((c - nx).^2) / (m - 1));
    done = abs(nx - xs) <= 1e-10 * max(abs(xs), 1e-300) && abs(ns - ss) <= 1e-10 * ss;
    xs = nx; ss = ns;
    if done || ss == 0, break; end
end
u = 1.25 * ss / sqrt(m);
end

function [y, uy] = largest_consistent_subset(x, u)
% sequential exclusion by chi-square at 0.05 (after Cox, 2007)
keep = 1:numel(x);
while true
    w = 1 ./ u(keep).^2;
    sw = sum(w);
    y = sum(w .* x(keep)) / sw;
    chi2 = sum(w .* (x(keep) - y).^2);
    df = numel(keep) - 1;
    a = 2 / (9 * df);
    q95 = df * (1 - a + 1.6448536269514722 * sqrt(a))^3;   % Wilson-Hilferty
    if numel(keep) <= 2 || chi2 <= q95
        uy = 1 / sqrt(sw);
        return
    end
    d = abs(x(keep) - y) ./ sqrt(max(u(keep).^2 - 1 / sw, 1e-300));
    [~, worst] = max(d);
    keep(worst) = [];
end
end

function t = mixture_mode(x, u)
% mode of the mixture of N(x_k, u_k^2) with equal weights (Ciarlini, Cox, Pavese, Regoliosi, 2004; Duewer, 2004)
f = @(t) sum(exp(-0.5 * ((t - x) ./ u).^2) ./ u);
points = 2001;
lo = min(x - 3 * u);
hi = max(x + 3 * u);
step = (hi - lo) / (points - 1);
fv = arrayfun(f, lo + step * (0:points-1));
[~, i] = max(fv);
i = i - 1;                                   % 0-based, first maximum
a = lo + step * max(i - 1, 0);
b = lo + step * min(i + 1, points - 1);
g = (sqrt(5) - 1) / 2;
c = b - g * (b - a);
d = a + g * (b - a);
for it = 1:100
    if f(c) >= f(d), b = d; else, a = c; end
    c = b - g * (b - a);
    d = a + g * (b - a);
end
t = (a + b) / 2;
end