from types import MethodType
import numpy as np
from scipy.stats import norm
from . import common_args
from ..util import ResultDict, read_param_file
from ..plotting.shapley import plot as shapley_plot
[docs]
def analyze(
problem: dict,
X: np.ndarray,
Y: np.ndarray,
conf_level: float = 0.95,
print_to_console: bool = False,
) -> ResultDict:
"""Estimate Shapley effects with Goda's Monte Carlo algorithm.
Returns Shapley effects in output-variance units (``shapley``, summing to
the overall output variance) alongside normalized shares that sum to one
(``shapley_normalized``), each with a normal-approximation confidence
interval half-width (``shapley_conf``, ``shapley_normalized_conf``).
Notes
-----
Compatible with:
Compatible with :func:`SALib.sample.shapley.sample`.
The estimator assumes mutually independent inputs and does not support
grouped parameters.
Parameters
----------
problem : dict
The problem definition.
X : numpy.ndarray
Model inputs generated by :func:`SALib.sample.shapley.sample`.
Y : numpy.ndarray
One-dimensional array of model outputs.
conf_level : float, default=0.95
Confidence level used for normal-approximation interval half-widths.
print_to_console : bool, default=False
Print results directly to the console.
Returns
-------
ResultDict
Raw and normalized Shapley effect estimates, their confidence
interval half-widths, and factor names.
References
----------
1. Goda, T. (2021). A simple algorithm for global sensitivity analysis
with Shapley effects. Reliability Engineering & System Safety, 213,
107702. https://doi.org/10.1016/j.ress.2021.107702
"""
_validate_ungrouped(problem)
X, Y, num_trajectories = _validate_inputs(problem, X, Y, conf_level)
num_vars = problem["num_vars"]
trajectory_size = num_vars + 1
trajectories = X.reshape(num_trajectories, trajectory_size, num_vars)
outputs = Y.reshape(num_trajectories, trajectory_size)
# Each trajectory step changes exactly one factor; recover which one by
# diffing consecutive rows, then scatter the step-ordered contributions
# into factor-ordered columns with a single vectorized reorder instead of
# a per-trajectory, per-step Python loop.
changed = trajectories[:, 1:, :] != trajectories[:, :-1, :]
change_counts = changed.sum(axis=2)
mismatched_steps = change_counts != 1
if np.any(mismatched_steps):
# Identify location of mismatch
trajectory_index, step = (
int(i)
for i in np.unravel_index(np.argmax(mismatched_steps), change_counts.shape)
)
raise ValueError(
"Each Shapley trajectory step must change exactly one factor; "
f"trajectory {trajectory_index}, step {step + 1} changed "
f"{int(change_counts[trajectory_index, step])}."
)
factors = changed.argmax(axis=2)
incomplete = np.any(np.sort(factors, axis=1) != np.arange(num_vars), axis=1)
if np.any(incomplete):
trajectory_index = int(np.argmax(incomplete))
raise ValueError(
"Each Shapley trajectory must change every factor exactly once; "
f"trajectory {trajectory_index} does not."
)
base_output = outputs[:, :1]
before = outputs[:, :-1]
after = outputs[:, 1:]
step_contributions = (base_output - (before + after) / 2.0) * (before - after)
contributions = np.empty_like(step_contributions)
np.put_along_axis(contributions, factors, step_contributions, axis=1)
effects = contributions.mean(axis=0)
estimator_variance = contributions.var(axis=0, ddof=1) / num_trajectories
z_score = norm.ppf(0.5 + conf_level / 2.0)
confidence = z_score * np.sqrt(estimator_variance)
normalized, normalized_confidence = _normalize_with_confidence(
contributions, effects, z_score
)
result = ResultDict(
shapley_normalized=normalized,
shapley_normalized_conf=normalized_confidence,
shapley=effects,
shapley_conf=confidence,
names=problem["names"],
)
result.plot = MethodType(shapley_plot, result)
if print_to_console:
print(result.to_df())
return result
def _normalize_with_confidence(contributions, effects, z_score):
"""Estimate normalized Shapley shares and their confidence half-widths.
Each trajectory's per-factor contributions sum to that trajectory's own
estimate of ``Var[Y]``, so the normalized share of a factor is a ratio of two
correlated trajectory-averaged quantities: the mean contribution to that
factor, and the mean row total. Its sampling variance is approximated
with the standard delta-method formula for a ratio of means (e.g.
Cochran, 1977, on ratio estimators), reusing the per-trajectory
contributions already computed for the raw ``shapley_conf`` interval.
"""
num_trajectories = contributions.shape[0]
row_sums = contributions.sum(axis=1)
total = row_sums.mean()
if total == 0.0:
undefined = np.full_like(effects, np.nan)
return undefined, undefined.copy()
normalized = effects / total
centered = contributions - effects
centered_totals = row_sums - total
var_contributions = np.sum(centered**2, axis=0) / (num_trajectories - 1)
covar_with_total = np.sum(centered * centered_totals[:, None], axis=0) / (
num_trajectories - 1
)
var_total = np.sum(centered_totals**2) / (num_trajectories - 1)
ratio_variance = (
var_contributions
- 2.0 * normalized * covar_with_total
+ normalized**2 * var_total
) / (num_trajectories * total**2)
# Guard against tiny negative values from floating-point round-off.
ratio_variance = np.clip(ratio_variance, a_min=0.0, a_max=None)
return normalized, z_score * np.sqrt(ratio_variance)
def _validate_inputs(problem, X, Y, conf_level):
if not 0 < conf_level < 1:
raise ValueError("Confidence level must be between 0 and 1.")
num_vars = problem["num_vars"]
if num_vars <= 0:
raise ValueError("The problem must contain at least one variable.")
X = np.asarray(X)
Y = np.asarray(Y)
if X.ndim != 2 or X.shape[1] != num_vars:
raise ValueError(f"X must have shape (N * (D + 1), {num_vars}).")
if Y.ndim != 1:
raise ValueError("Y must be a one-dimensional array.")
if X.shape[0] != Y.size:
raise ValueError("X and Y must contain the same number of rows.")
if Y.size % (num_vars + 1):
raise ValueError("The number of model outputs must be divisible by D + 1.")
num_trajectories = Y.size // (num_vars + 1)
if num_trajectories < 2:
raise ValueError("At least two Shapley trajectories are required.")
if not np.all(np.isfinite(X)) or not np.all(np.isfinite(Y)):
raise ValueError("X and Y must contain only finite values.")
return X, Y, num_trajectories
def _validate_ungrouped(problem: dict) -> None:
groups = problem.get("groups")
if groups and list(groups) != list(problem.get("names", [])):
raise ValueError("Goda's Shapley estimator does not support groups.")
[docs]
def cli_parse(parser):
parser.add_argument(
"-X",
"--model-input-file",
type=str,
required=True,
help="Model input file",
)
parser.add_argument(
"--conf-level",
type=float,
required=False,
default=0.95,
help="Confidence interval level",
)
return parser
[docs]
def cli_action(args):
"""Analyze Shapley trajectories from command-line arguments."""
problem = read_param_file(args.paramfile)
X = np.loadtxt(args.model_input_file, delimiter=args.delimiter, ndmin=2)
Y = np.loadtxt(
args.model_output_file,
delimiter=args.delimiter,
usecols=(args.column,),
ndmin=1,
)
analyze(
problem,
X,
Y,
conf_level=args.conf_level,
print_to_console=True,
)
if __name__ == "__main__":
common_args.run_cli(cli_parse, cli_action)
