pod.compute_kriging_loo_residuals
compute_kriging_loo_residuals(gpr, X_2d, y)Computes exact LOO predictions, variances, standardized residuals e_i, and the outlier factor gamma, following Malkiel et al. (2026).
With B the top-left block of S^-1, S = [[K + alpha*I, F], [F^T, 0]]:
- LOO mean ([M26] Eq. 10): mu_{-i} = -sum_{j != i} (B_ij / B_ii) y_j
- LOO variance ([M26] Eq. 11): sigma_{-i}^2 = 1 / B_ii
- Standardised residual ([M26] Eq. 19): e_i = (y_i - mu_{-i}) / sigma_{-i}
- Outlier factor ([M26] Sec. 2.2.5.2): gamma = max(1, max|e_i| / 3)
As in [M26], the hyperparameters are those fitted to the full data set; they are not re-optimised with point i removed.
K is the fitted kernel evaluated at the training inputs, so for the models built by fit_all_robust_mean_models it already includes the learned noise (WhiteKernel) variance on its diagonal. The LOO variance is therefore the predictive variance of a new noisy observation, which is the correct scale for standardizing observed residuals.
In [M26], gamma scales the Kriging interpolation standard deviation used to draw the GP sample paths for the PoD uncertainty bound (Fig. 10b shows the residuals divided by gamma, so the worst one lands on 3). digiqual computes its PoD bounds by bootstrap ([M25]) instead, so gamma is reported as a diagnostic and is not applied to the PoD [digiqual].
If the model was fitted with normalize_y=True, the kernel hyperparameters are in normalised units. The calculation is then done on the normalised response and the LOO means and standard deviations are converted back to the original units. Inputs are standardised first if the model does so (ScaledGaussianProcessRegressor).
Parameters
| Name | Type | Description | Default |
|---|---|---|---|
| gpr | GaussianProcessRegressor | A fitted scikit-learn Gaussian Process model. | required |
| X_2d | np.ndarray | 2D matrix of input training coordinates (N_samples, N_features). | required |
| y | np.ndarray | 1D array of observed responses (N_samples,). | required |
Returns
| Name | Type | Description |
|---|---|---|
| tuple[np.ndarray, np.ndarray, np.ndarray, float] | tuple[np.ndarray, np.ndarray, np.ndarray, float]: - loo_means: Array of Leave-One-Out predicted mean responses. - loo_stds: Array of Leave-One-Out predicted standard deviations. - std_residuals: Array of standardized LOO residuals e_i = (y_i - mu_{-i}) / sigma_{-i}. - gamma: Outlier factor gamma = max(1.0, max|e_i| / 3.0). |
Examples
import numpy as np
from sklearn.gaussian_process import GaussianProcessRegressor
from digiqual.pod import compute_kriging_loo_residuals
X = np.linspace(0, 5, 20).reshape(-1, 1)
y = 2.0 * X.flatten() + np.random.normal(0, 0.1, 20)
gpr = GaussianProcessRegressor()
gpr.fit(X, y)
loo_means, loo_stds, std_res, gamma = compute_kriging_loo_residuals(gpr, X, y)
print(f"Outlier factor gamma: {gamma:.3f}")