stable_cart.linear_instability

stable_cart.linear_instability(X, X_eval, sigma=None, y=None, robust=False)[source]

Analytic prediction instability of least squares, conditional on the design.

This is the closed form of what bootstrap_instability() estimates by resampling. With \(\hat\beta \sim N(\beta, \sigma^2 (X'X)^{-1})\), the prediction at a point x has variance \(\sigma^2 x'(X'X)^{-1}x\), and every instability measure follows from it.

Parameters

X

Training design matrix of shape (n_samples, n_features). Include a column of ones if the model has an intercept; this function takes the design as given.

X_eval

Points at which to evaluate, shape (n_eval, n_features).

sigma

Standard deviation of the noise, assumed constant across observations. This is the true value, not an estimate — pass \(\hat\sigma = \sqrt{\mathrm{RSS}/(n-p)}\) for the plug-in version. Required unless robust=True.

y

Training targets. Required when robust=True, which needs residuals.

robust

Drop the constant-variance assumption and use the HC0 plug-in estimate heteroskedasticity-consistent form \((X'X)^{-1}\big(\sum_i x_i x_i' \hat e_i^2\big)(X'X)^{-1}\). This branch is an estimated asymptotic covariance, not an exact finite-sample result. The constant-variance assumption is not a technicality. Writing \(a = (X'X)^{-1}x\) and \(w_i = (x_i'a)^2\), the true prediction variance is \(\sum_i \sigma_i^2 w_i\) while the constant-variance form returns \(\bar\sigma^2 \sum_i w_i\). Their ratio is therefore the \(w\)-weighted mean of \(\sigma_i^2\) divided by its unweighted mean, so the constant-variance form is exact when the noise is uncorrelated with \(w_i\), too small when observations that move this prediction are the noisy ones, and too large when they are the quiet ones. The error has no fixed sign and no characteristic size: it is a property of the design and the noise pattern together. Under genuine homoskedasticity the two forms agree, so the cost of using this one is only the loss of a known sigma.

Returns

dict[str, Any]

variance — per-point prediction variance; s1 — per-point \(E|f_D(x)-f_{D'}(x)|\), the pairwise measure; centered_mad — per-point \(E|f_D(x)-E f_D(x)|\) under Gaussian sampling; s2 — per-point squared pairwise instability, exactly twice variance; variance_mean, s1_mean, centered_mad_mean — integrated versions.

Raises

ValueError

If inputs are nonfinite, if sigma is negative, if X and X_eval disagree on width, or if the arguments needed for the requested variance form are missing. The robust form also requires more observations than columns so residual variation can be estimated.

Notes

A RuntimeWarning is issued when X'X is so ill-conditioned that the result is not meaningful. A silently finite answer from a near-singular design is the dangerous case; an exactly singular one already raises.

s1_mean is the mean of the per-point s1, which is proportional to \(E\sqrt{v(x)}\) and not to \(\sqrt{E v(x)}\). The two differ by Jensen’s inequality whenever the variance is not constant across evaluation points, and the first is what a resampling protocol reports.

Examples

>>> import numpy as np
>>> from stable_cart import linear_instability
>>> rng = np.random.default_rng(0)
>>> X = rng.normal(size=(200, 4))
>>> out = linear_instability(X, rng.normal(size=(50, 4)), sigma=2.0)
>>> bool(np.allclose(out["s2"], 2 * out["variance"]))
True
Parameters:

  • X (NDArray[floating])

  • X_eval (NDArray[floating])

  • sigma (float | None)

  • y (NDArray[floating] | None)

  • robust (bool)

Return type:

dict[str, Any]