How a Synthetic Control Manufactured Shared Error

Technical note · precip-disagreement project · 15 October 2026 · reproducible on request

Summary. A synthetic control for triple collocation on structural rain-field coordinates simulated coarse observers by widening Gaussian rain cells while holding their peaks fixed. That operation does not conserve mass, and the mass gain depends on a property of the truth (cell width) that every observer shares. Two observers with independent planted errors therefore acquired correlated coordinate errors (about 0.16 on centroid y in the original 15 × 15 windows), and the collocation estimator mis-ranked observers on dispersion. A finding built on this baseline correlation was later withdrawn. A mass-conserving Gaussian convolution removes the spurious centroid correlation on the full canvas and restores the correct rankings; a small dispersion correlation remains, of unresolved origin. The note gives the mechanism, the invariant checks that would have caught it, a minimal reproduction, and a short checklist.

1. The setup

Triple collocation estimates each observer's noise fraction from the covariances of three observers of the same quantity, under the assumption that their errors are mutually independent. In the project the quantities are structural coordinates of a rain field: the intensity-weighted centroid (y, x) and the intensity-weighted dispersion about it. The synthetic control (T1) existed to check that the estimator recovers planted noise when every assumption holds.

Each truth field is a sum of 1 + Poisson(1.5) Gaussian cells, with centres uniform in ±5 cells, widths w ~ U(0.8, 2.5) and peaks lognormal(1.5, 0.8). Three observers re-render the same cells with independent per-cell position and intensity errors. The radar-like observer R is sharp; the microwave-like S1 and infrared-like S2 are coarser, blurred by b = 1 and b = 2 cells. No shared random perturbation was planted between observers.

2. The symptom

With nothing shared planted, the residuals of S1 and S2 on truth were correlated: 0.16 on centroid y and x, 0.34 on dispersion. The collocation estimate for dispersion ranked the sharp radar-like observer as noisier than S1 (estimate 0.121 vs 0.058; truth 0.046 vs 0.137). The baseline correlation was read at the time as a finding: smoothing itself creates shared structural error between coarse observers (F6). It had a plausible physical reading, and no invariant check on the renderer was in place to contradict it.

3. The mechanism

The coarse observers were produced by re-rendering each cell with width w + b and the same peak I. A two-dimensional Gaussian cell has mass

M = ∫ I exp(−r² / 2w²) dA = 2π I w²,

so widening multiplies each cell's mass by

M′ / M = ((w + b) / w)².

For b = 2 this factor is 12.3 for the narrowest cell (w = 0.8) and 3.2 for the widest (w = 2.5). A field-level centroid is the mass-weighted mean of cell centres, Σ Mk ck / Σ Mk. Widening reweights it toward narrow cells. Because cell width belongs to the truth, both coarse observers are pulled toward the same cells in the same window, by amounts that depend on the configuration of cells rather than linearly on the true centroid. The pull survives regression on truth and appears as correlated residuals. It needs at least two cells: a lone symmetric cell keeps its centre when widened. In minimal.py, single-cell fields give a centroid-y correlation of −0.015 under widening and multi-cell fields +0.099. Dispersion is affected twice: through the same reweighting, and because widths add linearly (w + b) where a convolution adds them in quadrature (√(w² + b²)).

Panel a: one-dimensional profiles of a narrow and a wide Gaussian cell, the widened versions keep peak 1 and grow, the convolved versions flatten. Panel b: mass multiplier against cell width, falling from 12.3 to 3.2 for b equal 2 and from 5.1 to 2.0 for b equal 1, flat at 1 for convolution.
Figure 1. (a) A narrow (solid) and a wide (dashed) cell under each blur, b = 2. Widening keeps the peak and adds mass; convolution keeps the mass and lowers the peak. (b) The widening mass multiplier falls with cell width, so it reweights every coarse observer toward the same narrow cells.

4. The invariant checks

