Linear Regression Diagnostics and Metrics in C++

July 2026 · 46 min read

Five posts in this sub-series have each estimated $\beta$ a different way: the closed-form normal equations, the numerically steadier QR and SVD solvers, iterative gradient descent, penalized ridge, lasso, and elastic net, and loss-robust Huber and quantile regression. This post asks what comes after every one of those fits produces a $\hat\beta$, and it is useful to separate that question into three genuinely different layers rather than treating "how good is this model" as a single number. A predictive metric (MSE, RMSE, MAE, $R^2$) answers how close predictions land to observed targets, on data the fit has or has not seen, without any reference to how the fit was produced. A diagnostic (leverage, a studentized residual, Cook's distance, PRESS) answers a structural question about the fit itself: which rows disproportionately shaped it, and would the predictions look different with any single row removed. Inference (a standard error, a $t$-statistic, a confidence interval) answers a probabilistic question that only makes sense once a noise model is assumed: how much would $\hat\beta$ plausibly move under repeated sampling from the same data-generating process. The three layers use overlapping ingredients, the same residuals and the same $(X_b^\top X_b)^{-1}$ turn up everywhere, but they answer different questions and fail in different ways, and conflating them is the most common way a fitted model looks fine on every number that was checked and is still wrong.

Predictive Metrics: MSE, RMSE, MAE, and R-Squared

A predictive metric takes only two vectors, the observed targets $y$ and the fitted model's predictions $\hat y$, and returns a single number; it does not care whether $\hat y$ came from OLS, ridge, a Huber fit, or a lookup table. Three of the four in routine use share one definition of error, $e_i = y_i - \hat y_i$, and differ only in how they aggregate it:

MSE squares every error before averaging, so a handful of large errors dominates the number even if most predictions are close; RMSE rescales MSE back into the original units of $y$ by taking a square root, which is why RMSE rather than raw MSE is usually the number reported alongside a target's natural units; MAE averages absolute error directly and grows linearly rather than quadratically with a single bad prediction, so a dataset with a few extreme errors produces an MAE that looks much better, relative to RMSE, than a dataset whose errors are uniformly moderate. None of the three has an intrinsic scale: an RMSE of $2.0$ is excellent for a target ranging over thousands of units and useless for a target ranging over hundredths, so all three are typically read next to the target's own standard deviation or range rather than in isolation.

Train and Test R-Squared, and Why the Test Value Can Go Negative

$R^2$ rescales MSE against a fixed reference point instead of reporting it directly: $R^2 = 1 - \text{SS}_{\text{res}}/\text{SS}_{\text{tot}}$, where $\text{SS}_{\text{res}} = \sum_i(y_i-\hat y_i)^2$ and $\text{SS}_{\text{tot}} = \sum_i(y_i-\bar y)^2$ uses the mean of whichever $y$-vector is being scored. On the training set this comparison is guaranteed favorable in one direction: OLS chooses $\hat\beta$ specifically to minimize $\text{SS}_{\text{res}}$, and the constant predictor $\hat y_i = \bar y$ is one of the models it is minimizing against, so training $R^2$ is bounded below by $0$ and typically comfortably above it. On a held-out test set, no such guarantee exists. $\text{SS}_{\text{tot}}$ in the test-$R^2$ formula uses the test set's own mean, not the training mean, so the reference point the fit is compared against is a model that was never actually fit, the constant equal to the test-set average; a model whose test predictions are worse than that constant baseline scores $R^2 < 0$. This is not a corner case reserved for badly broken code: a model trained on too few rows relative to its noise level, or a model that overfits a small or unrepresentative training sample, routinely lands a negative test $R^2$ the moment the training-set noise pattern fails to describe the test set. Scenario D in the complete program below constructs exactly this case deliberately: an OLS fit on five noisy training rows with two features leaves only two residual degrees of freedom, the fitted line chases sampling noise rather than signal, and the resulting test $R^2$ prints as a negative number, not because the diagnostics code has a bug but because a genuinely poor out-of-sample fit is precisely what a negative test $R^2$ is supposed to report.

Adjusted R-Squared and Residual Standard Error

Training $R^2$ has a one-directional flaw that adjusted $R^2$ exists to correct: adding any feature, even pure noise uncorrelated with $y$, cannot increase $\text{SS}_{\text{res}}$ on the training set, since OLS is always free to set that feature's coefficient to exactly zero and any nonzero coefficient that lowers $\text{SS}_{\text{res}}$ even slightly is preferred over zero. Adjusted $R^2$ explicitly penalizes the parameter count instead of leaving that penalty implicit: $\bar R^2 = 1-(1-R^2)\dfrac{n-1}{n-d-1}$, where $d$ is the number of features (not counting the intercept). The denominator $n-d-1$ is the residual degrees of freedom after the intercept and every slope are estimated; getting this degree-of-freedom count wrong, using $n-d$ instead of $n-d-1$, or forgetting the intercept entirely, is the most common implementation bug in adjusted $R^2$, and it silently under-penalizes every additional feature. Unlike plain $R^2$, $\bar R^2$ can decrease when an added feature's improvement in fit is too small to offset the lost degree of freedom. The residual standard error reports the same residual variance on the scale of $y$ rather than as a dimensionless fraction: $s = \sqrt{\text{SSE}/(n-p)}$ with $p=d+1$ total estimated parameters and $\text{SSE}=\sum_i e_i^2$; $s$ is the quantity every standard error, $t$-statistic, and confidence interval below is built from, so getting its degrees-of-freedom count right matters beyond adjusted $R^2$ alone.

Where This Diagnostic Layer Stops Applying

Every formula from this point forward assumes an unpenalized OLS fit with an intercept, estimated by the closed-form normal equations, with approximately Gaussian, homoskedastic, independent errors. That assumption is not a formality; it is load-bearing in three separate ways, and each of the other four posts in this sub-series breaks it differently.

None of this means the other four posts' fits go undiagnosed; it means their diagnostic layer is not this one. A penalized fit's honest uncertainty estimate typically comes from resampling (the bootstrap) rather than a closed-form covariance formula that assumes zero bias; a Huber or quantile fit's own leverage-and-influence analogue exists but is built from that fit's own weighting or its own loss gradient, not from the unweighted hat matrix below. Every diagnostic in this post is scoped to the one estimator whose assumptions it was derived under.

Fitted Values and Residuals

Every quantity below is built from the same intercept-augmented design matrix $X_b=[\mathbf 1\mid X]$ and the same fitted coefficient vector $\hat\beta=(X_b^\top X_b)^{-1}X_b^\top y$ that the first post already solves for by Gaussian elimination rather than by explicit inversion. The fitted values are $\hat y = X_b\hat\beta$ and the raw residuals are $e_i=y_i-\hat y_i$; the residual sum of squares $\text{SSE}=\sum_i e_i^2$ feeds directly into MSE, RMSE, the residual standard error $s$, and every diagnostic that follows. What is new in this post is not the fit itself but $(X_b^\top X_b)^{-1}$ kept around explicitly rather than discarded once $\hat\beta$ is solved for: leverage, coefficient standard errors, and PRESS residuals all reuse that one matrix directly, so it is computed once per fit and cached on the fitted model rather than recomputed per diagnostic.

OLS Coefficient Covariance, Standard Errors, and Confidence Intervals

Under the assumptions in the scope section above, the sampling distribution of $\hat\beta$ has a closed form: $\widehat{\text{Cov}}(\hat\beta) = s^2(X_b^\top X_b)^{-1}$, where $s^2$ is the residual variance from the fitted model. The diagonal of this $p\times p$ matrix gives each coefficient's own variance; its square root is the coefficient's standard error, $\text{SE}_j = s\sqrt{[(X_b^\top X_b)^{-1}]_{jj}}$. Dividing the coefficient by its own standard error gives a $t$-statistic, $t_j=\hat\beta_j/\text{SE}_j$, testing whether that coefficient is distinguishable from zero once its own estimation uncertainty is accounted for. A normal-approximation 95% interval is then $\hat\beta_j \pm 1.95996\cdot\text{SE}_j$, using the $z_{0.975}$ quantile of the standard normal rather than a $t$-distribution quantile; this is an approximation, not an exact small-sample interval, and it is a good one exactly when the residual degrees of freedom $n-p$ are large enough that the $t$-distribution and the normal distribution are nearly indistinguishable (a common rule of thumb is $n-p \gtrsim 30$), which is why the complete program below reports it only on the $150$-row scenario and not on the eight-point hand-checkable one.

This entire machinery rests on four assumptions, and each has a direct diagnostic consequence introduced later in this post:

• Linearity in the parameters. $E[y\mid X]=X_b\beta$ for the true $\beta$; violated by an omitted nonlinear term, not by anything this post's diagnostics can detect directly.
• Homoskedasticity. $\text{Var}(e_i\mid X)=\sigma^2$ for every row, the same $\sigma^2$ that $s^2$ estimates as a single pooled number; the Breusch-Pagan check below tests this assumption directly, because every standard error above silently assumes it holds.
• Independence. errors are uncorrelated across rows; violated by, among other things, repeated measurements of the same unit or serially correlated time-ordered data.
• Approximate normality. needed for the $t$-statistic and the confidence interval to have their advertised meaning at finite $n$; less critical as $n-p$ grows, by the same large-sample reasoning that makes the normal-approximation interval reasonable in the first place.

Leverage Without Forming the Full Hat Matrix

The hat matrix $H=X_b(X_b^\top X_b)^{-1}X_b^\top$ is an $n\times n$ object, and forming it explicitly costs $O(n^2p)$ time and $O(n^2)$ memory, both wasted, because every diagnostic below needs only its diagonal, $h_{ii}=\partial\hat y_i/\partial y_i$, the amount a row's own target value is allowed to pull its own fitted value. Each diagonal entry is a single row of $X_b$ multiplied against the already-cached $p\times p$ inverse and against itself, $h_{ii}=x_{b,i}^\top(X_b^\top X_b)^{-1}x_{b,i}$, an $O(p^2)$ computation with no $n\times n$ object ever allocated; computing every row's leverage this way costs $O(np^2)$ total, linear in $n$ rather than quadratic. Three properties make leverage worth computing for every fit at that low cost:

Leverage and Cook's Distance by Hand on Eight Points

For a single feature, $h_{ii}$ reduces to a formula with no matrix inversion at all: $h_{ii}=\tfrac1n+\dfrac{(x_i-\bar x)^2}{S_{xx}}$, where $S_{xx}=\sum_i(x_i-\bar x)^2$. Take $x=[1,2,\dots,8]$ and $y=[5.1,8.0,10.7,14.2,16.8,19.5,23.9,25.6]$, close to but not exactly the line $y=2+3x$: $\bar x=4.5$, $S_{xx}=\sum_i(x_i-4.5)^2=42$, and the OLS fit is $\hat\beta_0=1.9750$, $\hat\beta_1=3.0000$ exactly (the complete program below reproduces both to four decimal places). Leverage for the last point ($x_8=8$): $h_{88}=\tfrac18+\dfrac{(8-4.5)^2}{42}=0.1250+0.2917=0.4167$, matching point $1$ by the symmetry of $(x_i-\bar x)^2$ around $\bar x$. The residual standard error over $n-p=6$ degrees of freedom is $s=0.4822$, and every remaining column below follows directly from $h_{ii}$, $e_i$, and $s$:

Point $7$'s residual, $0.9250$, is the largest in the set by a wide margin, and even though its leverage, $0.2738$, is only the third highest, the combination pushes its Cook's distance to $0.9554$, more than double the next highest value. Points $1$ and $8$ carry the highest leverage, $0.4167$ each, yet their Cook's distances, $0.0411$ and $0.3703$, stay well below point $7$'s: leverage alone does not predict influence, and the next two sections make that distinction precise.

Internally Studentized Residuals, Cook's Distance, and DFFITS

A raw residual is not directly comparable across rows: high leverage mechanically shrinks a row's own residual relative to a low-leverage row with identical underlying noise, since the fit is pulled toward exactly that row's own target. The internally studentized residual corrects for this by dividing by the residual standard error and a leverage-dependent factor, $r_i = \dfrac{e_i}{s\sqrt{1-h_{ii}}}$, putting every row on a comparable scale before it is judged as ordinary or extreme. Cook's distance combines this internally studentized residual with leverage into a single measure of how much every fitted value would change if row $i$ were dropped and the model refit: $D_i = \dfrac{r_i^2}{p}\cdot\dfrac{h_{ii}}{1-h_{ii}}$; $D_i$ is only large when a row carries both a large residual and non-trivial leverage together, since a huge residual at a low-leverage point, or large leverage with a tiny residual, both keep the product small.

DFFITS answers a related but distinct question, how many standard errors a row's own fitted value moves if that row is genuinely dropped and the model refit without it, and answering it correctly needs a residual that is actually computed against a refit-without-row-$i$ variance, not the internally studentized $r_i$ above. That is the externally studentized (or deleted) residual $t_i = \dfrac{e_i}{s_{(i)}\sqrt{1-h_{ii}}}$, where $s_{(i)}^2=\dfrac{\text{SSE}-e_i^2/(1-h_{ii})}{n-p-1}$ is the residual variance a genuine refit with row $i$ held out would produce, obtained here from an algebraic identity rather than by ever actually refitting: $\text{DFFITS}_i = t_i\sqrt{\dfrac{h_{ii}}{1-h_{ii}}}$. Substituting the internally studentized $r_i$ for $t_i$ is a common implementation mistake, not a harmless shortcut, and it is seductive precisely because it is algebraically clean: doing so makes $\text{DFFITS}_i$ identically equal to $\sqrt{p\cdot D_i}$ up to sign, since both then reduce to the same $r_i$, hiding a real distinction behind a numerical coincidence that only looks like confirmation. The gap between the two is largest exactly where DFFITS is meant to matter most, on a genuinely influential row: point $7$ from the eight-point example above has $\sqrt{p\cdot D_7}=1.3823$ under that shortcut, but $t_7$ properly excludes point $7$'s own contribution from the deleted-fit variance and gives the very different $\text{DFFITS}_7=3.2012$, as the complete program's output further below confirms. Computing $t_i$ costs nothing beyond arithmetic already on hand, since $\text{SSE}$, $e_i$, and $h_{ii}$ are all already cached from the single full fit, so no row is ever actually refit; the one case that needs an explicit check rather than a silent clip is a row so extreme that $e_i^2/(1-h_{ii})$ meets or exceeds the whole fit's $\text{SSE}$, which would make $s_{(i)}^2$ non-positive and $\text{DFFITS}_i$ undefined for that row, and a fit with $n-p-1\le 0$, which leaves no degrees of freedom for $s_{(i)}^2$ at all.

PRESS Residuals and Leave-One-Out Cross-Validation

Recomputing $n$ separate leave-one-out fits to see how well row $i$'s own target is predicted by a model that never saw it would cost $n$ extra solves of the normal equations, but an exact algebraic identity gets the same answer from quantities already on hand: the PRESS residual, the leave-one-out prediction error for row $i$ without ever refitting, is

This is not an approximation; it is the exact value a genuine leave-one-out refit would produce, a consequence of how the hat matrix updates under a single row's removal. Squaring and averaging every PRESS residual gives a leave-one-out estimate of predictive error without ever leaving the training set aside, the leave-one-out cross-validation RMSE: $\text{RMSE}_{\text{LOOCV}} = \sqrt{\tfrac1n\sum_i \text{PRESS}_i^2}$. Because PRESS divides by $1-h_{ii}$, a high-leverage row's leave-one-out error is inflated relative to its in-sample residual exactly in proportion to how much that row was allowed to pull its own fit, which is why $\text{RMSE}_{\text{LOOCV}}$ is consistently at least as large as the in-sample RMSE and is a materially more honest estimate of out-of-sample error whenever any row's leverage is far from the $p/n$ average.

The PRESS Identity, Checked Against an Explicit Refit

Point $7$ from the eight-point table above has $e_7=0.9250$ and $h_{77}=0.2738$, so $\text{PRESS}_7 = 0.9250/(1-0.2738) = 1.2738$. Refitting the same simple linear regression on the other seven points (excluding point $7$ entirely, so $\bar x$ and $S_{xx}$ change too) and predicting at $x_7=7$ with that refit line gives a residual of $1.273770$, matching $\text{PRESS}_7=1.273770$ to six decimal places with no refitting shortcut taken in the identity's derivation, only in this particular check of it. The complete program below repeats this explicit-refit comparison for every one of the eight points and again for five representative rows of the $150$-row scenario, in every case to within floating-point rounding of $10^{-14}$ or smaller, confirming the identity holds exactly rather than approximately.

Heuristic Thresholds at a Glance

Every threshold in this table is triage, not a hypothesis test with a controlled error rate; the leverage and Cook's-distance thresholds in particular are conventions that took hold because they scale sensibly with $n$ and $p$, not because they were derived from a sampling distribution. A row above one of these thresholds deserves a second look, and a row below one is not thereby proven trustworthy.

A Model-Level Diagnostic: the Breusch-Pagan Test

Every diagnostic up to this point is row-level: it says something about one observation. The Breusch-Pagan test is model-level: it asks whether the homoskedasticity assumption underlying every standard error, $t$-statistic, and confidence interval above holds across the fit as a whole, rather than whether any single row is unusual. The test regresses the squared residuals from the original fit on the same features, with their own intercept, exactly like any other OLS fit: $\hat e_i^2 = \gamma_0+\gamma_1 x_{i1}+\dots+\gamma_d x_{id}+\text{noise}$. If residual variance genuinely does not depend on $X$, this auxiliary regression should explain almost none of the variation in $e_i^2$, so its own $R^2_{\text{aux}}$ should be small. The test statistic scales that auxiliary $R^2$ by the sample size, $\text{LM}=n\cdot R^2_{\text{aux}}$, and under the null of homoskedasticity, $\text{LM}$ follows a $\chi^2_d$ distribution asymptotically; the heuristic threshold $2d$ used in the complete program below is a rough stand-in for that $\chi^2_d$ comparison rather than a substitute for a proper $p$-value, since a from-scratch, standard-library-only implementation has no chi-squared quantile function to call. A small $\text{LM}$, well under the heuristic threshold, is exactly what the synthetic homoskedastic data in the complete program below produces; a large $\text{LM}$ would indicate that the standard errors and confidence intervals from the inference section above are unreliable, understated in the region of $X$ where residual variance is actually largest.

A Clean Baseline, an Injected Influential Point, and a Refit Without It

Every formula above is easiest to trust once it is shown separating a genuinely influential row from a large body of ordinary ones at a realistic scale, not just on eight hand-picked points. The complete program below builds a clean $150$-row training set and a held-out $60$-row test set from $y=5+1.5x_1-2x_2+\varepsilon$, fits the baseline OLS model, then appends a single additional row with extreme feature values and a target that deviates sharply from the trend, refits on all $151$ rows, and finally refits again with that one row removed. Three comparisons come out of this:

The predictive-metric row of that table is the entire point: adding one row out of $151$ moves test $R^2$ by about a percentage point, from $0.9448$ to $0.9347$, a change easy to dismiss as ordinary sampling noise if it were the only number examined. The diagnostic rows tell a completely different story: the injected row's Cook's distance is $5.8764$, roughly $197\times$ the largest Cook's distance among every one of the other $150$ rows, and its DFFITS, leverage, and studentized residual (computed live in the animation below) are equally extreme. Refitting with that one row removed recovers the baseline coefficients exactly, to every printed digit, confirming that the single row was solely responsible for the shift. Predictive metrics answer "does this model still predict reasonably well"; diagnostics answer "is this fit's structure being quietly dictated by one row", and this scenario is deliberately built so both questions get different-looking answers from the same fit.

A Negative Test R-Squared, and Explicit Failure Modes

Scenario D in the complete program below fits an OLS model on only five noisy rows with two features (three degrees of freedom used, two left over) and evaluates it on thirty held-out rows drawn from the same generating process but a different noise realization. Training $R^2$ already looks unimpressive, $0.3948$, but test $R^2$ crosses zero entirely, printing $-0.8345$: the fitted line, chasing five points' worth of sampling noise rather than the underlying signal, predicts the test set worse than the test set's own mean would. This is the train/test $R^2$ discussion from the predictive-metrics section made concrete rather than hypothetical, and it is a useful sanity check in its own right: a training pipeline that never sees a negative test $R^2$ across any of its models, ever, is either evaluating exclusively on training data by mistake or has never yet been asked to fit too little data.

Scenario E demonstrates eleven explicit, intentional failure modes rather than letting them appear only as abstract possibilities, none of them caught by a broad catch-all or papered over with a clipped or default value. Fitting with more rows than targets throws immediately with a message identifying the exact row/target mismatch; fitting with fewer rows than parameters (specifically, without enough rows to leave at least one residual degree of freedom) throws before any linear algebra is attempted, naming both the required and actual row counts; and fitting two exactly collinear features (one column an exact multiple of another, engineered so $X_b^\top X_b$ is exactly singular) throws from inside solveLinearSystem itself, naming the exact pivot column where no usable pivot could be found. Beyond those three, a non-finite value anywhere in $X$ or $y$ is rejected by OLSRegressor::fit at the exact row and column responsible, rather than being allowed to propagate into $X_b^\top X_b$ and surface later as an opaque singular-matrix message; calling predict with the wrong number of features throws instead of silently reading past the end of the coefficient vector; diagnose refuses a design matrix whose row count or column count does not match the model it was fit on, and refuses a saturated or otherwise degenerate fit whose residual standard error is zero or non-finite, rather than dividing by it; and r2Score, mse, adjustedR2, and residualStandardError each refuse the input that would otherwise force a silent out-of-bounds read, an empty-vector division by zero, or an undefined value, be it a size mismatch between predictions and targets, an empty pair of vectors, a target vector with exactly zero variance, or non-positive residual degrees of freedom. Each of these is caught and its exact message printed in the byte-exact output below, so the error text is verified output, not a design claim.

Implementation

The linear algebra underneath every diagnostic in this post is two functions: partial-pivoting Gaussian elimination to solve $X_b^\top X_b\,\beta = X_b^\top y$, and a Gauss-Jordan inverse to obtain $(X_b^\top X_b)^{-1}$ explicitly, since leverage and the coefficient covariance both need that inverse directly rather than just a solved system. Both functions validate their own preconditions before touching a single pivot: a non-empty, square matrix with every entry finite, so a malformed or non-finite input is rejected with a specific message rather than producing a nonsensical inverse or an out-of-bounds access deeper in the elimination.

vector<double> solveLinearSystem(Matrix A, vector<double> b) {
    int n = (int)A.size();
    if (n == 0) throw invalid_argument("solveLinearSystem: A must be non-empty");
    if (b.size() != A.size())
        throw invalid_argument("solveLinearSystem: A and b dimensions do not match");
    double matrixScale = 0.0;
    for (int i = 0; i < n; ++i) {
        if (A[i].size() != A.size())
            throw invalid_argument("solveLinearSystem: A must be square; row " + to_string(i) +
                                    " has " + to_string(A[i].size()) + " entries, expected " + to_string(n));
        for (double v : A[i]) {
            if (!isfinite(v)) throw invalid_argument("solveLinearSystem: A contains a non-finite value");
            matrixScale = max(matrixScale, fabs(v));
        }
        if (!isfinite(b[i])) throw invalid_argument("solveLinearSystem: b contains a non-finite value");
    }
    const double tol = numeric_limits<double>::epsilon() * max(1.0, matrixScale) * n * 100.0;
    for (int k = 0; k < n; ++k) {
        int piv = k;
        double best = fabs(A[k][k]);
        for (int i = k + 1; i < n; ++i)
            if (fabs(A[i][k]) > best) { best = fabs(A[i][k]); piv = i; }
        if (best <= tol)
            throw runtime_error("solveLinearSystem: matrix is singular or near-singular "
                                "(no usable pivot in column " + to_string(k) + "); "
                                "check for duplicated/collinear features");
        if (piv != k) { swap(A[k], A[piv]); swap(b[k], b[piv]); }
        for (int i = k + 1; i < n; ++i) {
            double factor = A[i][k] / A[k][k];
            if (factor == 0.0) continue;
            for (int j = k; j < n; ++j) A[i][j] -= factor * A[k][j];
            b[i] -= factor * b[k];
        }
    }
    vector<double> beta(n, 0.0);
    for (int k = n - 1; k >= 0; --k) {
        double s = b[k];
        for (int j = k + 1; j < n; ++j) s -= A[k][j] * beta[j];
        beta[k] = s / A[k][k];
    }
    return beta;
}

Matrix invertMatrix(Matrix A) {
    // Same preconditions as solveLinearSystem, since XtXinv is built from the same
    // XtX matrix: non-empty, square, and every entry finite before elimination starts.
    int n = (int)A.size();
    if (n == 0) throw invalid_argument("invertMatrix: A must be non-empty");
    double matrixScale = 0.0;
    for (int i = 0; i < n; ++i) {
        if (A[i].size() != A.size())
            throw invalid_argument("invertMatrix: A must be square; row " + to_string(i) +
                                    " has " + to_string(A[i].size()) + " entries, expected " + to_string(n));
        for (double v : A[i]) {
            if (!isfinite(v)) throw invalid_argument("invertMatrix: A contains a non-finite value");
            matrixScale = max(matrixScale, fabs(v));
        }
    }
    Matrix I(n, vector<double>(n, 0.0));
    for (int i = 0; i < n; ++i) I[i][i] = 1.0;
    const double tol = numeric_limits<double>::epsilon() * max(1.0, matrixScale) * n * 100.0;
    for (int k = 0; k < n; ++k) {
        int piv = k;
        double best = fabs(A[k][k]);
        for (int i = k + 1; i < n; ++i)
            if (fabs(A[i][k]) > best) { best = fabs(A[i][k]); piv = i; }
        if (best <= tol)
            throw runtime_error("invertMatrix: matrix is singular or near-singular "
                                 "(no usable pivot in column " + to_string(k) + ")");
        if (piv != k) { swap(A[k], A[piv]); swap(I[k], I[piv]); }
        double d = A[k][k];
        for (int j = 0; j < n; ++j) { A[k][j] /= d; I[k][j] /= d; }
        for (int i = 0; i < n; ++i) {
            if (i == k) continue;
            double factor = A[i][k];
            if (factor == 0.0) continue;
            for (int j = 0; j < n; ++j) { A[i][j] -= factor * A[k][j]; I[i][j] -= factor * I[k][j]; }
        }
    }
    return I;
}

OLSRegressor keeps $(X_b^\top X_b)^{-1}$ and the residual standard error $s$ on the fitted object, since both are reused by every diagnostic below rather than recomputed per call. fit validates every row of $X$ and every entry of $y$ for a matching shape and a finite value before any arithmetic touches them, and requires strictly more rows than parameters so that $n-p$, the residual degrees of freedom $s$ is built from, is always positive rather than silently falling back to an $s$ of zero; predict checks that a query row has exactly $d$ features before reading past whichever end of the coefficient vector a shorter or longer row would otherwise walk off:

struct OLSRegressor {
    bool fitted = false;
    int n = 0, d = 0, p = 0;
    double intercept = 0.0;
    vector<double> coef;
    Matrix XtXinv;      // (X_b^T X_b)^{-1}: leverage and coefficient covariance both reuse this
    double sse = 0.0;
    double s = 0.0;      // residual standard error, sqrt(sse / (n - p))

    void fit(const Matrix& X, const vector<double>& y) {
        if (X.size() != y.size())
            throw invalid_argument("OLSRegressor::fit: X has " + to_string(X.size()) +
                                    " rows but y has " + to_string(y.size()) + " entries");
        n = (int)X.size();
        if (n == 0) throw invalid_argument("OLSRegressor::fit: empty training set");
        d = (int)X[0].size();
        if (d == 0) throw invalid_argument("OLSRegressor::fit: zero feature columns");
        for (int i = 0; i < n; ++i) {
            if ((int)X[i].size() != d)
                throw invalid_argument("OLSRegressor::fit: row " + to_string(i) + " has " +
                                        to_string(X[i].size()) + " columns, expected " + to_string(d));
            // Reject non-finite X/y up front, at the exact row/column that caused it,
            // rather than letting a NaN or Inf propagate into X_b^T X_b and surface
            // later as an opaque "matrix contains a non-finite value".
            for (int j = 0; j < d; ++j)
                if (!isfinite(X[i][j]))
                    throw invalid_argument("OLSRegressor::fit: X[" + to_string(i) + "][" + to_string(j) +
                                            "] is not finite");
            if (!isfinite(y[i]))
                throw invalid_argument("OLSRegressor::fit: y[" + to_string(i) + "] is not finite");
        }
        p = d + 1;
        // n > p (not just n >= p) is required: with n == p the fit is saturated and
        // has zero residual degrees of freedom, so s below would have no honest value
        // to fall back on other than a silently wrong 0.
        if (n <= p)
            throw invalid_argument("OLSRegressor::fit: need at least " + to_string(p + 1) +
                                    " rows for " + to_string(d) + " features plus intercept " +
                                    "(residual degrees of freedom must be positive), got " + to_string(n));

        Matrix Xb = augment(X);
        Matrix XtX(p, vector<double>(p, 0.0));
        vector<double> Xty(p, 0.0);
        for (int a = 0; a < p; ++a) {
            for (int i = 0; i < n; ++i) Xty[a] += Xb[i][a] * y[i];
            for (int b = 0; b < p; ++b) {
                double acc = 0.0;
                for (int i = 0; i < n; ++i) acc += Xb[i][a] * Xb[i][b];
                XtX[a][b] = acc;
            }
        }
        vector<double> beta = solveLinearSystem(XtX, Xty);
        XtXinv = invertMatrix(XtX);
        intercept = beta[0];
        coef.assign(beta.begin() + 1, beta.end());
        fitted = true;

        sse = 0.0;
        for (int i = 0; i < n; ++i) { double e = y[i] - predict(X[i]); sse += e * e; }
        int df = n - p;
        s = sqrt(sse / df);
    }

    double predict(const vector<double>& x) const {
        if (!fitted) throw logic_error("OLSRegressor::predict: model is not fitted yet");
        if ((int)x.size() != d)
            throw invalid_argument("OLSRegressor::predict: expected " + to_string(d) +
                                    " features, got " + to_string(x.size()));
        double out = intercept;
        for (int j = 0; j < d; ++j) {
            if (!isfinite(x[j]))
                throw invalid_argument("OLSRegressor::predict: feature " + to_string(j) +
                                       " is not finite");
            out += coef[j] * x[j];
        }
        if (!isfinite(out))
            throw runtime_error("OLSRegressor::predict: prediction became non-finite");
        return out;
    }
    vector<double> predict(const Matrix& X) const {
        vector<double> out; out.reserve(X.size());
        for (auto& row : X) out.push_back(predict(row));
        return out;
    }
};

Predictive metrics are free functions of $(y,\hat y)$ alone, unaware of how $\hat y$ was produced. Each validates that the two vectors it is given agree in length, are non-empty, and contain only finite values before indexing or accumulating them, since a silent length mismatch is an out-of-bounds read and a single NaN would silently poison every aggregate. r2Score additionally refuses a constant $y$, whose $\text{SS}_{\text{tot}}$ is exactly zero and would otherwise make $R^2$ a division by zero rather than a well-defined number, and adjustedR2 and residualStandardError both refuse non-positive residual degrees of freedom rather than returning a value computed from a division by zero or a negative count:

// Shared precondition check: every metric below needs the same two vectors to
// agree in length and be non-empty, so validate both once rather than repeating
// an ad hoc check (or, worse, an out-of-bounds read) in each function.
void checkMetricInputs(const char* who, const vector<double>& yTrue, const vector<double>& yPred) {
    if (yTrue.size() != yPred.size())
        throw invalid_argument(string(who) + ": yTrue has " + to_string(yTrue.size()) +
                                " entries but yPred has " + to_string(yPred.size()));
    if (yTrue.empty()) throw invalid_argument(string(who) + ": yTrue and yPred must be non-empty");
    for (size_t i = 0; i < yTrue.size(); ++i)
        if (!isfinite(yTrue[i]) || !isfinite(yPred[i]))
            throw invalid_argument(string(who) + ": non-finite value at index " +
                                   to_string(i));
}
double mse(const vector<double>& yTrue, const vector<double>& yPred) {
    checkMetricInputs("mse", yTrue, yPred);
    double s = 0.0;
    for (size_t i = 0; i < yTrue.size(); ++i) { double e = yTrue[i] - yPred[i]; s += e * e; }
    return s / yTrue.size();
}
double rmseOf(const vector<double>& yTrue, const vector<double>& yPred) { return sqrt(mse(yTrue, yPred)); }
double maeOf(const vector<double>& yTrue, const vector<double>& yPred) {
    checkMetricInputs("maeOf", yTrue, yPred);
    double s = 0.0;
    for (size_t i = 0; i < yTrue.size(); ++i) s += fabs(yTrue[i] - yPred[i]);
    return s / yTrue.size();
}
// R^2 against the mean of whichever y-vector is passed in, so a test-set call
// compares against the test-set mean and can legitimately return a negative value.
double r2Score(const vector<double>& yTrue, const vector<double>& yPred) {
    checkMetricInputs("r2Score", yTrue, yPred);
    double m = meanOf(yTrue), ssRes = 0.0, ssTot = 0.0;
    for (size_t i = 0; i < yTrue.size(); ++i) {
        ssRes += (yTrue[i] - yPred[i]) * (yTrue[i] - yPred[i]);
        ssTot += (yTrue[i] - m) * (yTrue[i] - m);
    }
    // A constant yTrue makes SS_tot exactly zero, so R^2 has no defined value
    // (not even "1.0"); refuse instead of returning the -inf/NaN division would produce.
    if (!(ssTot > 0.0))
        throw invalid_argument("r2Score: yTrue is constant (SS_tot = 0); R^2 is undefined");
    return 1.0 - ssRes / ssTot;
}
double adjustedR2(double r2v, int n, int d) {
    if (n - d - 1 <= 0)
        throw invalid_argument("adjustedR2: need n - d - 1 > 0 (n=" + to_string(n) + ", d=" +
                                to_string(d) + "); residual degrees of freedom must be positive");
    return 1.0 - (1.0 - r2v) * (double)(n - 1) / (double)(n - d - 1);
}
double residualStandardError(double sseIn, int n, int p) {
    if (n - p <= 0)
        throw invalid_argument("residualStandardError: need n - p > 0 (n=" + to_string(n) + ", p=" +
                                to_string(p) + "); residual degrees of freedom must be positive");
    return sqrt(sseIn / (double)(n - p));
}

Row-level diagnostics form one struct, computed in a single pass that never allocates an $n\times n$ object. diagnose first checks that $X$, yTrue, and the fitted model agree on row count, column count, and that the model has actually been fitted, since every formula below is only meaningful when $X$ is the exact design matrix $(X_b^\top X_b)^{-1}$ was built from; it then checks that the model's residual standard error is positive and finite, since a saturated or otherwise degenerate fit has no defined studentized residual to compute:

struct Diagnostics {
    vector<double> h;         // leverage, hat-matrix diagonal only
    vector<double> resid;     // raw residuals
    vector<double> stdResid;  // internally studentized residuals
    vector<double> cooksD;    // Cook's distance
    vector<double> dffits;    // DFFITS
    vector<double> press;     // PRESS (leave-one-out) residuals
};

Diagnostics diagnose(const Matrix& X, const vector<double>& yTrue, const OLSRegressor& fit) {
    // Every quantity below assumes X, yTrue, and fit describe the exact same fitted
    // data: h_ii is only "leverage" (and PRESS only "leave-one-out") when the row
    // actually took part in the fit whose (X_b^T X_b)^{-1} is being reused here.
    if (!fit.fitted) throw invalid_argument("diagnose: model is not fitted yet");
    int n = (int)X.size();
    if (n != (int)yTrue.size())
        throw invalid_argument("diagnose: X has " + to_string(n) + " rows but yTrue has " +
                                to_string(yTrue.size()) + " entries");
    if (n != fit.n)
        throw invalid_argument("diagnose: X has " + to_string(n) + " rows but the fitted model was " +
                                "trained on " + to_string(fit.n) + " rows; diagnose requires the exact " +
                                "training data the model was fit on");
    for (int i = 0; i < n; ++i)
        if ((int)X[i].size() != fit.d)
            throw invalid_argument("diagnose: row " + to_string(i) + " has " + to_string(X[i].size()) +
                                    " columns, expected " + to_string(fit.d) + " to match the fitted model");
    int p = fit.p;
    // DFFITS needs a residual-variance estimate with row i deleted, which itself
    // needs n - p - 1 residual degrees of freedom; refuse rather than silently
    // producing a nonsensical or infinite value.
    if (n - p - 1 <= 0)
        throw invalid_argument("diagnose: DFFITS requires n - p - 1 > 0 (need at least " +
                                to_string(p + 2) + " rows for p=" + to_string(p) + ", got " +
                                to_string(n) + ")");
    // A saturated or otherwise perfect fit has s == 0 (or, degenerately, non-finite);
    // every studentized residual below divides by s, so refuse rather than emitting
    // a silent NaN or Inf for a residual variance that is not actually estimable.
    if (!isfinite(fit.s) || !(fit.s > 0.0))
        throw runtime_error("diagnose: residual standard error must be positive and finite, got " +
                             to_string(fit.s) + "; studentized residuals are undefined for a " +
                             "saturated or degenerate fit");
    Matrix Xb = augment(X);
    Diagnostics diag;
    diag.h.resize(n);
    diag.resid.resize(n);
    for (int i = 0; i < n; ++i) {
        // h_ii = x_bi^T (X_b^T X_b)^{-1} x_bi: one length-p row against a
        // p x p inverse, so the full n x n hat matrix is never materialized.
        double hii = 0.0;
        for (int a = 0; a < p; ++a) {
            double acc = 0.0;
            for (int b = 0; b < p; ++b) acc += fit.XtXinv[a][b] * Xb[i][b];
            hii += Xb[i][a] * acc;
        }
        diag.h[i] = hii;
        diag.resid[i] = yTrue[i] - fit.predict(X[i]);
    }
    diag.stdResid.resize(n);
    diag.cooksD.resize(n);
    diag.dffits.resize(n);
    diag.press.resize(n);
    for (int i = 0; i < n; ++i) {
        double oneMinusH = 1.0 - diag.h[i];
        if (!(oneMinusH > 0.0))
            throw runtime_error("diagnose: leverage must be below 1 at row " +
                                to_string(i) + ", got " + to_string(diag.h[i]));
        // Internally studentized residual r_i: divides by the one fit's own s,
        // computed with row i still included. Cook's distance uses this one.
        diag.stdResid[i] = diag.resid[i] / (fit.s * sqrt(oneMinusH));
        diag.cooksD[i] = (diag.stdResid[i] * diag.stdResid[i] / p) * (diag.h[i] / oneMinusH);
        // DFFITS needs the externally studentized (deleted) residual t_i instead:
        // s_(i)^2 = (SSE - e_i^2/(1-h_ii)) / (n-p-1) is the residual variance a
        // genuine refit without row i would produce, obtained here from an
        // algebraic identity rather than by actually refitting. A row extreme
        // enough that its own squared residual, inflated by 1/(1-h_ii), meets or
        // exceeds the whole fit's SSE makes this deleted variance non-positive;
        // that is a genuinely undefined DFFITS, not a value to clip or fall back on.
        double deletedVar = (fit.sse - diag.resid[i] * diag.resid[i] / oneMinusH) / (double)(n - p - 1);
        if (!(deletedVar > 0.0))
            throw runtime_error("diagnose: non-positive deleted residual variance at row " +
                                 to_string(i) + " (s_(i)^2=" + to_string(deletedVar) +
                                 "); DFFITS is undefined for this row");
        double sDeleted = sqrt(deletedVar);
        double tExt = diag.resid[i] / (sDeleted * sqrt(oneMinusH));
        diag.dffits[i] = tExt * sqrt(diag.h[i] / oneMinusH);
        diag.press[i] = diag.resid[i] / oneMinusH;
    }
    return diag;
}

The Breusch-Pagan check is a single free function that reuses OLSRegressor itself on the squared residuals as its auxiliary regression:

// LM = n * R^2_aux from regressing squared residuals on the same features.
double breuschPaganLM(const Matrix& X, const vector<double>& resid) {
    int n = (int)X.size();
    vector<double> e2(n);
    for (int i = 0; i < n; ++i) e2[i] = resid[i] * resid[i];
    OLSRegressor aux;
    aux.fit(X, e2);
    double r2aux = r2Score(e2, aux.predict(X));
    return n * r2aux;
}

Complete Program and Output

The program below is self-contained: solveLinearSystem, invertMatrix, OLSRegressor, every predictive-metric function, Diagnostics/diagnose, and breuschPaganLM exactly as shown above, followed by five scenarios in main. Scenario A is the eight-point hand-checked example, including the explicit leave-one-out refit check against every PRESS residual. Scenario B fits the $150$-row/$60$-row train/test split, reports every predictive metric, the OLS inference table, the Breusch-Pagan check, and the LOOCV RMSE with its own explicit-refit agreement check on five representative rows. Scenario C injects the single high-leverage, high-residual row, refits, reports the clean-versus-injected comparison, and refits again with the row removed. Scenario D is the five-row negative-test-$R^2$ demonstration. Scenario E triggers and catches eleven explicit failure modes, one per invalid input each function above was hardened against.

#include <algorithm>
#include <cmath>
#include <cstdio>
#include <limits>
#include <numeric>
#include <random>
#include <stdexcept>
#include <string>
#include <vector>
using namespace std;

using Matrix = vector<vector<double>>;

vector<double> solveLinearSystem(Matrix A, vector<double> b) {
    int n = (int)A.size();
    if (n == 0) throw invalid_argument("solveLinearSystem: A must be non-empty");
    if (b.size() != A.size())
        throw invalid_argument("solveLinearSystem: A and b dimensions do not match");
    double matrixScale = 0.0;
    for (int i = 0; i < n; ++i) {
        if (A[i].size() != A.size())
            throw invalid_argument("solveLinearSystem: A must be square; row " + to_string(i) +
                                    " has " + to_string(A[i].size()) + " entries, expected " + to_string(n));
        for (double v : A[i]) {
            if (!isfinite(v)) throw invalid_argument("solveLinearSystem: A contains a non-finite value");
            matrixScale = max(matrixScale, fabs(v));
        }
        if (!isfinite(b[i])) throw invalid_argument("solveLinearSystem: b contains a non-finite value");
    }
    const double tol = numeric_limits<double>::epsilon() * max(1.0, matrixScale) * n * 100.0;
    for (int k = 0; k < n; ++k) {
        int piv = k;
        double best = fabs(A[k][k]);
        for (int i = k + 1; i < n; ++i)
            if (fabs(A[i][k]) > best) { best = fabs(A[i][k]); piv = i; }
        if (best <= tol)
            throw runtime_error("solveLinearSystem: matrix is singular or near-singular "
                                "(no usable pivot in column " + to_string(k) + "); "
                                "check for duplicated/collinear features");
        if (piv != k) { swap(A[k], A[piv]); swap(b[k], b[piv]); }
        for (int i = k + 1; i < n; ++i) {
            double factor = A[i][k] / A[k][k];
            if (factor == 0.0) continue;
            for (int j = k; j < n; ++j) A[i][j] -= factor * A[k][j];
            b[i] -= factor * b[k];
        }
    }
    vector<double> beta(n, 0.0);
    for (int k = n - 1; k >= 0; --k) {
        double s = b[k];
        for (int j = k + 1; j < n; ++j) s -= A[k][j] * beta[j];
        beta[k] = s / A[k][k];
    }
    return beta;
}

Matrix invertMatrix(Matrix A) {
    int n = (int)A.size();
    if (n == 0) throw invalid_argument("invertMatrix: A must be non-empty");
    double matrixScale = 0.0;
    for (int i = 0; i < n; ++i) {
        if (A[i].size() != A.size())
            throw invalid_argument("invertMatrix: A must be square; row " + to_string(i) +
                                    " has " + to_string(A[i].size()) + " entries, expected " + to_string(n));
        for (double v : A[i]) {
            if (!isfinite(v)) throw invalid_argument("invertMatrix: A contains a non-finite value");
            matrixScale = max(matrixScale, fabs(v));
        }
    }
    Matrix I(n, vector<double>(n, 0.0));
    for (int i = 0; i < n; ++i) I[i][i] = 1.0;
    const double tol = numeric_limits<double>::epsilon() * max(1.0, matrixScale) * n * 100.0;
    for (int k = 0; k < n; ++k) {
        int piv = k;
        double best = fabs(A[k][k]);
        for (int i = k + 1; i < n; ++i)
            if (fabs(A[i][k]) > best) { best = fabs(A[i][k]); piv = i; }
        if (best <= tol)
            throw runtime_error("invertMatrix: matrix is singular or near-singular "
                                 "(no usable pivot in column " + to_string(k) + ")");
        if (piv != k) { swap(A[k], A[piv]); swap(I[k], I[piv]); }
        double d = A[k][k];
        for (int j = 0; j < n; ++j) { A[k][j] /= d; I[k][j] /= d; }
        for (int i = 0; i < n; ++i) {
            if (i == k) continue;
            double factor = A[i][k];
            if (factor == 0.0) continue;
            for (int j = 0; j < n; ++j) { A[i][j] -= factor * A[k][j]; I[i][j] -= factor * I[k][j]; }
        }
    }
    return I;
}

Matrix augment(const Matrix& X) {
    int n = (int)X.size(), d = (int)X[0].size();
    Matrix Xb(n, vector<double>(d + 1));
    for (int i = 0; i < n; ++i) {
        Xb[i][0] = 1.0;
        for (int j = 0; j < d; ++j) Xb[i][j + 1] = X[i][j];
    }
    return Xb;
}

struct OLSRegressor {
    bool fitted = false;
    int n = 0, d = 0, p = 0;
    double intercept = 0.0;
    vector<double> coef;
    Matrix XtXinv;
    double sse = 0.0;
    double s = 0.0;

    void fit(const Matrix& X, const vector<double>& y) {
        if (X.size() != y.size())
            throw invalid_argument("OLSRegressor::fit: X has " + to_string(X.size()) +
                                    " rows but y has " + to_string(y.size()) + " entries");
        n = (int)X.size();
        if (n == 0) throw invalid_argument("OLSRegressor::fit: empty training set");
        d = (int)X[0].size();
        if (d == 0) throw invalid_argument("OLSRegressor::fit: zero feature columns");
        for (int i = 0; i < n; ++i) {
            if ((int)X[i].size() != d)
                throw invalid_argument("OLSRegressor::fit: row " + to_string(i) + " has " +
                                        to_string(X[i].size()) + " columns, expected " + to_string(d));
            for (int j = 0; j < d; ++j)
                if (!isfinite(X[i][j]))
                    throw invalid_argument("OLSRegressor::fit: X[" + to_string(i) + "][" + to_string(j) +
                                            "] is not finite");
            if (!isfinite(y[i]))
                throw invalid_argument("OLSRegressor::fit: y[" + to_string(i) + "] is not finite");
        }
        p = d + 1;
        if (n <= p)
            throw invalid_argument("OLSRegressor::fit: need at least " + to_string(p + 1) +
                                    " rows for " + to_string(d) + " features plus intercept " +
                                    "(residual degrees of freedom must be positive), got " + to_string(n));

        Matrix Xb = augment(X);
        Matrix XtX(p, vector<double>(p, 0.0));
        vector<double> Xty(p, 0.0);
        for (int a = 0; a < p; ++a) {
            for (int i = 0; i < n; ++i) Xty[a] += Xb[i][a] * y[i];
            for (int b = 0; b < p; ++b) {
                double acc = 0.0;
                for (int i = 0; i < n; ++i) acc += Xb[i][a] * Xb[i][b];
                XtX[a][b] = acc;
            }
        }
        vector<double> beta = solveLinearSystem(XtX, Xty);
        XtXinv = invertMatrix(XtX);
        intercept = beta[0];
        coef.assign(beta.begin() + 1, beta.end());
        fitted = true;

        sse = 0.0;
        for (int i = 0; i < n; ++i) { double e = y[i] - predict(X[i]); sse += e * e; }
        int df = n - p;
        s = sqrt(sse / df);
    }

    double predict(const vector<double>& x) const {
        if (!fitted) throw logic_error("OLSRegressor::predict: model is not fitted yet");
        if ((int)x.size() != d)
            throw invalid_argument("OLSRegressor::predict: expected " + to_string(d) +
                                    " features, got " + to_string(x.size()));
        double out = intercept;
        for (int j = 0; j < d; ++j) {
            if (!isfinite(x[j]))
                throw invalid_argument("OLSRegressor::predict: feature " + to_string(j) +
                                       " is not finite");
            out += coef[j] * x[j];
        }
        if (!isfinite(out))
            throw runtime_error("OLSRegressor::predict: prediction became non-finite");
        return out;
    }
    vector<double> predict(const Matrix& X) const {
        vector<double> out; out.reserve(X.size());
        for (auto& row : X) out.push_back(predict(row));
        return out;
    }
};

void checkMetricInputs(const char* who, const vector<double>& yTrue, const vector<double>& yPred) {
    if (yTrue.size() != yPred.size())
        throw invalid_argument(string(who) + ": yTrue has " + to_string(yTrue.size()) +
                                " entries but yPred has " + to_string(yPred.size()));
    if (yTrue.empty()) throw invalid_argument(string(who) + ": yTrue and yPred must be non-empty");
    for (size_t i = 0; i < yTrue.size(); ++i)
        if (!isfinite(yTrue[i]) || !isfinite(yPred[i]))
            throw invalid_argument(string(who) + ": non-finite value at index " +
                                   to_string(i));
}
double meanOf(const vector<double>& v) {
    if (v.empty()) throw invalid_argument("meanOf: v must be non-empty");
    return accumulate(v.begin(), v.end(), 0.0) / v.size();
}
double mse(const vector<double>& yTrue, const vector<double>& yPred) {
    checkMetricInputs("mse", yTrue, yPred);
    double s = 0.0;
    for (size_t i = 0; i < yTrue.size(); ++i) { double e = yTrue[i] - yPred[i]; s += e * e; }
    return s / yTrue.size();
}
double rmseOf(const vector<double>& yTrue, const vector<double>& yPred) { return sqrt(mse(yTrue, yPred)); }
double maeOf(const vector<double>& yTrue, const vector<double>& yPred) {
    checkMetricInputs("maeOf", yTrue, yPred);
    double s = 0.0;
    for (size_t i = 0; i < yTrue.size(); ++i) s += fabs(yTrue[i] - yPred[i]);
    return s / yTrue.size();
}
double r2Score(const vector<double>& yTrue, const vector<double>& yPred) {
    checkMetricInputs("r2Score", yTrue, yPred);
    double m = meanOf(yTrue), ssRes = 0.0, ssTot = 0.0;
    for (size_t i = 0; i < yTrue.size(); ++i) {
        ssRes += (yTrue[i] - yPred[i]) * (yTrue[i] - yPred[i]);
        ssTot += (yTrue[i] - m) * (yTrue[i] - m);
    }
    if (!(ssTot > 0.0))
        throw invalid_argument("r2Score: yTrue is constant (SS_tot = 0); R^2 is undefined");
    return 1.0 - ssRes / ssTot;
}
double adjustedR2(double r2v, int n, int d) {
    if (n - d - 1 <= 0)
        throw invalid_argument("adjustedR2: need n - d - 1 > 0 (n=" + to_string(n) + ", d=" +
                                to_string(d) + "); residual degrees of freedom must be positive");
    return 1.0 - (1.0 - r2v) * (double)(n - 1) / (double)(n - d - 1);
}
double residualStandardError(double sseIn, int n, int p) {
    if (n - p <= 0)
        throw invalid_argument("residualStandardError: need n - p > 0 (n=" + to_string(n) + ", p=" +
                                to_string(p) + "); residual degrees of freedom must be positive");
    return sqrt(sseIn / (double)(n - p));
}

struct Diagnostics {
    vector<double> h, resid, stdResid, cooksD, dffits, press;
};

Diagnostics diagnose(const Matrix& X, const vector<double>& yTrue, const OLSRegressor& fit) {
    if (!fit.fitted) throw invalid_argument("diagnose: model is not fitted yet");
    int n = (int)X.size();
    if (n != (int)yTrue.size())
        throw invalid_argument("diagnose: X has " + to_string(n) + " rows but yTrue has " +
                                to_string(yTrue.size()) + " entries");
    if (n != fit.n)
        throw invalid_argument("diagnose: X has " + to_string(n) + " rows but the fitted model was " +
                                "trained on " + to_string(fit.n) + " rows; diagnose requires the exact " +
                                "training data the model was fit on");
    for (int i = 0; i < n; ++i)
        if ((int)X[i].size() != fit.d)
            throw invalid_argument("diagnose: row " + to_string(i) + " has " + to_string(X[i].size()) +
                                    " columns, expected " + to_string(fit.d) + " to match the fitted model");
    int p = fit.p;
    if (n - p - 1 <= 0)
        throw invalid_argument("diagnose: DFFITS requires n - p - 1 > 0 (need at least " +
                                to_string(p + 2) + " rows for p=" + to_string(p) + ", got " +
                                to_string(n) + ")");
    if (!isfinite(fit.s) || !(fit.s > 0.0))
        throw runtime_error("diagnose: residual standard error must be positive and finite, got " +
                             to_string(fit.s) + "; studentized residuals are undefined for a " +
                             "saturated or degenerate fit");
    Matrix Xb = augment(X);
    Diagnostics diag;
    diag.h.resize(n);
    diag.resid.resize(n);
    for (int i = 0; i < n; ++i) {
        double hii = 0.0;
        for (int a = 0; a < p; ++a) {
            double acc = 0.0;
            for (int b = 0; b < p; ++b) acc += fit.XtXinv[a][b] * Xb[i][b];
            hii += Xb[i][a] * acc;
        }
        diag.h[i] = hii;
        diag.resid[i] = yTrue[i] - fit.predict(X[i]);
    }
    diag.stdResid.resize(n);
    diag.cooksD.resize(n);
    diag.dffits.resize(n);
    diag.press.resize(n);
    for (int i = 0; i < n; ++i) {
        double oneMinusH = 1.0 - diag.h[i];
        if (!(oneMinusH > 0.0))
            throw runtime_error("diagnose: leverage must be below 1 at row " +
                                to_string(i) + ", got " + to_string(diag.h[i]));
        diag.stdResid[i] = diag.resid[i] / (fit.s * sqrt(oneMinusH));
        diag.cooksD[i] = (diag.stdResid[i] * diag.stdResid[i] / p) * (diag.h[i] / oneMinusH);
        double deletedVar = (fit.sse - diag.resid[i] * diag.resid[i] / oneMinusH) / (double)(n - p - 1);
        if (!(deletedVar > 0.0))
            throw runtime_error("diagnose: non-positive deleted residual variance at row " +
                                 to_string(i) + " (s_(i)^2=" + to_string(deletedVar) +
                                 "); DFFITS is undefined for this row");
        double sDeleted = sqrt(deletedVar);
        double tExt = diag.resid[i] / (sDeleted * sqrt(oneMinusH));
        diag.dffits[i] = tExt * sqrt(diag.h[i] / oneMinusH);
        diag.press[i] = diag.resid[i] / oneMinusH;
    }
    return diag;
}

double breuschPaganLM(const Matrix& X, const vector<double>& resid) {
    int n = (int)X.size();
    vector<double> e2(n);
    for (int i = 0; i < n; ++i) e2[i] = resid[i] * resid[i];
    OLSRegressor aux;
    aux.fit(X, e2);
    double r2aux = r2Score(e2, aux.predict(X));
    return n * r2aux;
}

int main() {
    printf("=== Scenario A: predictive metrics, leverage, Cook's D, and PRESS on eight hand-checkable points ===\n");
    {
        vector<double> x = {1, 2, 3, 4, 5, 6, 7, 8};
        vector<double> y = {5.1, 8.0, 10.7, 14.2, 16.8, 19.5, 23.9, 25.6};
        int n = (int)x.size();
        Matrix X;
        for (double xi : x) X.push_back({xi});

        OLSRegressor ols;
        ols.fit(X, y);
        auto pred = ols.predict(X);
        printf("OLS fit: intercept = %.4f, coef = %.4f\n", ols.intercept, ols.coef[0]);
        printf("MSE = %.4f  RMSE = %.4f  MAE = %.4f  R2 = %.4f  adj.R2 = %.4f  RSE = %.4f (df=%d)\n",
               mse(y, pred), rmseOf(y, pred), maeOf(y, pred), r2Score(y, pred),
               adjustedR2(r2Score(y, pred), n, 1), residualStandardError(ols.sse, n, ols.p), n - ols.p);

        Diagnostics diag = diagnose(X, y, ols);
        double sumH = 0.0;
        for (double v : diag.h) sumH += v;
        printf("%-6s %6s %10s %10s %10s %10s %10s %10s\n",
               "point", "x", "resid", "leverage", "std.res", "Cook's D", "DFFITS", "PRESS");
        for (int i = 0; i < n; ++i) {
            printf("%-6d %6.1f %10.4f %10.4f %10.4f %10.4f %10.4f %10.4f\n",
                   i + 1, x[i], diag.resid[i], diag.h[i], diag.stdResid[i], diag.cooksD[i],
                   diag.dffits[i], diag.press[i]);
        }
        printf("numerical check: sum(h_ii) = %.6f, expected p = %d\n", sumH, ols.p);

        double pressSum = 0.0;
        for (double pr : diag.press) pressSum += pr * pr;
        double loocvRmse = sqrt(pressSum / n);
        printf("LOOCV RMSE (via PRESS identity) = %.4f (vs in-sample RMSE %.4f)\n", loocvRmse, rmseOf(y, pred));

        printf("explicit leave-one-out agreement check (refit without row i, predict row i):\n");
        double maxDiff = 0.0;
        for (int skip = 0; skip < n; ++skip) {
            Matrix Xloo;
            vector<double> yloo;
            for (int i = 0; i < n; ++i) {
                if (i == skip) continue;
                Xloo.push_back(X[i]);
                yloo.push_back(y[i]);
            }
            OLSRegressor loo;
            loo.fit(Xloo, yloo);
            double predSkip = loo.predict(X[skip]);
            double explicitResid = y[skip] - predSkip;
            double diff = fabs(explicitResid - diag.press[skip]);
            maxDiff = max(maxDiff, diff);
            printf("  row %d: PRESS=%.6f explicit-LOO-residual=%.6f |diff|=%.2e\n",
                   skip + 1, diag.press[skip], explicitResid, diff);
        }
        printf("max |PRESS - explicit LOO residual| over all %d rows = %.2e\n", n, maxDiff);

        bool allFinite = true;
        for (int i = 0; i < n; ++i)
            if (!isfinite(diag.cooksD[i]) || !isfinite(diag.dffits[i])) allFinite = false;
        printf("numerical check: Cook's D and DFFITS finite for all rows: %s\n", allFinite ? "true" : "false");

        double leverageThresh = 2.0 * ols.p / n;
        double cookThresh = 4.0 / n;
        printf("leverage threshold 2p/n = %.4f, Cook's D threshold 4/n = %.4f\n", leverageThresh, cookThresh);
    }

    printf("\n=== Scenario B: predictive metrics, inference, and a model-level diagnostic on a 150/60 train/test split ===\n");
    OLSRegressor baseline;
    Matrix Xtr, Xte;
    vector<double> ytr, yte;
    int Ntrain = 150, Ntest = 60;
    {
        mt19937 rng(7);
        normal_distribution<double> noise(0.0, 1.0);
        uniform_real_distribution<double> ux(-3.0, 3.0);
        for (int i = 0; i < Ntrain; ++i) {
            double x1 = ux(rng), x2 = ux(rng);
            Xtr.push_back({x1, x2});
            ytr.push_back(5.0 + 1.5 * x1 - 2.0 * x2 + noise(rng));
        }
        for (int i = 0; i < Ntest; ++i) {
            double x1 = ux(rng), x2 = ux(rng);
            Xte.push_back({x1, x2});
            yte.push_back(5.0 + 1.5 * x1 - 2.0 * x2 + noise(rng));
        }

        baseline.fit(Xtr, ytr);
        auto predTr = baseline.predict(Xtr);
        auto predTe = baseline.predict(Xte);
        printf("baseline fit: intercept=%.4f coef=[%.4f, %.4f]\n", baseline.intercept, baseline.coef[0], baseline.coef[1]);
        printf("train: MSE=%.4f RMSE=%.4f MAE=%.4f R2=%.4f\n", mse(ytr, predTr), rmseOf(ytr, predTr), maeOf(ytr, predTr), r2Score(ytr, predTr));
        printf("test:  MSE=%.4f RMSE=%.4f MAE=%.4f R2=%.4f\n", mse(yte, predTe), rmseOf(yte, predTe), maeOf(yte, predTe), r2Score(yte, predTe));
        printf("adj.R2(train)=%.4f  RSE=%.4f (df=%d)\n",
               adjustedR2(r2Score(ytr, predTr), Ntrain, baseline.d), residualStandardError(baseline.sse, Ntrain, baseline.p), Ntrain - baseline.p);

        printf("\n-- OLS inference on the baseline fit: SE, t, normal-approx 95%% CI --\n");
        double s2 = baseline.s * baseline.s;
        vector<double> beta = {baseline.intercept, baseline.coef[0], baseline.coef[1]};
        const char* names[3] = {"intercept", "x1", "x2"};
        const double z975 = 1.959964;
        for (int j = 0; j < baseline.p; ++j) {
            double se = sqrt(s2 * baseline.XtXinv[j][j]);
            double t = beta[j] / se;
            double lo = beta[j] - z975 * se, hi = beta[j] + z975 * se;
            printf("%-10s beta=%.4f SE=%.4f t=%.4f 95%%CI=[%.4f, %.4f]\n", names[j], beta[j], se, t, lo, hi);
        }

        printf("\n-- Breusch-Pagan heteroskedasticity check --\n");
        vector<double> residTr(Ntrain);
        for (int i = 0; i < Ntrain; ++i) residTr[i] = ytr[i] - predTr[i];
        double lm = breuschPaganLM(Xtr, residTr);
        double bpThresh = 2.0 * baseline.d;
        printf("LM = %.4f, heuristic threshold 2d = %.4f -> %s\n", lm, bpThresh, (lm > bpThresh ? "flagged" : "not flagged"));

        printf("\n-- PRESS / LOOCV RMSE on the baseline fit --\n");
        Diagnostics diagB = diagnose(Xtr, ytr, baseline);
        double pressSum = 0.0;
        for (double pr : diagB.press) pressSum += pr * pr;
        printf("LOOCV RMSE (via PRESS identity) = %.4f (vs in-sample RMSE %.4f)\n", sqrt(pressSum / Ntrain), rmseOf(ytr, predTr));
        vector<int> checkIdx = {0, 10, 50, 100, 149};
        double maxDiff = 0.0;
        for (int idx : checkIdx) {
            Matrix Xloo;
            vector<double> yloo;
            for (int i = 0; i < Ntrain; ++i) {
                if (i == idx) continue;
                Xloo.push_back(Xtr[i]);
                yloo.push_back(ytr[i]);
            }
            OLSRegressor loo;
            loo.fit(Xloo, yloo);
            double explicitResid = ytr[idx] - loo.predict(Xtr[idx]);
            double diff = fabs(explicitResid - diagB.press[idx]);
            maxDiff = max(maxDiff, diff);
            printf("  row %3d: PRESS=%.6f explicit-LOO-residual=%.6f |diff|=%.2e\n", idx, diagB.press[idx], explicitResid, diff);
        }
        printf("max |PRESS - explicit LOO residual| over checked rows = %.2e\n", maxDiff);
    }

    printf("\n=== Scenario C: injecting one high-leverage, high-residual row and refitting without it ===\n");
    {
        Matrix XtrExt = Xtr;
        vector<double> ytrExt = ytr;
        XtrExt.push_back({7.0, 7.0});
        ytrExt.push_back(5.0 + 1.5 * 7.0 - 2.0 * 7.0 + 14.0);
        OLSRegressor withInj;
        withInj.fit(XtrExt, ytrExt);
        printf("with-injection fit: intercept=%.4f coef=[%.4f, %.4f]\n", withInj.intercept, withInj.coef[0], withInj.coef[1]);

        auto predTrBase = baseline.predict(Xtr);
        auto predTrInj = withInj.predict(Xtr);
        auto predTeBase = baseline.predict(Xte);
        auto predTeInj = withInj.predict(Xte);
        printf("clean 150 train rows, baseline fit:      MSE=%.4f RMSE=%.4f MAE=%.4f R2=%.4f\n",
               mse(ytr, predTrBase), rmseOf(ytr, predTrBase), maeOf(ytr, predTrBase), r2Score(ytr, predTrBase));
        printf("clean 150 train rows, with-injection fit: MSE=%.4f RMSE=%.4f MAE=%.4f R2=%.4f\n",
               mse(ytr, predTrInj), rmseOf(ytr, predTrInj), maeOf(ytr, predTrInj), r2Score(ytr, predTrInj));
        printf("60 test rows, baseline fit:      MSE=%.4f RMSE=%.4f MAE=%.4f R2=%.4f\n",
               mse(yte, predTeBase), rmseOf(yte, predTeBase), maeOf(yte, predTeBase), r2Score(yte, predTeBase));
        printf("60 test rows, with-injection fit: MSE=%.4f RMSE=%.4f MAE=%.4f R2=%.4f\n",
               mse(yte, predTeInj), rmseOf(yte, predTeInj), maeOf(yte, predTeInj), r2Score(yte, predTeInj));

        int Next = (int)XtrExt.size();
        Diagnostics diagExt = diagnose(XtrExt, ytrExt, withInj);
        double sumH = 0.0;
        for (double v : diagExt.h) sumH += v;
        printf("numerical check: sum(h_ii) = %.6f, expected p = %d\n", sumH, withInj.p);
        int last = Next - 1;
        printf("injected row (index %d): leverage=%.4f resid=%.4f std.res=%.4f Cook's D=%.4f DFFITS=%.4f PRESS=%.4f\n",
               last, diagExt.h[last], diagExt.resid[last], diagExt.stdResid[last], diagExt.cooksD[last],
               diagExt.dffits[last], diagExt.press[last]);
        double maxCookOther = 0.0;
        for (int i = 0; i < Next; ++i) { if (i != last) maxCookOther = max(maxCookOther, diagExt.cooksD[i]); }
        printf("largest Cook's D among the other %d rows = %.4f (vs %.4f for the injected row)\n", Next - 1, maxCookOther, diagExt.cooksD[last]);
        double leverageThresh = 2.0 * withInj.p / Next, cookThresh = 4.0 / Next, dffitsThresh = 2.0 * sqrt((double)withInj.p / Next);
        printf("thresholds: leverage 2p/n=%.4f, Cook's D 4/n=%.4f, |DFFITS| 2*sqrt(p/n)=%.4f\n", leverageThresh, cookThresh, dffitsThresh);

        OLSRegressor removed;
        removed.fit(Xtr, ytr);
        printf("refit with the injected row removed: intercept=%.4f coef=[%.4f, %.4f] (recovers the baseline fit exactly)\n",
               removed.intercept, removed.coef[0], removed.coef[1]);
    }

    printf("\n=== Scenario D: a small, high-variance fit with a negative test R^2 ===\n");
    {
        mt19937 rng(3);
        normal_distribution<double> bigNoise(0.0, 6.0);
        uniform_real_distribution<double> ux(-3.0, 3.0);
        int nTrainSmall = 5, nTestSmall = 30;
        Matrix XtrSmall, XteSmall;
        vector<double> ytrSmall, yteSmall;
        for (int i = 0; i < nTrainSmall; ++i) {
            double x1 = ux(rng), x2 = ux(rng);
            XtrSmall.push_back({x1, x2});
            ytrSmall.push_back(5.0 + 1.5 * x1 - 2.0 * x2 + bigNoise(rng));
        }
        for (int i = 0; i < nTestSmall; ++i) {
            double x1 = ux(rng), x2 = ux(rng);
            XteSmall.push_back({x1, x2});
            yteSmall.push_back(5.0 + 1.5 * x1 - 2.0 * x2 + bigNoise(rng));
        }
        OLSRegressor tiny;
        tiny.fit(XtrSmall, ytrSmall);
        auto predTrSmall = tiny.predict(XtrSmall);
        auto predTeSmall = tiny.predict(XteSmall);
        printf("fit on %d noisy rows (p=%d, df=%d): intercept=%.4f coef=[%.4f, %.4f]\n",
               nTrainSmall, tiny.p, nTrainSmall - tiny.p, tiny.intercept, tiny.coef[0], tiny.coef[1]);
        printf("train: MSE=%.4f RMSE=%.4f MAE=%.4f R2=%.4f\n",
               mse(ytrSmall, predTrSmall), rmseOf(ytrSmall, predTrSmall), maeOf(ytrSmall, predTrSmall), r2Score(ytrSmall, predTrSmall));
        printf("test:  MSE=%.4f RMSE=%.4f MAE=%.4f R2=%.4f\n",
               mse(yteSmall, predTeSmall), rmseOf(yteSmall, predTeSmall), maeOf(yteSmall, predTeSmall), r2Score(yteSmall, predTeSmall));
    }

    printf("\n=== Scenario E: explicit failure modes ===\n");
    {
        try {
            Matrix Xbad = {{1.0}, {2.0}, {3.0}};
            vector<double> ybad = {1.0, 2.0};
            OLSRegressor bad;
            bad.fit(Xbad, ybad);
        } catch (const invalid_argument& ex) { printf("caught: %s\n", ex.what()); }
        try {
            Matrix Xbad = {{1.0, 2.0}, {2.0, 4.0}};
            vector<double> ybad = {1.0, 2.0};
            OLSRegressor bad;
            bad.fit(Xbad, ybad);
        } catch (const invalid_argument& ex) { printf("caught: %s\n", ex.what()); }
        try {
            Matrix Xcol = {{1.0, 2.0}, {2.0, 4.0}, {3.0, 6.0}, {4.0, 8.0}};
            vector<double> ycol = {1.0, 2.0, 3.5, 3.9};
            OLSRegressor bad;
            bad.fit(Xcol, ycol);
        } catch (const runtime_error& ex) { printf("caught: %s\n", ex.what()); }

        printf("\n-- non-finite input and dimension-mismatch guards --\n");
        try {
            Matrix Xnan = {{1.0}, {numeric_limits<double>::quiet_NaN()}, {3.0}, {4.0}};
            vector<double> yok = {1.0, 2.0, 3.0, 4.0};
            OLSRegressor bad;
            bad.fit(Xnan, yok);
        } catch (const invalid_argument& ex) { printf("caught: %s\n", ex.what()); }
        try {
            Matrix Xok = {{1.0}, {2.0}, {3.0}, {4.0}};
            vector<double> yinf = {1.0, 2.0, numeric_limits<double>::infinity(), 4.0};
            OLSRegressor bad;
            bad.fit(Xok, yinf);
        } catch (const invalid_argument& ex) { printf("caught: %s\n", ex.what()); }

        Matrix Xgood = {{1.0, 2.0}, {2.0, 3.0}, {3.0, 5.0}, {4.0, 6.0}};
        vector<double> ygood = {1.0, 2.0, 3.0, 4.0};
        OLSRegressor good;
        good.fit(Xgood, ygood);
        try {
            vector<double> wrongX = {1.0};
            good.predict(wrongX);
        } catch (const invalid_argument& ex) { printf("caught: %s\n", ex.what()); }
        try {
            Matrix Xwrong = {{1.0, 2.0}, {2.0, 3.0}};
            vector<double> ywrong = {1.0, 2.0};
            diagnose(Xwrong, ywrong, good);
        } catch (const invalid_argument& ex) { printf("caught: %s\n", ex.what()); }

        printf("\n-- predictive-metric and degrees-of-freedom guards --\n");
        try {
            vector<double> yt = {1.0, 2.0, 3.0}, yp = {1.0, 2.0};
            r2Score(yt, yp);
        } catch (const invalid_argument& ex) { printf("caught: %s\n", ex.what()); }
        try {
            vector<double> yt, yp;
            mse(yt, yp);
        } catch (const invalid_argument& ex) { printf("caught: %s\n", ex.what()); }
        try {
            vector<double> yt = {4.0, 4.0, 4.0}, yp = {4.0, 4.1, 3.9};
            r2Score(yt, yp);
        } catch (const invalid_argument& ex) { printf("caught: %s\n", ex.what()); }
        try {
            adjustedR2(0.9, 5, 4);
        } catch (const invalid_argument& ex) { printf("caught: %s\n", ex.what()); }
        try {
            residualStandardError(12.0, 5, 5);
        } catch (const invalid_argument& ex) { printf("caught: %s\n", ex.what()); }
    }
    return 0;
}

Compiled with g++ -O2 -std=c++17 (also verified clean under -Wall -Wextra -Wpedantic -Wshadow), this prints, byte for byte:

A Practical Diagnostic Workflow

• Predictive metrics answer whether the model works; diagnostics answer why. Scenario C's own table makes the distinction unavoidable: test $R^2$ moves by about one percentage point when a single influential row is added, an amount easy to write off as noise, while that same row's Cook's distance sits nearly $200\times$ above every other row's. A workflow that only tracks RMSE or $R^2$ over time can miss a structural problem for a long stretch of retraining runs before it becomes large enough to show up in the aggregate number.
• Leverage is available before a single target value is collected, so an early check on $X$ alone is inexpensive. Because $h_{ii}$ depends only on the design matrix, auditing it the moment features are assembled catches rows sitting unusually far out in feature space independent of whether their eventual target values turn out to be well-behaved.
• Cook's distance and DFFITS, not leverage or a residual in isolation, settle which rows actually move the fit. A high-leverage row with a residual consistent with the trend stays low on both; only the combination the multiplicative structure of $D_i$ and $\text{DFFITS}_i$ captures is worth escalating.
• $\text{RMSE}_{\text{LOOCV}}$ is a more honest error estimate than in-sample RMSE whenever leverage is uneven. Scenario B's LOOCV RMSE, $1.0008$, sits modestly above its in-sample RMSE, $0.9823$, because a few rows carry above-average leverage; the PRESS identity makes that comparison available at the cost of one pass over already-cached quantities, not $n$ extra model fits.
• The Breusch-Pagan check is worth running once per fit, since every standard error above assumes its null hypothesis. A homoskedastic fit, like Scenario B's synthetic data, produces a small LM comfortably under the heuristic threshold; a large LM is a direct warning that the confidence intervals reported alongside it understate the true uncertainty in some region of $X$.

Production Guidance

• Leverage and Cook's distance belong inside the training pipeline itself, not in an offline notebook step. Both are $O(np^2)$ given the already-cached $(X_b^\top X_b)^{-1}$, cheap relative to the fit itself, and logging the top few rows by Cook's distance on every retraining run turns a one-off audit into a continuously monitored signal.
• LOOCV RMSE alongside in-sample RMSE is worth reporting whenever the training set is small enough that a handful of rows carry meaningfully above-average leverage. The gap between the two numbers is itself informative: a large gap is a warning that in-sample error is optimistic in a way plain train/test splitting on a small dataset may not fully reveal.
• A flagged row is a prompt for investigation, not an automatic deletion candidate. Scenario C's injected row is a genuine data problem by construction, but a real flagged row might instead be a rare, correctly recorded, and entirely legitimate observation; removing it without checking would quietly bias the model against exactly the cases it was built to handle.
• The Breusch-Pagan check is worth re-running after any feature-engineering change, not only once at initial model build. Adding, removing, or transforming a feature changes the auxiliary regression entirely, and a fit that was homoskedastic under one feature set is not guaranteed to remain so under another.
• The normal-approximation confidence intervals stay honest only when their own assumptions travel with them into any report that surfaces them. They assume $n-p$ is large enough for the normal and $t$ distributions to coincide and that the Breusch-Pagan check passed; presenting them without either caveat overstates the certainty of a coefficient estimate to whoever reads the report next.

Complexity and Memory

The dominant one-time cost is the same $O(np^2+p^3)$ every earlier post's normal-equation solve already pays; leverage, the studentized residual, Cook's distance, DFFITS, and PRESS all reuse that single cached inverse and add only linear-in-$n$ work per diagnostic. The explicit leave-one-out refit is the one operation in this post whose cost scales the way a naive (non-PRESS) implementation would: it exists only to verify the PRESS identity against ground truth, not as part of the production diagnostic path, and the complete program above deliberately checks it on a handful of representative rows in the $150$-row scenario rather than on all $150$, since the identity itself already provides the exact answer at $O(1)$ per row.

Assumptions and Limitations

• Every formula in this post assumes the unpenalized OLS estimator under approximately Gaussian, homoskedastic, independent errors. The scope section above is not optional context: applying the coefficient covariance formula to a ridge or lasso fit, or the leverage/hat-matrix machinery to a Huber or quantile fit, produces numbers that look plausible and are not valid for those estimators.
• Rule-of-thumb thresholds are triage, not hypothesis tests with a controlled error rate. The leverage threshold $2p/n$, the Cook's-distance threshold $4/n$, the DFFITS threshold $2\sqrt{p/n}$, and the Breusch-Pagan heuristic threshold $2d$ are all useful starting points for investigation, not proof that a flagged row is erroneous or that an unflagged row is safe.
• The normal-approximation confidence interval is an approximation that degrades as $n-p$ shrinks. It is a reasonable stand-in for a $t$-distribution interval once residual degrees of freedom are comfortably large, which is why this post reports it only for the $150$-row scenario, not the eight-point one, where $n-p=6$ is too small for the normal and $t$ distributions to be close.
• The Breusch-Pagan heuristic threshold stands in for a proper $\chi^2_d$ comparison. A from-scratch, standard-library-only implementation has no chi-squared quantile function available; $2d$ is a rough guide to the same idea, not a $p$-value, and a borderline LM value deserves a proper statistical test before being treated as conclusive.
• PRESS and $\text{RMSE}_{\text{LOOCV}}$ are exact only for this unpenalized OLS estimator. The identity $\text{PRESS}_i=e_i/(1-h_{ii})$ is derived from how the hat matrix updates under single-row removal; a regularized or reweighted fit does not update the same way under row removal, so the identity does not carry over without rederiving it for that estimator's own influence function.
• DFFITS's deleted variance $s_{(i)}^2$ needs $n-p-1>0$ residual degrees of freedom and is not guaranteed positive for every row. Both are checked explicitly rather than clipped or silently defaulted: a fit with $n-p-1\le 0$ cannot support DFFITS at all, and a row extreme enough that $e_i^2/(1-h_{ii})$ meets or exceeds the whole fit's $\text{SSE}$ leaves $s_{(i)}^2$ non-positive, at which point DFFITS is genuinely undefined for that row rather than a number a fallback could paper over.
• None of this replaces looking at the data directly. Every diagnostic here is a summary statistic; a residual-versus-fitted plot, the leverage-versus-Cook's-distance view animated above, or simply sorting rows by Cook's distance and reading the raw values remains the fastest way to catch a genuine data problem a summary number can understate.

Next Step

This closes the linear-regression sub-series: closed-form, numerically robust, iterative, penalized, loss-robust, and now measured and diagnosed on the same footing. The next model family in the from-scratch series keeps the same $\hat y=\beta_0+\sum_j\beta_j x_j$ linear form but changes what it predicts: logistic regression replaces a continuous target with a binary one, trading the normal-equations closed form for a sigmoid link and log-loss, and reuses several ideas from this series directly, IRLS from the robust regression post reappears nearly unchanged as the standard logistic-regression solver, and this post's leverage and influence machinery has a direct classification-side analogue once the loss changes from squared error to log-loss.