wehub-resource-sync
228 lines
9.8 KiB
Python
228 lines
9.8 KiB
Python
import logging
|
|
from typing import Any
|
|
|
|
import numpy as np
|
|
import numpy.typing as npt
|
|
import pandas as pd
|
|
|
|
from .._explanation import Explanation
|
|
from ..utils._exceptions import ExplainerError
|
|
from ..utils._legacy import convert_to_instance, match_instance_to_data
|
|
from ._kernel import KernelExplainer
|
|
|
|
log = logging.getLogger("shap")
|
|
|
|
|
|
class SamplingExplainer(KernelExplainer):
|
|
"""Computes SHAP values using an extension of the Shapley sampling values explanation method
|
|
(also known as IME).
|
|
|
|
SamplingExplainer computes SHAP values under the assumption of feature independence and is an
|
|
extension of the algorithm proposed in "An Efficient Explanation of Individual Classifications
|
|
using Game Theory", Erik Strumbelj, Igor Kononenko, JMLR 2010. It is a good alternative to
|
|
KernelExplainer when you want to use a large background set (as opposed to a single reference
|
|
value for example).
|
|
|
|
Parameters
|
|
----------
|
|
model : function
|
|
User supplied function that takes a matrix of samples (# samples x # features) and
|
|
computes the output of the model for those samples. The output can be a vector
|
|
(# samples) or a matrix (# samples x # model outputs).
|
|
|
|
data : numpy.array or pandas.DataFrame
|
|
The background dataset to use for integrating out features. To determine the impact
|
|
of a feature, that feature is set to "missing" and the change in the model output
|
|
is observed. Since most models aren't designed to handle arbitrary missing data at test
|
|
time, we simulate "missing" by replacing the feature with the values it takes in the
|
|
background dataset. So if the background dataset is a simple sample of all zeros, then
|
|
we would approximate a feature being missing by setting it to zero. Unlike the
|
|
KernelExplainer, this data can be the whole training set, even if that is a large set. This
|
|
is because SamplingExplainer only samples from this background dataset.
|
|
|
|
"""
|
|
|
|
def __init__(
|
|
self,
|
|
model: Any,
|
|
data: npt.NDArray[Any] | pd.DataFrame,
|
|
**kwargs: Any,
|
|
) -> None:
|
|
# silence warning about large datasets
|
|
level = log.level
|
|
log.setLevel(logging.ERROR)
|
|
super().__init__(model, data, **kwargs)
|
|
log.setLevel(level)
|
|
|
|
if str(self.link) != "identity":
|
|
emsg = f"SamplingExplainer only supports the identity link, not {self.link}"
|
|
raise ValueError(emsg)
|
|
|
|
def __call__( # type: ignore[override]
|
|
self,
|
|
X: npt.NDArray[Any] | pd.DataFrame,
|
|
y: Any = None,
|
|
nsamples: int = 2000,
|
|
) -> Explanation:
|
|
if isinstance(X, pd.DataFrame):
|
|
feature_names = list(X.columns)
|
|
X = X.values
|
|
else:
|
|
feature_names = None # we can make self.feature_names from background data eventually if we have it
|
|
|
|
v = self.shap_values(X, nsamples=nsamples)
|
|
if isinstance(v, list):
|
|
v = np.stack(v, axis=-1) # put outputs at the end
|
|
e = Explanation(v, self.expected_value, X, feature_names=feature_names)
|
|
return e
|
|
|
|
def explain(
|
|
self,
|
|
incoming_instance: Any,
|
|
**kwargs: Any,
|
|
) -> npt.NDArray[np.floating[Any]]:
|
|
# convert incoming input to a standardized iml object
|
|
instance = convert_to_instance(incoming_instance)
|
|
match_instance_to_data(instance, self.data)
|
|
|
|
if len(self.data.groups) != self.P: # type: ignore[arg-type]
|
|
emsg = "SamplingExplainer does not support feature groups!"
|
|
raise ExplainerError(emsg)
|
|
|
|
# find the feature groups we will test. If a feature does not change from its
|
|
# current value then we know it doesn't impact the model
|
|
self.varyingInds = self.varying_groups(instance.x)
|
|
# self.varyingFeatureGroups = [self.data.groups[i] for i in self.varyingInds]
|
|
self.M = len(self.varyingInds)
|
|
|
|
# find f(x)
|
|
if self.keep_index:
|
|
model_out = self.model.f(instance.convert_to_df())
|
|
else:
|
|
model_out = self.model.f(instance.x)
|
|
if isinstance(model_out, pd.DataFrame):
|
|
model_out = model_out.values[0]
|
|
elif isinstance(model_out, pd.Series):
|
|
model_out = model_out.values
|
|
self.fx = model_out[0]
|
|
|
|
if not self.vector_out:
|
|
self.fx = np.array([self.fx])
|
|
|
|
# if no features vary then there no feature has an effect
|
|
if self.M == 0:
|
|
phi = np.zeros((len(self.data.groups), self.D)) # type: ignore[arg-type]
|
|
phi_var = np.zeros((len(self.data.groups), self.D)) # type: ignore[arg-type]
|
|
|
|
# if only one feature varies then it has all the effect
|
|
elif self.M == 1:
|
|
phi = np.zeros((len(self.data.groups), self.D)) # type: ignore[arg-type]
|
|
phi_var = np.zeros((len(self.data.groups), self.D)) # type: ignore[arg-type]
|
|
diff = self.fx - self.fnull
|
|
for d in range(self.D):
|
|
phi[self.varyingInds[0], d] = diff[d]
|
|
|
|
# if more than one feature varies then we have to do real work
|
|
else:
|
|
# pick a reasonable number of samples if the user didn't specify how many they wanted
|
|
self.nsamples = kwargs.get("nsamples", "auto")
|
|
if self.nsamples == "auto":
|
|
self.nsamples = 1000 * self.M
|
|
|
|
min_samples_per_feature = kwargs.get("min_samples_per_feature", 100)
|
|
round1_samples = self.nsamples
|
|
round2_samples = 0
|
|
if round1_samples > self.M * min_samples_per_feature:
|
|
round2_samples = round1_samples - self.M * min_samples_per_feature
|
|
round1_samples -= round2_samples
|
|
|
|
# divide up the samples among the features for round 1
|
|
nsamples_each1 = np.ones(self.M, dtype=np.int64) * 2 * (round1_samples // (self.M * 2))
|
|
for i in range((round1_samples % (self.M * 2)) // 2):
|
|
nsamples_each1[i] += 2
|
|
|
|
# explain every feature in round 1
|
|
phi = np.zeros((self.P, self.D))
|
|
phi_var = np.zeros((self.P, self.D))
|
|
self.X_masked = np.zeros((nsamples_each1.max() * 2, self.data.data.shape[1]))
|
|
for i, ind in enumerate(self.varyingInds):
|
|
phi[ind, :], phi_var[ind, :] = self.sampling_estimate(
|
|
ind, self.model.f, instance.x, self.data.data, nsamples=nsamples_each1[i]
|
|
)
|
|
|
|
# optimally allocate samples according to the variance
|
|
if phi_var.sum() == 0:
|
|
phi_var += 1 # spread samples uniformally if we found no variability
|
|
phi_var /= phi_var.sum(0)[np.newaxis, :]
|
|
nsamples_each2 = (phi_var[self.varyingInds, :].mean(1) * round2_samples).astype(int)
|
|
for i in range(len(nsamples_each2)):
|
|
if nsamples_each2[i] % 2 == 1:
|
|
nsamples_each2[i] += 1
|
|
for i in range(len(nsamples_each2)):
|
|
if nsamples_each2.sum() > round2_samples:
|
|
nsamples_each2[i] -= 2
|
|
elif nsamples_each2.sum() < round2_samples:
|
|
nsamples_each2[i] += 2
|
|
else:
|
|
break
|
|
|
|
self.X_masked = np.zeros((nsamples_each2.max() * 2, self.data.data.shape[1]))
|
|
for i, ind in enumerate(self.varyingInds):
|
|
if nsamples_each2[i] > 0:
|
|
val, var = self.sampling_estimate(
|
|
ind, self.model.f, instance.x, self.data.data, nsamples=nsamples_each2[i]
|
|
)
|
|
|
|
total_samples = nsamples_each1[i] + nsamples_each2[i]
|
|
phi[ind, :] = (phi[ind, :] * nsamples_each1[i] + val * nsamples_each2[i]) / total_samples
|
|
phi_var[ind, :] = (phi_var[ind, :] * nsamples_each1[i] + var * nsamples_each2[i]) / total_samples
|
|
|
|
# convert from the variance of the differences to the variance of the mean (phi)
|
|
for i, ind in enumerate(self.varyingInds):
|
|
phi_var[ind, :] /= np.sqrt(nsamples_each1[i] + nsamples_each2[i])
|
|
|
|
# correct the sum of the SHAP values to equal the output of the model using a linear
|
|
# regression model with priors of the coefficients equal to the estimated variances for each
|
|
# SHAP value (note that 1e6 is designed to increase the weight of the sample and so closely
|
|
# match the correct sum)
|
|
sum_error = self.fx - phi.sum(0) - self.fnull
|
|
for i in range(self.D):
|
|
# this is a ridge regression with one sample of all ones with sum_error[i] as the label
|
|
# and 1/v as the ridge penalties. This simplified (and stable) form comes from the
|
|
# Sherman-Morrison formula
|
|
v = (phi_var[:, i] / phi_var[:, i].max()) * 1e6
|
|
adj = sum_error[i] * (v - (v * v.sum()) / (1 + v.sum()))
|
|
phi[:, i] += adj
|
|
|
|
if phi.shape[1] == 1:
|
|
phi = phi[:, 0]
|
|
|
|
return phi
|
|
|
|
def sampling_estimate(
|
|
self,
|
|
j: int,
|
|
f: Any,
|
|
x: npt.NDArray[Any],
|
|
X: npt.NDArray[Any],
|
|
nsamples: int = 10,
|
|
) -> tuple[npt.NDArray[np.floating[Any]], npt.NDArray[np.floating[Any]]]:
|
|
X_masked = self.X_masked[: nsamples * 2, :]
|
|
inds = np.arange(X.shape[1])
|
|
|
|
for i in range(nsamples):
|
|
np.random.shuffle(inds)
|
|
pos = np.where(inds == j)[0][0]
|
|
rind = np.random.randint(X.shape[0])
|
|
X_masked[i, :] = x
|
|
X_masked[i, inds[pos + 1 :]] = X[rind, inds[pos + 1 :]]
|
|
X_masked[-(i + 1), :] = x
|
|
X_masked[-(i + 1), inds[pos:]] = X[rind, inds[pos:]]
|
|
|
|
evals = f(X_masked)
|
|
evals_on = evals[:nsamples]
|
|
evals_off = evals[nsamples:][::-1]
|
|
d = evals_on - evals_off
|
|
|
|
return np.mean(d, 0), np.var(d, 0)
|