estimators.py
Checks: RULER: The knife: a consistent point inside an inconsistent interval observation; RULER: More data, a worse held-out interval observation; RULER: One bound at two statistics theorem; RULER: The spanning tree is not identified property.
Cited on: RULER.
To run it: save estimators.py; standard library only.
"""estimators.py -- standard Bayes-error estimators scored against an
exact Bayes error, on the threshold family of ceiling.py.
QUESTION. ceiling.py derives the Bayes error of [x >= t] read through
x mod M as a Fraction. That makes the family a RULER: any estimator of
the Bayes error can be scored on its BIAS against the truth and on the
COVERAGE of the 95% interval it states for itself, the fraction of
trials whose interval contains the truth. Where do the textbook
estimators fail, is the failure the ring's or the estimator's, and
what do the distribution-free brackets actually certify?
THE ESTIMATORS. A sample is n pairs (bin r, label y).
RESUBSTITUTION (1/n) sum_r min(n_r0, n_r1).
HELD-OUT fit the majority label per bin on half the sample
(an empty or tied bin guesses 1), score its 0-1
error on the other half.
1-NN the nearest-neighbour error R_NN; on a categorical
channel with every bin populated, a test point's
neighbours are its own bin, so R_NN is the rate at
which a random same-bin training label disagrees.
COVER-HART the bracket [(1 - sqrt(1 - 2 R_NN))/2, R_NN].
DIVERGENCE the Friedman-Rafsky count R of cross-label edges of a
minimum spanning tree of the pooled sample, u =
1 - 2R/n, and the bracket [(1 - sqrt u)/2, (1 - u)/2]
of Berisha, Wisler, Hero and Spanias (IEEE Trans.
Signal Process. 2016, Theorem 2, whose D~_p is u).
Each states the binomial bar 1.96 sqrt(p(1 - p)/m) on the m points
its rate is computed from; a bracket pushes each end out by it
(section_mst prints each tree's raw bracket).
THE ARGUMENT (written before the engine).
(A) THE RING IS A MIXTURE OF BERNOULLIS (property). Drawing x
uniform on Z/N and reading (x mod M, [x >= t]) draws r uniform
and y ~ Bernoulli(eta_r), eta_r = (c - J(r))/c with J(r) the
below-count. So whatever an estimator does on the ring it does
on bare Bernoulli bins with those eta_r; the ring supplies
the placement, eta_r = 1/2 exactly at a knife slice (a modulus
pair (N, M) of ceiling.py with its threshold t, here c = N/M even
and t = N/2) with the truth in exact rational arithmetic.
(B) THE KNIFE. At eta_r = 1/2 on every bin, E min(n_r0, n_r1) =
n_r/2 - E|n_r1 - n_r/2|, and the mean absolute deviation of a
fair binomial is sqrt(n_r/2pi) asymptotically; with n_r = n/M
the resubstitution deficit is (1/2) sqrt(2M/(pi n)), while the
stated half-bar is 1.96/(2 sqrt n). Both are 1/sqrt n, their
ratio sqrt(2M/pi)/1.96 is free of n (1.5766 at M = 15), and
the interval never covers at any n.
(C) THE TRANSIENT. Held-out scores a fitted RULE, whose excess over
the optimum is sum_r P(r)|2 eta_r - 1| P(flip at n_r/2 training
points); at |eta_r - 1/2| = delta that decays like
exp(-2 delta^2 n_r) only once n_r passes 1/delta^2, while the
bar decays as 1/sqrt n from the start. So coverage FALLS with n
until n_r ~ 1/delta^2 and recovers after: at M = 105,
eta = 6/11, 1/delta^2 = 484 per bin, n ~ 50,000.
(D) ONE BOUND AT TWO STATISTICS. On a categorical channel with n >>
M, R_NN -> E[2 eta (1 - eta)] and 2R/n -> the same quantity for
any spanning tree that ignores the labels (C). The population u
is E[(2 eta - 1)^2], and (1 - u)/2 = E[2 eta (1 - eta)]
identically, so the Cover-Hart and divergence brackets are one
bound read through two statistics. R* = 1/2 - E|2 eta - 1|/2, so
the lower end's gap is the Jensen gap of |2 eta - 1|: zero iff
|c - 2J(r)| is constant over r, iff s = 0 or (c odd and
Q = (c-1)/2), the odd-c middle range of ceiling.py.
(E) THE BOUND IS EMPTY ON THE RANGE. At N = 105, M = 15 the range
holds t = 46..59 with R* = 3/7 flat, the exact bracket is
[3/7, 24/49], and the no-evidence value min(t, N - t)/N exceeds
24/49 only at t = 52, 53: the exact bound certifies a lift at 2
of the 14 thresholds, and with a sampling bar at fewer.
(F) THE MST IS NOT IDENTIFIED. Two points of one bin sit at distance
0, so any spanning tree of a bin's clique is part of a minimum
spanning tree, every one genuine, and the source's estimator
assumes a unique tree. A class-ordered chain puts one cross edge
in each two-label bin, so 2R/n <= 2(2M - 1)/n -> 0 and the
bracket collapses to [0, 0]; a spanning tree of the bipartite
graph between a bin's two label classes makes all n_r - 1 edges
cross, 2R/n -> 2 and u -> -1, the bracket at [1/2, 1/2] after
clipping. A tree whose shape ignores the labels has the labels
exchangeable on its n_r - 1 edges, expected cross count
2 n_r0 n_r1 / n_r, the population quantity: label-blindness is
the repair, randomness buys nothing beyond it.
DESIGN. Slices of the threshold family (ceiling.py's notation):
KNIFE N = 210, M = 15, t = 105, R* = 1/2
ODD N = 105, M = 15, t = 45, R* = 3/7
WIDE N = 1155, M = 105, t = 525, R* = 5/11
n in {500, 2000, 8000, 32000}, 40 trials (the 1-NN leg 5 draws at
32,000), fixed seeds. Ring arms draw
x uniform; bare arms draw per-bin binomial counts. The transient's
bare arm runs on to n = 512,000. The exact legs are Fractions over
every threshold cell with N <= 500 and every t.
PREDICTIONS (fixed before the engine).
P0 positive control, two equiprobable bins at eta = (0.2, 0.9),
R* = 3/20, n = 8000: resubstitution and held-out cover in at
least 34 of 40 trials with |bias| < 0.01; the Cover-Hart bracket
excludes the no-evidence 0.45 in at least 34 of 40.
P1 (B) resubstitution at KNIFE, ring and bare: coverage <= 0.10 at
every n, the mean deficit over the mean half-bar in
[-2.0, -1.2] at every n.
P2 (C) held-out at WIDE, ring and bare: coverage >= 0.70 at n = 500,
<= 0.35 at n = 32000; the bare arm >= 0.85 at n = 512,000.
P3 (D) exact: (1 - u)/2 = E[2 eta (1 - eta)] on every slice; the
lower end equals R* iff s = 0 or t in the odd-c middle range.
Sampled: |R_NN - E[2 eta (1 - eta)]| < 0.005 at n = 32000 on
KNIFE, ODD, WIDE.
P4 (E) the exact bracket contains the no-evidence value at 12 of
the 14 thresholds in the range, excluding t = 52, 53.
P5 (F) at n = 8000 on ODD: the class-ordered chain gives
2R/n < 4M/n, the bipartite trees 2R/n > 1.9, a chain in draw
order within 0.02 of 24/49 over 40 trials.
KILL. Any print off its prediction. Only the printed value decides.
FINDINGS. No kill fired. One frozen statement was wrong in its
normalization: (D) and P5 said 2R/n tends to E[2 eta (1 - eta)] under
a label-blind tree; it is R/n (expected cross count 2 n_r0 n_r1 / n_r
per bin sums to n E[2 eta (1 - eta)]), and 2R/n = 1 - u tends to
1 - E[(2 eta - 1)^2]. Also in (C): the rule is fit on n_r/2
points, so the flip probability decays like exp(-delta^2 n_r), not
exp(-2 delta^2 n_r); the crossover scale 1/delta^2 stands. The
tree leg counts within-bin edges only; the at most M - 1 edges
between bins move R/n by at most M/n. The first run printed 2R/n = 0.9769 and FAILED
that check; the check was rewritten to R/n, which printed 0.4885
against 24/49 = 0.4898, and the bracket from the same count,
[0.4241, 0.4885], sits on the exact [3/7, 24/49].
controls eta = (0.2, 0.9): bias +0.0010 and -0.0006, coverage 0.95
and 0.93; Cover-Hart covers 40/40 and excludes 0.45 40/40.
knife ratio -1.605, -1.606, -1.623, -1.612 on the ring and
-1.682, -1.644, -1.462, -1.544 bare, against 1.5766;
coverage 0.00 to 0.03 at every n in both ((B)'s never).
transient held-out at WIDE: ring 0.72, 0.60, 0.12, 0.30 over the
four n, bare 0.88, 0.57, 0.17, 0.17, then 0.95 at 128,000
and 0.90 at 512,000. Resubstitution at M = 15, n = 8000
covers 0.03 at delta = 0, 0.90 at 1/64 (offset 0.72), 0.97 at
1/32 and 0.93 at 1/8.
bound the identity on 40,914 slices; 7,296 tight slices, exactly
s = 0 or the range; the exact range bracket excludes the
no-evidence value only at t = 52, 53; R_NN within 0.0003
of its limit at KNIFE, ODD and WIDE.
tree class-ordered 2R/n = 0.0038 (bracket [0.0009, 0.0019]),
bipartite 1.9963 ([1/2, 1/2]), draw order R/n = 0.4885.
Coverage over 40 trials reads to about a tenth; the directions and the
ratio's constancy are the finding, never the rates.
[Read later, on a code read: the knife dial's check (none at
delta = 0, nominal by delta = 1/8) is gated but was added after the
frozen predictions, and no P states it. (B)'s "never covers" is
false as worded: coverage 0.00 to 0.03.]
RUN RECORD. python estimators.py: 17 of 17 checks, 9.2 s, peak commit
9.2 MB; after an audit (the two tree checks that held by construction
cut, the range bracket's ends checked, the control gating), 16 of 16,
9.1 s, 9.3 MB.
"""
import os
os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
import random
import sys
import time
from fractions import Fraction
from math import sqrt, prod
from itertools import combinations
Z = 1.96
TRIALS = 40
NS = (500, 2000, 8000, 32000)
CHECKS = []
def check(name, ok, detail=""):
CHECKS.append(bool(ok))
tag = "ok " if ok else "FAIL"
print(f" [{tag}] {name}" + (f" ({detail})" if detail else ""))
def section(title):
print()
print(title)
# ------------------------------------------------------------ samplers
def ring_etas(N, M, t):
"""eta_r = P(x >= t | x mod M = r), exact."""
c = N // M
return [Fraction(sum(1 for j in range(c) if r + j * M >= t), c)
for r in range(M)]
def ring_sample(rng, N, M, t, n):
"""Per-bin label counts [n_r0, n_r1] from n uniform draws of x."""
counts = [[0, 0] for _ in range(M)]
for _ in range(n):
x = rng.randrange(N)
counts[x % M][x >= t] += 1
return counts
def bare_sample(rng, etas, n):
"""Per-bin label counts: a multinomial over equiprobable bins, then
a binomial label count in each."""
M = len(etas)
counts = []
left = n
for r in range(M):
nr = left if r == M - 1 else rng.binomialvariate(left, 1 / (M - r))
left -= nr
n1 = rng.binomialvariate(nr, float(etas[r])) if nr else 0
counts.append([nr - n1, n1])
return counts
def bayes(etas):
return sum(min(e, 1 - e) for e in etas) / len(etas)
def nn_limit(etas):
return sum(2 * e * (1 - e) for e in etas) / len(etas)
def bar(p, m):
return Z * sqrt(max(p * (1 - p), 0.0) / m)
# ---------------------------------------------------------- estimators
def resub(counts, n):
est = sum(min(a, b) for a, b in counts) / n
return est, bar(est, n)
def held_out(train, test):
m = sum(a + b for a, b in test)
wrong = 0
for (ta, tb), (a, b) in zip(train, test):
guess_one = tb >= ta # empty or tied guesses 1
wrong += a if guess_one else b
est = wrong / m
return est, bar(est, m)
def one_nn(train, test):
"""Expected error of 'a random same-bin training label'."""
m = sum(a + b for a, b in test)
err = 0.0
for (ta, tb), (a, b) in zip(train, test):
tr = ta + tb
if tr == 0:
err += (a + b) / 2
continue
err += a * tb / tr + b * ta / tr
return err / m, m
def cover_hart(rnn, m):
b = bar(rnn, m)
lo_r = max(rnn - b, 0.0)
lo = (1 - sqrt(max(1 - 2 * lo_r, 0.0))) / 2
hi = min(rnn + b, 0.5)
return lo, hi
def split(rng, counts):
"""Split per-bin counts into two halves, each point a fair coin."""
train, test = [], []
for a, b in counts:
a1 = sum(1 for _ in range(a) if rng.random() < 0.5)
b1 = sum(1 for _ in range(b) if rng.random() < 0.5)
train.append([a1, b1])
test.append([a - a1, b - b1])
return train, test
def run(arm, n, trials, seed, estimator, truth):
rng = random.Random(seed)
est_sum = bar_sum = 0.0
hits = 0
for _ in range(trials):
est, b = estimator(rng, arm(rng, n), n)
est_sum += est
bar_sum += b
hits += abs(est - truth) <= b
return est_sum / trials - truth, bar_sum / trials, hits / trials
def est_resub(rng, counts, n):
return resub(counts, n)
def est_held(rng, counts, n):
train, test = split(rng, counts)
return held_out(train, test)
# ------------------------------------------------------------ sections
def section_controls():
section("CONTROLS")
etas = [Fraction(1, 5), Fraction(9, 10)]
truth = float(bayes(etas))
arm = lambda rng, n: bare_sample(rng, etas, n)
b1, _, c1 = run(arm, 8000, TRIALS, 11, est_resub, truth)
b2, _, c2 = run(arm, 8000, TRIALS, 12, est_held, truth)
check("two bins at eta = (0.2, 0.9): resubstitution and held-out "
"cover, unbiased",
c1 >= 0.85 and c2 >= 0.85 and abs(b1) < 0.01 and abs(b2) < 0.01,
f"bias {b1:+.4f}, {b2:+.4f}; coverage {c1:.2f}, {c2:.2f}")
rng = random.Random(13)
excl = cov = 0
for _ in range(TRIALS):
train, test = split(rng, bare_sample(rng, etas, 8000))
lo, hi = cover_hart(*one_nn(train, test))
excl += not (lo <= 0.45 <= hi)
cov += lo <= truth <= hi
check("Cover-Hart covers and excludes the no-evidence 0.45",
excl >= 34 and cov >= 34, f"{cov}/40 cover, {excl}/40 exclude")
etas = ring_etas(210, 15, 105)
check("the knife slice's eta_r is 1/2 at every bin",
set(etas) == {Fraction(1, 2)})
rng = random.Random(14)
counts = ring_sample(rng, 1155, 105, 525, 200000)
worst = max(abs(a + b - 200000 / 105) / sqrt(200000 / 105)
for a, b in counts)
rate = sum(b for _, b in counts) / 200000
check("ring sampler: bin frequencies and label rate",
worst < 4.5 and abs(rate - 630 / 1155) < 4 * sqrt(0.25 / 200000),
f"worst bin {worst:.2f} sd, label rate {rate:.5f} against "
f"{630 / 1155:.5f}")
def section_knife():
section("THE KNIFE: a consistent point inside an inconsistent interval")
etas = ring_etas(210, 15, 105)
const = sqrt(2 * 15 / 3.141592653589793) / Z
print(f" sqrt(2M/pi)/Z = {const:.4f}")
ok = True
for name, arm in (("ring", lambda rng, n: ring_sample(rng, 210, 15,
105, n)),
("bare", lambda rng, n: bare_sample(rng, etas, n))):
for i, n in enumerate(NS):
bias, hb, cov = run(arm, n, TRIALS, 100 + i, est_resub, 0.5)
ratio = bias / hb
print(f" {name} n = {n:5d}: bias {bias:+.4f}, half-bar "
f"{hb:.4f}, ratio {ratio:+.3f}, coverage {cov:.2f}")
ok &= cov <= 0.10 and -2.0 <= ratio <= -1.2
check("coverage <= 0.10 and ratio in [-2.0, -1.2] at every n, ring "
"and bare", ok)
def section_transient():
section("THE TRANSIENT: more data, a worse held-out interval, then not")
etas = ring_etas(1155, 105, 525)
truth = float(bayes(etas))
series = {}
for name, arm, ns in (
("ring", lambda rng, n: ring_sample(rng, 1155, 105, 525, n),
NS),
("bare", lambda rng, n: bare_sample(rng, etas, n),
NS + (128000, 512000))):
row = []
for i, n in enumerate(ns):
bias, hb, cov = run(arm, n, TRIALS, 200 + i, est_held, truth)
row.append(cov)
print(f" {name} n = {n:6d}: bias {bias:+.4f}, half-bar "
f"{hb:.4f}, coverage {cov:.2f}")
series[name] = row
check("coverage >= 0.70 at n = 500 and <= 0.35 at n = 32000, ring and "
"bare", all(s[0] >= 0.70 and s[3] <= 0.35
for s in series.values()))
check("the bare arm recovers, >= 0.85 at n = 512,000",
series["bare"][-1] >= 0.85)
rng_etas = lambda d: [Fraction(1, 2) + d] * 15
row = []
for d in (Fraction(0), Fraction(1, 64), Fraction(1, 32),
Fraction(1, 8)):
e = rng_etas(d)
arm = lambda rng, n, e=e: bare_sample(rng, e, n)
bias, hb, cov = run(arm, 8000, TRIALS, 300, est_resub,
float(bayes(e)))
u = 2 * float(d) * sqrt(8000 / 15)
row.append(cov)
print(f" knife dial delta = {str(d):5s}, scaled offset {u:.2f}: "
f"coverage {cov:.2f}")
check("the knife's failure is a band near 1/2: none at delta = 0, "
"nominal by delta = 1/8", row[0] <= 0.10 and row[-1] >= 0.85)
def threshold_family(cap):
PR = (2, 3, 5, 7, 11, 13)
out = []
for k in range(2, 7):
for ps in combinations(PR, k):
N = prod(ps)
if N > cap:
continue
for m in range(1, k):
for S in combinations(ps, m):
out.append((N, prod(S)))
return out
def section_bound():
section("ONE BOUND AT TWO STATISTICS")
bad_id = bad_tight = slices = tight = 0
for N, M in threshold_family(500):
c = N // M
for t in range(1, N):
etas = ring_etas(N, M, t)
u = sum((2 * e - 1) ** 2 for e in etas) / M
rstar = bayes(etas)
slices += 1
bad_id += (1 - u) / 2 != nn_limit(etas)
Q, s = divmod(t, M)
is_tight = (1 - 2 * rstar) ** 2 == u
tight += is_tight
bad_tight += is_tight != (s == 0 or (c % 2 == 1
and Q == (c - 1) // 2))
check("(1 - u)/2 = E[2 eta (1 - eta)] on every slice", bad_id == 0,
f"{slices} slices")
check("the lower end equals R* iff s = 0 or the odd-c middle range",
bad_tight == 0, f"{tight} tight slices")
inside, ends = [], []
for t in range(46, 60):
etas = ring_etas(105, 15, t)
u = sum((2 * e - 1) ** 2 for e in etas) / 15
hi = (1 - u) / 2
naive = Fraction(min(t, 105 - t), 105)
inside.append((t, naive <= hi))
# u = 1/49 puts the lower end (1 - sqrt u)/2 at 3/7 = R*
ends.append(u == Fraction(1, 49) and bayes(etas) == Fraction(3, 7))
check("the range's exact bracket is [3/7, 24/49] at every t in it",
all(ends))
excluded = [t for t, c in inside if not c]
check("the exact bracket [3/7, 24/49] contains the no-evidence value "
"at 12 of the 14 thresholds in the range", excluded == [52, 53],
f"excluded at {excluded}")
rng = random.Random(400)
for name, (N, M, t) in (("KNIFE", (210, 15, 105)),
("ODD", (105, 15, 45)),
("WIDE", (1155, 105, 525))):
etas = ring_etas(N, M, t)
gaps = []
for _ in range(5):
train, test = split(rng, bare_sample(rng, etas, 32000))
gaps.append(one_nn(train, test)[0] - float(nn_limit(etas)))
g = sum(gaps) / len(gaps)
print(f" {name}: R_NN - E[2 eta (1 - eta)] = {g:+.4f} "
f"(limit {nn_limit(etas)})")
check(f"1-NN reaches its limit at {name}", abs(g) < 0.005)
def chain_cross(labels):
return sum(1 for a, b in zip(labels, labels[1:]) if a != b)
def section_mst():
section("THE SPANNING TREE IS NOT IDENTIFIED")
etas = ring_etas(105, 15, 45)
target = float(nn_limit(etas))
rng = random.Random(500)
n = 8000
grouped = bipart = drawn = 0.0
for _ in range(TRIALS):
cells = [[] for _ in range(15)]
for _ in range(n):
r = rng.randrange(15)
cells[r].append(int(rng.random() < etas[r]))
g = b = d = 0
for lab in cells:
g += chain_cross(sorted(lab))
d += chain_cross(lab)
if 0 < sum(lab) < len(lab):
b += len(lab) - 1
grouped += 2 * g / n / TRIALS
bipart += 2 * b / n / TRIALS
drawn += 2 * d / n / TRIALS
def bracket(two_r):
u = 1 - two_r
return (1 - sqrt(max(u, 0.0))) / 2, (1 - max(u, 0.0)) / 2
for name, v in (("class-ordered chain", grouped),
("bipartite trees", bipart),
("chain in draw order", drawn)):
lo, hi = bracket(v)
print(f" {name}: 2R/n = {v:.4f}, bracket [{lo:.4f}, {hi:.4f}]")
check("draw order, label-blind: R/n within 0.02 of E[2 eta (1 - eta)]"
" = 24/49", abs(drawn / 2 - target) < 0.02, f"{drawn / 2:.4f}")
def main():
t0 = time.time()
section_controls()
if not all(CHECKS):
print("\nCONTROL FAILED: nothing below is read.")
return 1
section_knife()
section_transient()
section_bound()
section_mst()
print()
n, ok = len(CHECKS), sum(CHECKS)
print(f"{ok} of {n} checks pass, {time.time() - t0:.1f} s")
return 0 if ok == n else 1
if __name__ == "__main__":
sys.exit(main())