Approach: use nonparametric case bootstrap (resample rows (X,y) with replacement), refit or (if model expensive) use residual bootstrap when assumptions hold. For each bootstrap sample, refit the provided scikit-learn regressor and predict responses at the specified feature points X0. Collect predicted means across B resamples and compute the 2.5 and 97.5 percentiles → 95% CI.Assumptions:- Data rows are i.i.d.- Model class is appropriate for refitting on resamples.- If using residual bootstrap: model is approximately correctly specified and residuals are homoskedastic.Recommended resamples: B = 1000–2000 for stable 95% CIs; 500 may be acceptable for exploratory work.Computational trade-offs:- Refitting B times is expensive (O(B * fit_time)); consider parallelism (joblib), subsampling, or residual bootstrap if model is linear-ish.- Store only predictions to save memory.- Use a fixed random_state for reproducibility.Code (Python pseudocode / runnable pattern):python
import numpy as np
from sklearn.base import clone
from joblib import Parallel, delayed
import pandas as pd
def bootstrap_predict_ci(model, X, y, X0, B=1000, alpha=0.05, n_jobs=1, random_state=None):
"""
model: fitted scikit-learn regressor
X, y: training data (numpy or pandas)
X0: feature points to predict (array shape (m, p))
returns: DataFrame with mean_pred, ci_lower, ci_upper for each row in X0
"""
rng = np.random.RandomState(random_state)
n = len(X)
preds = np.zeros((B, len(X0)))
def one_boot(b):
idx = rng.randint(0, n, size=n) # resample indices with replacement
Xb = X[idx]
yb = y[idx]
m = clone(model)
m.fit(Xb, yb)
return m.predict(X0)
results = Parallel(n_jobs=n_jobs)(delayed(one_boot)(b) for b in range(B))
preds = np.vstack(results) # shape (B, m)
mean_pred = preds.mean(axis=0)
lower = np.percentile(preds, 100 * (alpha / 2), axis=0)
upper = np.percentile(preds, 100 * (1 - alpha / 2), axis=0)
return pd.DataFrame({
'mean_pred': mean_pred,
'ci_lower': lower,
'ci_upper': upper
}, index=range(len(X0)))