A blur that models a coarser sensor should satisfy three properties that can be computed on the truth field alone, before any observer error is added. For a normalised Gaussian kernel of standard deviation b on a canvas large enough to hold the field:

  1. mass is conserved (∫ f∗g = ∫ f · ∫ g = ∫ f);
  2. the centroid is unchanged (first moments add, and the kernel's is zero);
  3. squared dispersion grows by exactly 2b² (second central moments add, b² per axis).

On 1,000 truth fields from minimal.py:

blurbmedian mass ratiomedian centroid shift (cells)median |Δd² − 2b²| / 2b²verdict
widen (T1)12.520.15350%fail
widen (T1)24.750.22171%fail
convolve (T1-bis)11.0005 × 10⁻¹⁶0.01%pass
convolve (T1-bis)21.0005 × 10⁻¹⁶0.03%pass

These checks are the stronger reference. They do not depend on any statistical estimator, and a failure localises the problem to the renderer. In the project, T1-bis made them a pass/fail gate before any collocation was computed (centroid shift 5 × 10⁻¹⁶, squared-dispersion growth within 0.04%). The original flaw was identified by an external review of the project and confirmed against the simulation code.

5. The corrected result

With the same random cells and the mass-conserving blur, the S1–S2 centroid error correlation on the full canvas falls to the sampling-noise level (|r| ≤ 0.014, 95% field-bootstrap intervals including zero), and the collocation estimates for centroid y and dispersion lie close to the truth with the correct ranking (T1-bis showed the same for all four structural variables). A small dispersion correlation remains (sections 6 and 7).

Grouped bars of S1 to S2 error correlation. Full canvas: widen 0.07, 0.06, 0.29; convolve 0.00, minus 0.01, 0.03. Central window: widen 0.15, 0.13, 0.33; convolve 0.00, minus 0.01, 0.06.
Figure 2. Correlation between S1 and S2 coordinate errors with no shared error planted (10,000 fields). No shared perturbation was planted, so values whose intervals exclude zero come from the simulation itself.
Two scatter panels of collocation estimate against true noise fraction. Centroid y: both blurs near the identity line. Dispersion: under widening R is overestimated at 0.11 against 0.045 and S1 underestimated at 0.07 against 0.14; under convolution all three lie on the line.
Figure 3. Triple-collocation noise fraction against the true value, full canvas. Under widening, dispersion is mis-ranked (R estimated noisier than S1); under convolution all three observers are recovered.
full canvas, dispersionR est. / trueS1 est. / trueS2 est. / trueranking
widen0.110 / 0.0450.071 / 0.1380.250 / 0.304S1 < R < S2 (wrong)
convolve0.051 / 0.0450.102 / 0.1090.280 / 0.284R < S1 < S2 (correct)

The corrected control still does its job as a positive control: in T1-bis, a shared displacement planted in S1 and S2 (σ = 1.0 cell) inverts the centroid ranking as expected, and a five-observer tetrad audit identifies the planted pair in 10 of 10 blocks with no false alarms when nothing is planted.

6. The residual window effect

Real data are read in finite windows. A plausible reason windows add shared error is that they clip the tails of smoother fields more than sharp ones (not tested separately here). This is a genuine, smaller source of shared error that remains under a correct blur. In T1-bis, 15 × 15 windows gave S1–S2 error correlations up to 0.02 on centroids and up to 0.14 on dispersion, growing as the window was shifted off the storm (central window 0.06; shifted by (−3, +4) cells 0.14). minimal.py reads only the central window and shows the same pattern (centroids |r| ≤ 0.01, dispersion 0.06). Widening and windowing also interact: under widening the centroid correlation is 0.07 on the full canvas and 0.15 in the window, close to the original 0.16, which was measured in windows.

7. Differences between the minimal example and the project runs

8. A checklist for synthetic benchmarks

Two validations are separate tasks: that the rendering operation does what it claims (invariant checks), and that the extracted measurements meet the statistical assumptions of the estimator (independence and linearity of their errors). This note fails the first and, after the correction, shows a small residual in the second.

  1. Write down the invariants of each simulated operation (what each should preserve, and what it should change by a known amount: for a blur, mass and first moment preserved, second moment increased by 2b²) and test them on the truth alone, as a gate, before any downstream statistic is computed.
  2. Implement degradations as operators on the rendered field (convolution, resampling), not as edits to the generating parameters, unless the parameter edit is shown to equal the operator.
  3. Check that independence holds where it is assumed. Compute the true pairwise error correlations among simulated observers with nothing shared planted; they should be near zero. A non-zero value is a property of the simulator until shown otherwise.
  4. Look for common directional bias, not just common sensitivity. Two observers can both become noisier on small cells without their signed errors correlating. Shared error arises when both shift the estimate in the same direction as a function of the truth, as here, where both reweight the scene toward narrow cells.
  5. Separate the canvas from the read region. Simulate on a canvas large enough to contain the field, then apply windows, thresholds and clipping as separate, measured steps.
  6. Keep a planted positive control alongside the null, so a corrected simulator is shown to still detect what it should.
  7. Treat a surprising baseline result as a test of the simulator first. A physically plausible explanation is not evidence that the simulation produced it for that reason.

Reproduce

cd precip-disagreement/notes/synthetic-shared-error
../../.venv/bin/python minimal.py            # numbers, about 30 s
../../.venv/bin/python minimal.py --figures  # also writes fig1-3 (matplotlib)

Sources: T1_PLAN.md, T1_RESULTS.md, T1BIS_PLAN.md, T1BIS_RESULTS.md, scripts/t1_synthetic.py, scripts/t1bis.py, LABBOOK.html §7 and §12.

Appendix: minimal.py

Python 3 with numpy and scipy (matplotlib only for --figures; the figure code is in figures.py). Seed 20261015 for the fields, 7 for the bootstrap; 10,000 fields; about 30 s. The full source follows so that this file is self-contained.

"""Minimal reproduction: a non-mass-conserving blur manufactures shared error between simulated observers.

Same random rain cells, two ways to 'blur' a coarse observer:
  widen    : re-render each Gaussian cell with width w + b and the same peak (as in the original T1 control)
  convolve : render at width w, then convolve with a normalised Gaussian of sd b (mass-conserving, as in T1-bis)
Run:  python minimal.py            (numbers only)
      python minimal.py --figures  (also writes fig1-3 PNGs; needs matplotlib)
"""
import sys
import numpy as np
from scipy.ndimage import gaussian_filter

N, SEED, C = 10000, 20261015, 64
YY, XX = np.mgrid[0:C, 0:C] - (C - 1) / 2
OBS = {"R": (0.5, 0.0, 0.3), "S1": (0.8, 1.0, 0.4), "S2": (1.5, 2.0, 0.6)}  # position sd, blur b, log-intensity sd
WIN = (slice(C // 2 - 7, C // 2 + 8),) * 2  # central 15 x 15 window


def render(cy, cx, w, I):
    f = np.zeros((C, C))
    for a, b, s, i in zip(cy, cx, w, I):
        f += i * np.exp(-((YY - a) ** 2 + (XX - b) ** 2) / (2 * s ** 2))
    return f


def observe(cells, b, mode):
    cy, cx, w, I = cells
    if mode == "widen":
        return render(cy, cx, w + b, I)
    f = render(cy, cx, w, I)
    return gaussian_filter(f, b, mode="constant", truncate=4.0) if b > 0 else f


def reads(f, region=(slice(None), slice(None))):
    A, Y, X = f[region], YY[region], XX[region]
    m = A.sum()
    cy, cx = (A * Y).sum() / m, (A * X).sum() / m
    return m, cy, cx, np.sqrt((A * ((Y - cy) ** 2 + (X - cx) ** 2)).sum() / m)


def simulate():
    rng = np.random.default_rng(SEED)
    out = []
    for _ in range(N):
        K = 1 + rng.poisson(1.5)
        truth = (rng.uniform(-5, 5, K), rng.uniform(-5, 5, K), rng.uniform(0.8, 2.5, K), rng.lognormal(1.5, 0.8, K))
        obs = {o: (truth[0] + rng.normal(0, sd, K), truth[1] + rng.normal(0, sd, K), truth[2], truth[3] * rng.lognormal(0, li, K))
               for o, (sd, _, li) in OBS.items()}
        out.append((truth, obs))
    return out


def table(W, mode, region):
    """rows = windows; dict of observer -> array (N, 4) of (mass, cy, cx, d)."""
    T = np.array([reads(render(*t), region) for t, _ in W])
    return {"T": T, **{o: np.array([reads(observe(ob[o], OBS[o][1], mode), region) for _, ob in W]) for o in OBS}}


def errcorr(D, a, b, j):
    res = [D[o][:, j] - np.polyval(np.polyfit(D["T"][:, j], D[o][:, j], 1), D["T"][:, j]) for o in (a, b)]
    return np.corrcoef(*res)[0, 1]


def errcorr_ci(D, a, b, j, n_boot=1000, seed=7):
    """95% percentile interval for errcorr, resampling whole synthetic fields (rows) with replacement."""
    rng = np.random.default_rng(seed); n = len(D["T"])
    bs = []
    for _ in range(n_boot):
        k = rng.integers(0, n, n)
        bs.append(errcorr({o: v[k] for o, v in D.items()}, a, b, j))
    return np.percentile(bs, [2.5, 97.5])


def tc(D, j, trip=("R", "S1", "S2")):
    Cm = np.cov(np.array([D[o][:, j] for o in trip]))
    est = {o: 1 - Cm[i, (i + 1) % 3] * Cm[i, (i + 2) % 3] / (Cm[i, i] * Cm[(i + 1) % 3, (i + 2) % 3]) for i, o in enumerate(trip)}
    true = {o: 1 - np.corrcoef(D[o][:, j], D["T"][:, j])[0, 1] ** 2 for o in trip}
    return est, true


if __name__ == "__main__":
    W = simulate()
    print("1. Analytic checks on the truth field (first 1000 windows), blur b = 1 and 2")
    for mode in ("widen", "convolve"):
        for b in (1.0, 2.0):
            r = np.array([(*reads(observe(t, b, mode)), *reads(render(*t))) for t, _ in W[:1000]])
            mass = np.median(r[:, 0] / r[:, 4])
            shift = np.median(np.hypot(r[:, 1] - r[:, 5], r[:, 2] - r[:, 6]))
            d2 = np.median(np.abs(r[:, 3] ** 2 - r[:, 7] ** 2 - 2 * b ** 2) / (2 * b ** 2))
            ok = abs(mass - 1) < 1e-6 and shift < 0.01 and d2 < 0.01
            print(f"   {mode:8s} b={b:.0f}: mass ratio {mass:6.3f}  centroid shift {shift:.1e}  "
                  f"|d2 growth - 2b^2|/2b^2 {d2:6.2%}  -> {'PASS' if ok else 'FAIL'}")
    tabs = {(m, g): table(W, m, reg) for m in ("widen", "convolve") for g, reg in (("full", (slice(None),) * 2), ("15x15", WIN))}
    print("\n2. S1-S2 coordinate-error correlation with nothing shared planted (should be ~0)")
    print("   blur      region   centroid y  centroid x  dispersion   (95% field-bootstrap intervals, 1000 resamples)")
    for (m, g), D in tabs.items():
        cells = []
        for j in (1, 2, 3):
            lo, hi = errcorr_ci(D, "S1", "S2", j)
            cells.append(f"{errcorr(D, 'S1', 'S2', j):+.3f} [{lo:+.3f}, {hi:+.3f}]")
        print(f"   {m:8s}  {g:6s}   " + "  ".join(cells))
    print("\n3. Triple-collocation noise fractions (estimate / truth), full canvas")
    for m in ("widen", "convolve"):
        for j, v in ((1, "centroid y"), (3, "dispersion")):
            est, true = tc(tabs[(m, "full")], j)
            order = lambda d: "<".join(sorted(d, key=d.get))
            print(f"   {m:8s} {v:10s} " + "  ".join(f"{o} {est[o]:.3f}/{true[o]:.3f}" for o in OBS)
                  + f"   ranking est {order(est)}, true {order(true)}" + ("" if order(est) == order(true) else "  MIS-RANKED"))
    if "--figures" in sys.argv:
        import figures
        figures.make(W, tabs, observe, render, errcorr, tc, OBS)