1import json
  2from pathlib import Path
  3import numpy as np
  4
  5
  6SEED = 2161
  7
  8
  9def trim_mean(x, alpha):
 10    x = np.asarray(x)
 11    n = len(x)
 12    k = int(np.floor(alpha * n))
 13    if n - 2 * k <= 0:
 14        raise ValueError("alpha leaves no retained samples")
 15    s = np.sort(x, axis=0)
 16    return s[k:n-k].mean(axis=0)
 17
 18
 19def clustered_mixture(x, alpha, labels):
 20    x = np.asarray(x)
 21    labels = np.asarray(labels)
 22    comps, weights, retained = [], [], []
 23    for c in np.unique(labels):
 24        ids = np.flatnonzero(labels == c)
 25        z = x[ids]
 26        k = int(np.floor(alpha * len(z)))
 27        if len(z) - 2 * k <= 0:
 28            continue
 29        comps.append(trim_mean(z, alpha))
 30        weights.append(len(z) / len(x))
 31        # A conservative row-retention proxy for the covariance estimate.
 32        s = np.sort(z, axis=0)
 33        keep = np.all((z >= s[k]) & (z <= s[len(z)-k-1]), axis=1)
 34        retained.append(int(max(keep.sum(), 2)))
 35    weights = np.asarray(weights, dtype=float)
 36    weights /= weights.sum()
 37    return np.asarray(comps), weights, retained
 38
 39
 40def nearest_center_mse(centers, true_centers):
 41    return float(np.mean([np.min(np.sum((c - true_centers) ** 2, axis=1)) for c in centers]))
 42
 43
 44def mode_data(rng, separation, n_each=100, noise=.18, outlier_fraction=0.0, outlier_scale=4.0):
 45    true = np.array([[-separation / 2, 0.0], [separation / 2, 0.0]])
 46    x = np.vstack([true[0] + rng.normal(0, noise, (n_each, 2)),
 47                   true[1] + rng.normal(0, noise, (n_each, 2))])
 48    labels = np.repeat([0, 1], n_each)
 49    m = int(round(outlier_fraction * len(x)))
 50    if m:
 51        ids = rng.choice(len(x), m, replace=False)
 52        # One-sided adversarial/heavy-tailed contamination.
 53        x[ids] += np.array([outlier_scale, outlier_scale])
 54    return x, labels, true
 55
 56
 57def verify_math(rng):
 58    # Prediction 1: pi_j=n_j/N, hence weights sum to one for every partition.
 59    weight_errors = []
 60    for _ in range(30):
 61        n = rng.integers(2, 200, size=5)
 62        weight_errors.append(abs(np.sum(n / n.sum()) - 1.0))
 63    p1 = {"prediction": "mixture weights sum exactly to one", "max_abs_error": float(max(weight_errors))}
 64
 65    # Prediction 2: with no contamination and symmetric Gaussian modes, local
 66    # trimmed centers remain unbiased (error should not grow with separation),
 67    # while global averaging collapses to the midpoint.
 68    sep_rows = []
 69    for sep in [1., 2., 4., 8.]:
 70        rr = []
 71        for rep in range(15):
 72            x, labels, true = mode_data(rng, sep, n_each=100, noise=.18)
 73            global_center = trim_mean(x, .1)
 74            labels2 = np.argmin(np.stack([np.abs(x[:, 0] + sep / 2), np.abs(x[:, 0] - sep / 2)]), axis=0)
 75            comps, w, _ = clustered_mixture(x, .1, labels2)
 76            rr.append((float(np.mean((global_center - np.mean(true, axis=0)) ** 2)),
 77                       nearest_center_mse(comps, true)))
 78        a = np.mean(rr, axis=0)
 79        sep_rows.append({"separation": sep, "global_midpoint_mse": float(a[0]),
 80                         "clustered_mode_mse": float(a[1])})
 81    p2 = {"prediction": "global midpoint error scales as separation^2; local error is separation-invariant",
 82          "rows": sep_rows}
 83
 84    # Prediction 3: alpha trimming removes a one-sided contamination fraction q
 85    # only when q <= alpha (asymptotically); sweep q around alpha.
 86    trim_rows = []
 87    alpha = .10
 88    for q in [0., .05, .10, .15, .25]:
 89        errs_g, errs_c = [], []
 90        for _ in range(25):
 91            x, labels, true = mode_data(rng, 4., n_each=100, noise=.12,
 92                                         outlier_fraction=q, outlier_scale=5.)
 93            global_center = trim_mean(x, alpha)
 94            labels2 = np.argmin(np.stack([np.abs(x[:, 0] + 2.), np.abs(x[:, 0] - 2.)]), axis=0)
 95            comps, _, _ = clustered_mixture(x, alpha, labels2)
 96            errs_g.append(np.linalg.norm(global_center - np.mean(true, axis=0)))
 97            errs_c.append(nearest_center_mse(comps, true) ** .5)
 98        trim_rows.append({"contamination_q": q, "global_error": float(np.mean(errs_g)),
 99                          "clustered_error": float(np.mean(errs_c)),
100                          "predicted_local_threshold": alpha})
101    p3 = {"prediction": "one-sided trimming is effective through q approximately alpha, then residual outlier bias appears",
102          "rows": trim_rows}
103    return {"weight_normalization": p1, "separation_scaling": p2, "contamination_threshold": p3}
104
105
106def compare_methods(rng):
107    rows = []
108    for q in [0., .05, .10, .15, .25]:
109        g, t, c = [], [], []
110        for _ in range(30):
111            x, labels, true = mode_data(rng, 4., n_each=100, noise=.12,
112                                         outlier_fraction=q, outlier_scale=5.)
113            g.append(np.linalg.norm(x.mean(0) - np.mean(true, axis=0)))
114            t.append(np.linalg.norm(trim_mean(x, .1) - np.mean(true, axis=0)))
115            labels2 = np.argmin(np.stack([np.abs(x[:, 0] + 2.), np.abs(x[:, 0] - 2.)]), axis=0)
116            comps, weights, _ = clustered_mixture(x, .1, labels2)
117            # Report both mode error and mode coverage; the mixture preserves two modes.
118            c.append((nearest_center_mse(comps, true) ** .5, len(comps), float(weights.min())))
119        rows.append({"q": q, "ordinary_mean_error": float(np.mean(g)),
120                     "global_trim_error": float(np.mean(t)),
121                     "clustered_center_error": float(np.mean([z[0] for z in c])),
122                     "cluster_count": float(np.mean([z[1] for z in c])),
123                     "smallest_weight": float(np.mean([z[2] for z in c]))})
124    return rows
125
126
127def main():
128    rng = np.random.default_rng(SEED)
129    result = {"seed": SEED, "verification": verify_math(rng), "comparison": compare_methods(rng)}
130    Path("results.json").write_text(json.dumps(result, indent=2))
131    print(json.dumps(result, indent=2))
132
133
134if __name__ == "__main__":
135    main()