A Wrong Model with a Good \(r^2\)¶

def read_msa(path: Path) -> tuple[np.ndarray, np.ndarray]:
    """Read an EMSA/MAS ``.msa`` spectrum into ``(x, y)``.

    The format is a block of ``#KEYWORD : value`` header lines followed by the
    data, terminated by ``#ENDOFDATA``. This file is ``#DATATYPE : XY``, so each
    data row is a comma-separated ``x, y`` pair; skipping every ``#``-prefixed
    line handles the header and the terminator in one rule.
    """
    x: list[float] = []
    y: list[float] = []
    for line in path.read_text(encoding="latin-1").splitlines():
        if not line.strip() or line.startswith("#"):
            continue
        left, _, right = line.partition(",")
        x.append(float(left))
        y.append(float(right))
    return np.asarray(x), np.asarray(y)


def load_l3_window() -> tuple[np.ndarray, np.ndarray]:
    """Load the measured spectrum and cut to the L3 region.

    The window spans the L3 white line at 710.75 eV and the shoulder below it,
    and stops short of the L2 edge near 723 eV so the comparison stays about one
    absorption region rather than the whole spectrum. Intensities are raw
    detector counts, not normalised.
    """
    raw_x, raw_y = read_msa(SPECTRUM)
    lo, hi = WINDOW_EV
    keep = (raw_x >= lo) & (raw_x <= hi)
    return raw_x[keep], raw_y[keep]


x, y = load_l3_window()

def _peak(index: int, model_type: ModelType, amplitude: float, center: float) -> ModelNodeSpec:
    """One bounded peak component; Voigt gets its extra Lorentzian width."""
    parameters = {
        "amplitude": Parameter(value=amplitude, min=0.0),
        "center": Parameter(value=center, min=float(x.min()), max=float(x.max())),
        "sigma": Parameter(value=1.0, min=0.05, max=8.0),
    }
    if model_type is ModelType.TRUE_VOIGT:
        parameters["gamma"] = Parameter(value=0.5, min=1e-3, max=8.0)
    return ModelNodeSpec(id=f"p{index}", model_type=model_type, parameters=parameters)


def fit_ladder(x: np.ndarray, y: np.ndarray) -> dict[str, FitResult]:
    """Three models of increasing elaboration, fitted to the same window."""
    gaussian, voigt = ModelType.GAUSSIAN, ModelType.TRUE_VOIGT
    peak_height = float(y.max())
    ladder = {
        "1 Gaussian": [_peak(0, gaussian, peak_height, 710.75)],
        "2 Gaussians": [
            _peak(0, gaussian, peak_height, 710.75),
            _peak(1, gaussian, peak_height * 0.3, 708.5),
        ],
        "3 Voigts": [
            _peak(0, voigt, peak_height, 710.75),
            _peak(1, voigt, peak_height * 0.3, 708.5),
            _peak(2, voigt, peak_height * 0.2, 713.5),
        ],
    }
    data = MeasurementData(x=x.tolist(), y=y.tolist())
    return {
        name: fit(FitGraph(nodes=nodes), data, FitOptions(solver="lm"))
        for name, nodes in ladder.items()
    }


results = fit_ladder(x, y)

def runs_z(residuals: np.ndarray) -> tuple[int, float, float]:
    """Wald-Wolfowitz runs test on the residual signs.

    A *run* is a maximal stretch of same-signed residuals. Residuals that
    are independent noise around a correct model produce about
    ``2*n_pos*n_neg/n + 1`` runs; far fewer means the signs come in long
    stretches, i.e. the curve is systematically above the data in one place
    and below it in another. Returns ``(observed, expected, z)``.

    Caveat, stated because it changes how hard the number may be pushed:
    the test assumes independent residuals. This spectrum is heavily
    oversampled -- 0.05 eV channels against a 0.3 eV instrumental resolution,
    so roughly six channels per resolvable feature -- and genuinely correlated
    measurement noise would also depress the run count. Read a large negative z
    as strong corroboration of what the residual panel already shows, not as a
    calibrated p-value.
    """
    signs = np.sign(residuals)
    signs = signs[signs != 0]
    n = len(signs)
    n_pos = int((signs > 0).sum())
    n_neg = n - n_pos
    observed = int((signs[1:] != signs[:-1]).sum()) + 1
    expected = 2.0 * n_pos * n_neg / n + 1.0
    variance = (expected - 1.0) * (expected - 2.0) / (n - 1.0)
    return observed, expected, (observed - expected) / np.sqrt(variance)


def report(results: dict[str, FitResult]) -> None:
    """Print the ladder, and assert the point it is here to make."""
    print(f"{'model':13s} {'r^2':>9s} {'AIC':>9s} {'runs':>6s} {'expected':>9s} {'z':>7s}")
    stats = {}
    for name, result in results.items():
        observed, expected, z = runs_z(np.array(result.residuals))
        stats[name] = (result.r_squared, result.aic, z)
        print(
            f"{name:13s} {result.r_squared:9.5f} {result.aic:9.1f} "
            f"{observed:6d} {expected:9.1f} {z:7.1f}",
        )

    r2s = [s[0] for s in stats.values()]
    aics = [s[1] for s in stats.values()]
    zs = [s[2] for s in stats.values()]

    # Both scalar scores improve monotonically down the ladder ...
    assert r2s == sorted(r2s), "r^2 should climb at every rung"
    assert aics == sorted(aics, reverse=True), "AIC should fall at every rung"
    # ... while the residuals stay structured at every rung, including the last.
    assert max(zs) < -5.0, "every rung should remain far from random residuals"
    assert results["3 Voigts"].r_squared > 0.95, "the last rung should look excellent"


report(results)