Skip to content

[ENH] Bootstrap for ALE confidence intervals #703

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
abaziree wants to merge 2 commits into
base: 297-ale
Choose a base branch
from
Open

[ENH] Bootstrap for ALE confidence intervals #703

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
c3bfa5e
add bootstrap
abaziree Jun 29, 2026
485142b
update from base ALE
abaziree Jul 2, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
193 changes: 115 additions & 78 deletions src/hidimstat/visualization/accumulated_local_effects.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@
import numpy as np
import pandas as pd
import seaborn as sns
from joblib import Parallel, delayed
from scipy import stats
from sklearn.utils.validation import check_is_fitted

Expand Down Expand Up @@ -163,9 +164,11 @@ def compute_ale_1d(
feature_idx,
feature_type,
grid_resolution="auto",
percentiles=(5, 95),
confidence_interval=True,
confidence_level=0.95,
percentiles=(5, 95),
n_bootstraps=20,
n_jobs=-1,
):
"""Compute the 1D Accumulated Local Effect for a single numerical feature.

Expand Down Expand Up @@ -197,13 +200,17 @@ def compute_ale_1d(
if the data contains many duplicate values or fewer unique points
than requested.

percentiles : tuple of float, default=(5, 95)
The lower and upper percentile used to create the extreme values for the grid.
Must be in [0, 100].
confidence_interval : bool, default=True
Whether to compute the confidence intervals of the ALE curve.
confidence_level : float, default=0.95
The confidence level used to compute the confidence intervals (e.g., 0.95 for 95%).
percentiles : tuple of float, default=(5, 95)
The lower and upper percentile used to create the extreme values for the grid.
Must be in [0, 100].
n_bootstraps : int, default=20
Number of bootstrap samples to generate if `confidence_interval` is True.
n_jobs : int, default=-1
Number of jobs to run in parallel during bootstrapping. `-1` means using all processors.

Returns
-------
Expand Down Expand Up @@ -231,23 +238,10 @@ def compute_ale_1d(
f"Feature {feature_idx} has fewer than 2 unique quantile edges. Increase grid_resolution or check your data."
)

bin_idx = _bin_indices(x, grid_values)

# For each sample, evaluate the model at the lower and upper edge of its bin
X_low = X.copy()
X_high = X.copy()
X_low[:, feature_idx] = grid_values[bin_idx]
X_high[:, feature_idx] = grid_values[bin_idx + 1]

local_effects = _predict_fn(estimator, X_high) - _predict_fn(
estimator, X_low
) # shape (n_samples,)

combined_idx = bin_idx
combined_effects = local_effects
min_weight_bin = 0
X_bootstrap = X

elif feature_type == "categorical":
# Grid defined by unique values
if (
not isinstance(percentiles, tuple)
or len(percentiles) != 2
Expand All @@ -258,7 +252,6 @@ def compute_ale_1d(
"in [0, 1] in increasing order"
)

# Grid defined by unique values
low_bnd = np.percentile(x, percentiles[0])
high_bnd = np.percentile(x, percentiles[1])

Expand All @@ -274,76 +267,111 @@ def compute_ale_1d(
f"Feature {feature_idx} has fewer than 2 unique values. Check your data."
)

value_idx = np.digitize(x_filtered, grid_values) - 1
X_bootstrap = X_filtered

else:
raise ValueError(
"'feature_type' must be a string among 'continuous' and 'categorical'."
)

def _compute_ale_curve(X_sample):
"""Internal closure to compute the ALE curve on a fixed grid."""
if feature_type == "continuous":
x_sample = X_sample[:, feature_idx]
bin_idx = _bin_indices(x_sample, grid_values)

# For each sample, evaluate the model at the lower and upper edge of its bin
X_low = X_sample.copy()
X_high = X_sample.copy()
X_low[:, feature_idx] = grid_values[bin_idx]
X_high[:, feature_idx] = grid_values[bin_idx + 1]

# For each sample, evaluate the model at the lower and upper edge of its bin
X_low = X_filtered.copy()
X_high = X_filtered.copy()
mask_low = x_filtered != grid_values[0]
mask_high = x_filtered != grid_values[-1]
X_low[mask_low, feature_idx] = grid_values[value_idx[mask_low] - 1]
X_high[mask_high, feature_idx] = grid_values[value_idx[mask_high] + 1]
local_effects = _predict_fn(estimator, X_high) - _predict_fn(
estimator, X_low
) # shape (n_samples,)

local_effects_low = _predict_fn(
estimator, X_filtered[mask_low]
) - _predict_fn(estimator, X_low[mask_low])
local_effects_high = _predict_fn(
estimator, X_high[mask_high]
) - _predict_fn(estimator, X_filtered[mask_high])
combined_idx = bin_idx
combined_effects = local_effects
min_weight_bin = 0

else:
x_sample = X_sample[:, feature_idx]
value_idx = np.digitize(x_sample, grid_values) - 1

# For each sample, evaluate the model at the lower and upper edge of its bin
X_low = X_sample.copy()
X_high = X_sample.copy()
mask_low = x_sample != grid_values[0]
mask_high = x_sample != grid_values[-1]
X_low[mask_low, feature_idx] = grid_values[value_idx[mask_low] - 1]
X_high[mask_high, feature_idx] = grid_values[
value_idx[mask_high] + 1
]

local_effects_low = _predict_fn(
estimator, X_sample[mask_low]
) - _predict_fn(estimator, X_low[mask_low])
local_effects_high = _predict_fn(
estimator, X_high[mask_high]
) - _predict_fn(estimator, X_sample[mask_high])

# Average local effects within each bin
combined_idx = np.concatenate(
[value_idx[mask_low], value_idx[mask_high] + 1]
)
combined_effects = np.concatenate(
[local_effects_low, local_effects_high]
)
min_weight_bin = 1

# Average local effects within each bin
combined_idx = np.concatenate(
[value_idx[mask_low], value_idx[mask_high] + 1]
bin_counts = np.bincount(combined_idx, minlength=n_bins).astype(float)
bin_sums = np.bincount(
combined_idx, weights=combined_effects, minlength=n_bins
)
combined_effects = np.concatenate(
[local_effects_low, local_effects_high]
)
min_weight_bin = 1

else:
raise ValueError(
"'feature_type' must be a string among 'continuous' and 'categorical'."
)
mean_effects = np.zeros(n_bins, dtype=float)
non_zero = bin_counts > 0
mean_effects[non_zero] = bin_sums[non_zero] / bin_counts[non_zero]

# Average local effects within each bin
bin_counts = np.bincount(combined_idx, minlength=n_bins).astype(float)
bin_sums = np.bincount(
combined_idx, weights=combined_effects, minlength=n_bins
)
# Cumulative sum: uncentered ALE evaluated for each bin edge
if feature_type == "continuous":
ale_curve = np.array([0, *np.cumsum(mean_effects)])
else:
ale_curve = np.cumsum(mean_effects)

mean_effects = np.zeros(n_bins, dtype=float)
non_zero = bin_counts > 0
mean_effects[non_zero] = bin_sums[non_zero] / bin_counts[non_zero]
# Center: subtract the sample-weighted mean
ale_centers = (ale_curve[1:] + ale_curve[:-1]) / 2
ale_curve -= (
np.sum(ale_centers * bin_counts[min_weight_bin:])
/ bin_counts.sum()
)

# Cumulative sum: uncentered ALE evaluated for each bin edge
if feature_type == "continuous":
ale = np.array([0, *np.cumsum(mean_effects)])
else:
ale = np.cumsum(mean_effects)
return ale_curve

# Center: subtract the sample-weighted mean
ale_centers = (ale[1:] + ale[:-1]) / 2
ale -= np.sum(ale_centers * bin_counts[min_weight_bin:]) / bin_counts.sum()
# Base calculation
ale = _compute_ale_curve(X_bootstrap)
ale_err = None

# Confidence interval
ale_err = None
if confidence_interval:
sample_means = mean_effects[combined_idx]
squared_deviations = (combined_effects - sample_means) ** 2
sum_squared_deviations = np.bincount(
combined_idx, weights=squared_deviations, minlength=n_bins
)
if n_bootstraps < 1:
raise ValueError("'n_bootstrap' must be strictly greater than 0.")

var_of_mean = np.zeros(n_bins, dtype=float)
valid_bins = bin_counts > 1
var_of_mean[valid_bins] = sum_squared_deviations[valid_bins] / (
bin_counts[valid_bins] * (bin_counts[valid_bins] - 1)
def _bootstrap():
bootstrap_indices = np.random.default_rng().choice(
len(X_bootstrap), size=len(X_bootstrap), replace=True
)
return _compute_ale_curve(X_bootstrap[bootstrap_indices])

bootstrap_curves = Parallel(n_jobs=n_jobs)(
delayed(_bootstrap)() for _ in range(n_bootstraps)
)
bootstrap_curves = np.array(bootstrap_curves)

if feature_type == "continuous":
var_of_mean = np.array([0, *var_of_mean])
ale = np.mean(bootstrap_curves, axis=0)
z_score = stats.norm.ppf(1 - (1 - confidence_level) / 2)
ale_err = z_score * np.sqrt(var_of_mean)
ale_err = z_score * np.std(bootstrap_curves, axis=0)

return ale, grid_values, ale_err

Expand Down Expand Up @@ -654,9 +682,11 @@ def plot(
features,
feature_type="auto",
grid_resolution="auto",
percentiles=(5, 95),
confidence_interval=True,
confidence_level=0.95,
percentiles=(5, 95),
n_bootstraps=20,
n_jobs=-1,
cmap="viridis",
**kwargs,
):
Expand Down Expand Up @@ -686,13 +716,18 @@ def plot(
strictly less than `grid_resolution` (or the auto-calculated value)
if the data contains many duplicate values or fewer unique points
than requested.

percentiles : tuple of float, default=(5, 95)
The lower and upper percentile used to create the extreme values for the grid.
Must be in [0, 100].
confidence_interval : bool, default=True
Whether to compute and display confidence intervals around the 1D ALE curve.
confidence_level : float, default=0.95
The confidence level used to compute the confidence intervals (e.g., 0.95 for 95%).
percentiles : tuple of float, default=(5, 95)
The lower and upper percentile used to create the extreme values for the grid.
Must be in [0, 100].
n_bootstraps : int, default=20
Number of bootstrap samples to generate if `confidence_interval` is True.
n_jobs : int, default=-1
Number of jobs to run in parallel during bootstrapping. `-1` means using all processors.
cmap : str, default="viridis"
Matplotlib colormap used for the 2D mesh plot.
**kwargs
Expand Down Expand Up @@ -733,9 +768,11 @@ def plot(
feature_idx=features,
feature_type=feature_type,
grid_resolution=grid_resolution,
percentiles=percentiles,
confidence_interval=confidence_interval,
confidence_level=confidence_level,
percentiles=percentiles,
n_bootstraps=n_bootstraps,
n_jobs=n_jobs,
)
result = {
"ale": ale,
Expand Down
47 changes: 46 additions & 1 deletion test/visualization/test_accumulated_local_effects.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@
import matplotlib.pyplot as plt
import numpy as np
import pytest
from joblib import parallel_backend
from sklearn.base import BaseEstimator
from sklearn.ensemble import HistGradientBoostingRegressor
from sklearn.exceptions import NotFittedError
Expand Down Expand Up @@ -246,11 +247,26 @@ def test_compute_ale_1d_continuous(ale_test_data):
assert len(grid_values_ci) <= grid_resolution + 1
assert len(ale_err_ci) == len(grid_values_ci)

# To test _bootstrap()
with parallel_backend("threading"):
result_bootstrap = compute_ale_1d(
model,
X,
feature_idx=important_features[0],
feature_type="continuous",
confidence_interval=True,
n_bootstraps=2,
)


def test_compute_ale_1d_continuous_error(ale_test_data):
"""Check the raised error if the continuous feature does not have enough unique quantiles."""
data = ale_test_data["continuous"]
X, model = data["X"], data["model"]
X, model, important_features = (
data["X"],
data["model"],
data["important_features"],
)

X_const = X.copy()
X_const[:, 0] = 7
Expand All @@ -273,6 +289,15 @@ def test_compute_ale_1d_continuous_error(ale_test_data):
feature_type="invalid_feature_type",
grid_resolution=10,
)
with pytest.raises(ValueError, match="must be strictly greater than"):
compute_ale_1d(
model,
X,
feature_idx=important_features[0],
feature_type="continuous",
confidence_interval=True,
n_bootstraps=0,
)


def test_compute_ale_1d_discrete(ale_test_data):
Expand Down Expand Up @@ -310,6 +335,17 @@ def test_compute_ale_1d_discrete(ale_test_data):
assert len(grid_values_ci) <= 3
assert len(ale_err_ci) == len(grid_values_ci)

# To test _bootstrap()
with parallel_backend("threading"):
result_bootstrap = compute_ale_1d(
model,
X,
feature_idx=0,
feature_type="categorical",
confidence_interval=True,
n_bootstraps=2,
)


def test_compute_ale_1d_discrete_error(ale_test_data):
"""Check the raised error if the discrete feature does not have enough unique values."""
Expand All @@ -323,6 +359,15 @@ def test_compute_ale_1d_discrete_error(ale_test_data):
compute_ale_1d(
model, X_const, feature_idx=0, feature_type="categorical"
)
with pytest.raises(ValueError, match="must be strictly greater than"):
compute_ale_1d(
model,
X,
feature_idx=0,
feature_type="categorical",
confidence_interval=True,
n_bootstraps=0,
)

with pytest.raises(
ValueError, match="'percentiles' must be a tuple of 2 floats"
Expand Down
Loading