Skip to content

Cross-Validation & Diagnostics

This tutorial demonstrates K-fold cross-validated model selection and post-fit diagnostic inspection.

Generate data

import numpy as np

rng = np.random.default_rng(42)

n_features, n_samples, true_rank = 40, 150, 4
noise_std = 0.5

W = rng.standard_normal((n_features, true_rank))
S = rng.standard_normal((true_rank, n_samples))
X = W @ S + noise_std * rng.standard_normal((n_features, n_samples))

# 15% missing
mask = rng.random(X.shape) > 0.15
X_obs = np.where(mask, X, np.nan)

Cross-validate component count

cross_validate_components partitions the observed entries (not full rows) into folds. Each fold holds out a subset of entries, fits on the rest, and evaluates reconstruction on the held-out entries. The current implementation accepts dense inputs and guarantees that each training fold retains at least one observation in every non-empty row and column. CVConfig.seed controls both fold assignment and model initialization.

from vbpca_py import cross_validate_components, CVConfig

cfg = CVConfig(
    n_splits=5,
    metric="prms",       # probe RMS on held-out entries
    one_se_rule=True,    # prefer simpler model within 1 SE of best
    seed=42,
)

best_k, results = cross_validate_components(
    X_obs,
    mask=mask,
    components=range(0, 10),  # include the mean-only model
    config=cfg,
    maxiters=200,
)

print(f"Selected k = {best_k}  (true rank = {true_rank})")

Inspect CV results

for entry in results:
    k = entry["k"]
    mean_metric = entry["mean_prms"]
    se = entry["se_prms"]
    marker = " <-- selected" if k == best_k else ""
    print(f"  k={k:2d}  prms={mean_metric:.4f} ± {se:.4f}{marker}")

1-SE rule

When one_se_rule=True, the selected \(k\) is the smallest value whose mean metric is within one standard error of the overall best. This guards against overfitting by favouring parsimony. Rank zero denotes a mean-only model and is a valid result when the held-out data do not support a latent component. The one-standard-error rule is a model-selection heuristic rather than a hypothesis test.

Fit the final model

The simulated example has positive latent rank. If best_k == 0 on another dataset, report the mean-only result; there is no positive-rank VBPCA object to fit. Otherwise, fit the selected positive-rank model as below.

from vbpca_py import VBPCA, make_xprobe_mask

# Create a probe set for monitoring
X_train, X_probe = make_xprobe_mask(X_obs, fraction=0.10, rng=rng)

model = VBPCA(n_components=best_k, maxiters=300, verbose=1)
model.fit(X_train, mask=~np.isnan(X_train), xprobe=X_probe)

Convergence diagnostics

Learning curves

print(f"Final RMS:       {model.rms_:.6f}")
print(f"Final probe RMS: {model.prms_:.6f}")
print(f"Final cost:      {model.cost_:.2f}")

Explained variance

evr = model.explained_variance_ratio_
if evr is not None:
    total = 0.0
    for i, v in enumerate(evr):
        total += v
        print(f"  PC{i+1}: {v:.3f}  (cumulative: {total:.3f})")

Plotting

Install the plotting extra if you haven't already:

pip install vbpca_py[plot]

from vbpca_py.plotting import scree_plot, loadings_barplot, variance_explained_plot

# Scree plot — explained variance ratio per component
fig = scree_plot(model, cumulative=True)

# Loadings bar plot — feature importance for the first component
fig = loadings_barplot(model, component=0, top_n=15)

# Absolute variance per component
fig = variance_explained_plot(model)

Scree plot

The scree plot shows the explained variance ratio for each component. A sharp "elbow" suggests the number of meaningful components. The cumulative line shows the total variance explained.

Loadings bar plot

The loadings bar plot for a given component highlights which features contribute most. Use top_n to show only the highest-magnitude features.

Next steps