Updated

Author: RAMitchell

algorithms.py

@@ -1,7 +1,7 @@
 import numpy as np
 from numba import njit
 import sobol_seq
-import kernel_herding
+import kernel_methods
 from itertools import count
 from math import cos, gamma, pi, sin, sqrt
 from typing import Callable, Iterator, List
@@ -84,14 +84,16 @@ def owen_complement(X_background, X_foreground, predict_function, n_samples):
 
 
 @njit
-def _accumulate_samples_castro(phi, predictions, j):
+def _accumulate_samples_castro(phi, predictions, j, weights=None):
+    if weights == None:
+        weights = np.full(predictions.shape[1], 1 / predictions.shape[1])
     for foreground_idx in range(predictions.shape[0]):
         for sample_idx in range(predictions.shape[1]):
             phi[foreground_idx, j[sample_idx]] += predictions[foreground_idx][
-                                                      sample_idx] / predictions.shape[1]
+                                                      sample_idx] * weights[sample_idx]
 
 
-def estimate_shap_given_permutations(X_background, X_foreground, predict_function, p):
+def estimate_shap_given_permutations(X_background, X_foreground, predict_function, p, weights=None):
     n_features = X_background.shape[1]
     phi = np.zeros((X_foreground.shape[0], n_features))
     n_permutations = p.shape[0]
@@ -105,7 +107,7 @@ def estimate_shap_given_permutations(X_background, X_foreground, predict_functio
         predictions = (pred_on - pred_off).reshape(
             (X_foreground.shape[0], mask.shape[0], X_background.shape[0]))
         predictions = np.mean(predictions, axis=2)
-        _accumulate_samples_castro(phi, predictions, j)
+        _accumulate_samples_castro(phi, predictions, j, weights)
         pred_off = pred_on
 
     return phi
@@ -123,6 +125,28 @@ def monte_carlo(X_background, X_foreground, predict_function, n_samples):
     return estimate_shap_given_permutations(X_background, X_foreground, predict_function, p)
 
 
+def monte_carlo_weighted(X_background, X_foreground, predict_function, n_samples):
+    n_features = X_background.shape[1]
+    assert n_samples % (n_features + 1) == 0
+    # castro is allowed to take 2 * more samples than owen as it reuses predictions
+    samples_per_feature = 2 * (n_samples // (n_features + 1))
+    p = np.zeros((samples_per_feature, n_features), dtype=np.int64)
+    for i in range(samples_per_feature):
+        p[i] = np.random.permutation(n_features)
+    weights = kernel_methods.compute_bayesian_weights(p, kernel_methods.kt_kernel)
+    return estimate_shap_given_permutations(X_background, X_foreground, predict_function, p,
+                                            weights)
+
+
+def sbq(X_background, X_foreground, predict_function, n_samples):
+    n_features = X_background.shape[1]
+    assert n_samples % (n_features + 1) == 0
+    samples_per_feature = 2 * (n_samples // (n_features + 1))
+    p, w = kernel_methods.sequential_bayesian_quadrature(samples_per_feature, n_features)
+    return estimate_shap_given_permutations(X_background, X_foreground, predict_function, p,
+                                            w)
+
+
 def monte_carlo_antithetic(X_background, X_foreground, predict_function, n_samples):
     n_features = X_background.shape[1]
 
@@ -157,6 +181,7 @@ def sobol_sphere_permutations(n_samples, n_features):
 
     return np.argsort(sobol, axis=1)
 
+
 # sample with l ones and i off
 def draw_castro_stratified_samples(n_samples, n_features, i, l):
     mask = np.zeros((n_samples, n_features - 1), dtype=bool)
@@ -272,7 +297,7 @@ def kt_herding(X_background, X_foreground, predict_function, n_samples):
     assert n_samples % (n_features + 1) == 0
     # castro is allowed to take 2 * more samples than owen as it reuses predictions
     samples_per_feature = 2 * (n_samples // (n_features + 1))
-    p = kernel_herding.kt_herding_permutations(samples_per_feature, n_features)
+    p = kernel_methods.kt_herding_permutations(samples_per_feature, n_features)
     return estimate_shap_given_permutations(X_background, X_foreground, predict_function, p)
 
 
@@ -330,6 +355,16 @@ def orthogonal(X_background, X_foreground, predict_function, n_samples):
     return estimate_shap_given_permutations(X_background, X_foreground, predict_function, p)
 
 
+def orthogonal_weighted(X_background, X_foreground, predict_function, n_samples):
+    n_features = X_background.shape[1]
+    assert n_samples % (2 * (n_features + 1)) == 0
+    # castro is allowed to take 2 * more samples than owen as it reuses predictions
+    samples_per_feature = 2 * (n_samples // (n_features + 1))
+    p = _orthogonal_permutations(samples_per_feature, n_features)
+    w = kernel_methods.compute_bayesian_weights(p, kernel_methods.kt_kernel)
+    return estimate_shap_given_permutations(X_background, X_foreground, predict_function, p, w)
+
+
 def _int_sin_m(x: float, m: int) -> float:
     """Computes the integral of sin^m(t) dt from 0 to x recursively"""
     if m == 0:
@@ -438,6 +473,10 @@ def fibonacci(X_background, X_foreground, predict_function, n_samples):
 def min_sample_size(alg, n_features):
     if alg == monte_carlo:
         return n_features + 1
+    elif alg == monte_carlo_weighted:
+        return n_features + 1
+    elif alg == sbq:
+        return n_features + 1
     elif alg == qmc_sobol:
         return n_features + 1
     elif alg == fibonacci:
@@ -450,6 +489,8 @@ def min_sample_size(alg, n_features):
         return 2 * (n_features + 1)
     elif alg == orthogonal:
         return 2 * (n_features + 1)
+    elif alg == orthogonal_weighted:
+        return 2 * (n_features + 1)
     elif alg == owen or alg == owen_complement:
         return n_features * 4
     elif alg == castro_stratified:

experiments.py

@@ -21,23 +21,21 @@ def plot_experiments():
     repeats = 25
     foreground_examples = 10
     background_examples = 100
-    max_evals = 100000
+    max_evals = 5000
     datasets_set = {
         "make_regression": datasets.get_regression(foreground_examples, background_examples),
         "cal_housing": datasets.get_cal_housing(foreground_examples, background_examples),
         "adult": datasets.get_adult(foreground_examples, background_examples),
         "breast_cancer": datasets.get_breast_cancer(foreground_examples, background_examples),
     }
     algorithms_set = {
-        "Castro": algorithms.monte_carlo,
-        "Castro-Complement": algorithms.monte_carlo_antithetic,
-        # "Castro-LHS": algorithms.castro_lhs,
-    }
-    algorithms_set = {
-        # "Castro": algorithms.castro,
-        "Castro-Orthogonal": algorithms.orthogonal,
-        "Castro-Complement": algorithms.monte_carlo_antithetic,
-        "Fibonacci": algorithms.fibonacci,
+        # "MC": algorithms.monte_carlo,
+        # "Bayesian-MC": algorithms.monte_carlo_weighted,
+        "SBQ": algorithms.sbq,
+        # "MC-Orthogonal-Bayesian": algorithms.orthogonal_weighted,
+        "MC-Orthogonal": algorithms.orthogonal,
+        # "Castro-Complement": algorithms.monte_carlo_antithetic,
+        # "Fibonacci": algorithms.fibonacci,
         # "Castro-ControlVariate": algorithms.castro_control_variate,
         # "Castro-QMC": algorithms.castro_qmc,
         # "KT-Herding": algorithms.kt_herding,