Hands-on lab — Day 178: Interpreting Models: Importances and SHAP

import numpy as np
import itertools
from sklearn.base import clone

def compute_permutation_importance(estimator, X, y, metric_fn, n_repeats=5, random_state=42):
    """
    Compute out-of-sample Permutation Feature Importance.
    metric_fn: function(y_true, y_pred) -> score (higher is better)
    """
    rng = np.random.default_rng(random_state)
    X = np.asarray(X, dtype=float)
    y = np.asarray(y, dtype=float)
    
    baseline_pred = estimator.predict(X)
    baseline_score = metric_fn(y, baseline_pred)
    
    n_features = X.shape[1]
    importances = np.zeros((n_features, n_repeats))
    
    for j in range(n_features):
        for r in range(n_repeats):
            X_perm = X.copy()
            perm_indices = rng.permutation(len(X))
            X_perm[:, j] = X_perm[perm_indices, j]
            
            perm_pred = estimator.predict(X_perm)
            perm_score = metric_fn(y, perm_pred)
            importances[j, r] = baseline_score - perm_score
            
    return {
        "importances_mean": np.mean(importances, axis=1).tolist(),
        "importances_std": np.std(importances, axis=1).tolist(),
        "baseline_score": float(baseline_score)
    }

def compute_exact_shapley_values(predict_fn, x_instance, background_data):
    """
    Compute exact Shapley values for an instance across all 2^D feature coalitions.
    predict_fn: function(X) -> 1D array of predictions
    """
    x_instance = np.asarray(x_instance, dtype=float).ravel()
    background_data = np.asarray(background_data, dtype=float)
    d = len(x_instance)
    n_bg = len(background_data)
    
    # Base expected value across background dataset
    base_val = float(np.mean(predict_fn(background_data)))
    
    def evaluate_coalition(subset):
        # Construct synthetic background matrix with features in subset replaced by x_instance
        if len(subset) == 0:
            return base_val
        X_eval = background_data.copy()
        for idx in subset:
            X_eval[:, idx] = x_instance[idx]
        return float(np.mean(predict_fn(X_eval)))
    
    shapley_values = np.zeros(d)
    all_indices = set(range(d))
    
    import math
    for i in range(d):
        other_indices = list(all_indices - {i})
        phi_i = 0.0
        
        # Iterate over all subsets of other indices
        for s_len in range(d):
            subsets = list(itertools.combinations(other_indices, s_len))
            weight = (math.factorial(s_len) * math.factorial(d - s_len - 1)) / math.factorial(d)
            for S in subsets:
                v_with = evaluate_coalition(set(S) | {i})
                v_without = evaluate_coalition(set(S))
                phi_i += weight * (v_with - v_without)
                
        shapley_values[i] = phi_i
        
    instance_pred = float(predict_fn(x_instance[np.newaxis, :])[0])
    
    return {
        "base_value": base_val,
        "instance_prediction": instance_pred,
        "shapley_values": shapley_values.tolist(),
        "efficiency_check_diff": float(abs((base_val + np.sum(shapley_values)) - instance_pred))
    }

def compute_partial_dependence_1d(estimator, X, feature_idx, grid_resolution=20):
    """
    Compute 1D Partial Dependence curve for feature_idx.
    """
    X = np.asarray(X, dtype=float)
    feat_vals = X[:, feature_idx]
    grid = np.linspace(np.min(feat_vals), np.max(feat_vals), grid_resolution)
    
    pdp_values = []
    for val in grid:
        X_pdp = X.copy()
        X_pdp[:, feature_idx] = val
        preds = estimator.predict(X_pdp)
        pdp_values.append(float(np.mean(preds)))
        
    return {
        "grid": grid.tolist(),
        "pdp_values": pdp_values
    }

import pytest
import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score
from shap_lib import (
    compute_permutation_importance,
    compute_exact_shapley_values,
    compute_partial_dependence_1d
)

def test_permutation_importance():
    X = np.random.normal(size=(100, 3))
    # Feature 0 is primary driver, Feature 2 is pure noise
    y = 5.0 * X[:, 0] + 0.5 * X[:, 1] + np.random.normal(scale=0.1, size=100)
    
    model = LinearRegression().fit(X, y)
    res = compute_permutation_importance(model, X, y, r2_score, n_repeats=3)
    
    # Feature 0 must have strictly higher importance than Feature 2
    assert res["importances_mean"][0] > res["importances_mean"][2]
    assert res["importances_mean"][0] > 0.50

def test_exact_shapley_efficiency_axiom():
    # True linear model: y = 2*x0 + 3*x1 + 10
    X_bg = np.array([
        [1.0, 2.0],
        [2.0, 4.0],
        [3.0, 6.0],
        [4.0, 8.0]
    ])
    def predict_fn(X):
        return 2.0 * X[:, 0] + 3.0 * X[:, 1] + 10.0
        
    x_test = np.array([5.0, 10.0])
    res = compute_exact_shapley_values(predict_fn, x_test, X_bg)
    
    # Efficiency axiom: base_value + sum(shapley_values) == instance_prediction
    assert res["efficiency_check_diff"] < 1e-5
    # Feature 1 should have higher attribution than Feature 0
    assert res["shapley_values"][1] > res["shapley_values"][0]

def test_partial_dependence_1d():
    X = np.random.uniform(0, 10, size=(50, 2))
    y = 3.0 * X[:, 0] + 2.0
    model = LinearRegression().fit(X, y)
    
    pdp = compute_partial_dependence_1d(model, X, feature_idx=0, grid_resolution=10)
    assert len(pdp["grid"]) == 10
    assert len(pdp["pdp_values"]) == 10
    # PDP should be strictly monotonically increasing for feature 0
    assert pdp["pdp_values"][-1] > pdp["pdp_values"][0]

lesson_id: D178
day: 178
kind: applied-ml-interpretability
languages:
  - python
setup_commands:
  - python3 -m venv .venv
  - .venv/bin/pip install -r requirements/requirements.txt
run_commands:
  - .venv/bin/python examples/interpreting_models_importances_and_shap_lib.py
test_commands:
  - ./tests/run_tests.sh
cleanup_commands:
  - rm -rf .venv __pycache__ .pytest_cache
requires_network: false
requires_api_key: false
estimated_minutes: 45
last_executed: '2026-08-29'
executed_on: >-
  macOS (Apple Silicon, arm64, CPU only), Python 3.14.0, scikit-learn 1.9.0, pytest 9.1.1 -- bash tests/run_tests.sh -> 4 checks, 0 failure(s), exit 0. Verified Day 178 implementation.

numpy>=1.24.0
scipy>=1.10.0
scikit-learn>=1.3.0
pytest>=7.4.0

import numpy as np

def compute_permutation_importance(estimator, X, y, metric_fn, n_repeats=5, random_state=42):
    # TODO: Implement out-of-sample Permutation Importance
    pass

def compute_exact_shapley_values(predict_fn, x_instance, background_data):
    # TODO: Implement exact Shapley values across 2^D feature coalitions
    pass

def compute_partial_dependence_1d(estimator, X, feature_idx, grid_resolution=20):
    # TODO: Implement Partial Dependence curve computation
    pass

import pytest
import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score
from shap_lib import (
    compute_permutation_importance,
    compute_exact_shapley_values,
    compute_partial_dependence_1d
)

def test_permutation_importance():
    X = np.random.normal(size=(100, 3))
    # Feature 0 is primary driver, Feature 2 is pure noise
    y = 5.0 * X[:, 0] + 0.5 * X[:, 1] + np.random.normal(scale=0.1, size=100)
    
    model = LinearRegression().fit(X, y)
    res = compute_permutation_importance(model, X, y, r2_score, n_repeats=3)
    
    # Feature 0 must have strictly higher importance than Feature 2
    assert res["importances_mean"][0] > res["importances_mean"][2]
    assert res["importances_mean"][0] > 0.50

def test_exact_shapley_efficiency_axiom():
    # True linear model: y = 2*x0 + 3*x1 + 10
    X_bg = np.array([
        [1.0, 2.0],
        [2.0, 4.0],
        [3.0, 6.0],
        [4.0, 8.0]
    ])
    def predict_fn(X):
        return 2.0 * X[:, 0] + 3.0 * X[:, 1] + 10.0
        
    x_test = np.array([5.0, 10.0])
    res = compute_exact_shapley_values(predict_fn, x_test, X_bg)
    
    # Efficiency axiom: base_value + sum(shapley_values) == instance_prediction
    assert res["efficiency_check_diff"] < 1e-5
    # Feature 1 should have higher attribution than Feature 0
    assert res["shapley_values"][1] > res["shapley_values"][0]

def test_partial_dependence_1d():
    X = np.random.uniform(0, 10, size=(50, 2))
    y = 3.0 * X[:, 0] + 2.0
    model = LinearRegression().fit(X, y)
    
    pdp = compute_partial_dependence_1d(model, X, feature_idx=0, grid_resolution=10)
    assert len(pdp["grid"]) == 10
    assert len(pdp["pdp_values"]) == 10
    # PDP should be strictly monotonically increasing for feature 0
    assert pdp["pdp_values"][-1] > pdp["pdp_values"][0]

#!/bin/bash
set -e
pytest tests/ -v