Source code for SALib.analyze.pawn

from typing import Dict

import numpy as np
from scipy.stats import ks_2samp

from . import common_args
from ..util import read_param_file, ResultDict, extract_group_names, _check_groups

[docs] def analyze( problem: Dict, X: np.ndarray, Y: np.ndarray, S: int = 10, print_to_console: bool = False, seed: int = None, ): """Performs PAWN sensitivity analysis. The PAWN method [1] is a moment-independent approach to Global Sensitivity Analysis (GSA). It is described as producing robust results at relatively low sample sizes (see [2]) for the purpose of factor ranking and screening. The distribution of model outputs is examined rather than their variation as is typical in other common GSA approaches. The PAWN method further distinguishes itself from other moment-independent approaches by characterizing outputs by their cumulative distribution function (CDF) as opposed to their probability distribution function. As the CDF for a given random variable is typically normally distributed, PAWN can be more appropriately applied when outputs are highly-skewed or multi-modal, for which variance-based methods may produce unreliable results. PAWN characterizes the relationship between inputs and outputs by quantifying the variation in the output distributions after conditioning an input. A factor is deemed non-influential if distributions coincide at all ``S`` conditioning intervals. The Kolmogorov-Smirnov statistic is used as a measure of distance between the distributions. This implementation reports the PAWN index at the min, mean, median, and max across the slides/conditioning intervals as well as the coefficient of variation (``CV``). The median value is the typically reported value. As the ``CV`` is (standard deviation / mean), it indicates the level of variability across the slides, with values closer to zero indicating lower variation. Notes ----- Compatible with: all samplers This implementation ignores all NaNs. When applied to grouped factors, the analysis is conducted on each factor individually, and the mean of their results are reported. Examples -------- >>> X = latin.sample(problem, 1000) >>> Y = Ishigami.evaluate(X) >>> Si = pawn.analyze(problem, X, Y, S=10, print_to_console=False) Parameters ---------- problem : dict The problem definition X : numpy.array A NumPy array containing the model inputs Y : numpy.array A NumPy array containing the model outputs S : int Number of slides; the conditioning intervals (default 10) print_to_console : bool Print results directly to console (default False) seed : int Seed value to ensure deterministic results References ---------- 1. Pianosi, F., Wagener, T., 2015. A simple and efficient method for global sensitivity analysis based on cumulative distribution functions. Environmental Modelling & Software 67, 1-11. 2. Pianosi, F., Wagener, T., 2018. Distribution-based sensitivity analysis from a generic input-output sample. Environmental Modelling & Software 108, 197-207. 3. Baroni, G., Francke, T., 2020. An effective strategy for combining variance- and distribution-based global sensitivity analysis. Environmental Modelling & Software, 134, 104851. 4. Baroni, G., Francke, T., 2020. GSA-cvd Combining variance- and distribution-based global sensitivity analysis """ if seed: np.random.seed(seed) D = problem["num_vars"] var_names, _ = extract_group_names(problem) results = np.full((D, 5), np.nan) temp_pawn = np.full((S, D), np.nan) step = 1 / S for d_i in range(D): seq = np.arange(0, 1 + step, step) X_di = X[:, d_i] X_q = np.nanquantile(X_di, seq) for s in range(S): Y_sel = Y[(X_di >= X_q[s]) & (X_di < X_q[s + 1])] if len(Y_sel) == 0: # no available samples continue # KD value # Function returns a KS object which holds the KS statistic # and p-value # Note from scipy documentation: # if the K-S statistic is small or the p-value is high, then # we cannot reject the hypothesis that the distributions of # the two samples are the same. ks = ks_2samp(Y_sel, Y) temp_pawn[s, d_i] = ks.statistic p_ind = temp_pawn[:, d_i] mins = np.nanmin(p_ind) mean = np.nanmean(p_ind) med = np.nanmedian(p_ind) maxs = np.nanmax(p_ind) cv = np.nanstd(p_ind) / mean results[d_i, :] = [mins, mean, med, maxs, cv] groups = _check_groups(problem) if groups: unique_grps, n_groups = extract_group_names(problem) tmp = np.full((n_groups, 5), np.nan) # Take the mean of effects from parameters that are grouped together unique_grps = np.array(unique_grps) for grp_id, grp in enumerate(unique_grps): tmp[grp_id, :] = np.mean(results[unique_grps == grp, :], axis=0) results = tmp tmp = None Si = ResultDict( [ ("minimum", results[:, 0]), ("mean", results[:, 1]), ("median", results[:, 2]), ("maximum", results[:, 3]), ("CV", results[:, 4]), ] ) Si["names"] = var_names if print_to_console: print(Si.to_df()) return Si
[docs] def cli_parse(parser): parser.add_argument( "-X", "--model-input-file", type=str, required=True, help="Model input file" ) parser.add_argument( "-S", "--slices", type=int, required=False, help="Number of slices to take" ) return parser
[docs] def cli_action(args): problem = read_param_file(args.paramfile) X = np.loadtxt(args.model_input_file, delimiter=args.delimiter) Y = np.loadtxt( args.model_output_file, delimiter=args.delimiter, usecols=(args.column,) ) analyze(problem, X, Y, S=args.slices, print_to_console=True, seed=args.seed)
if __name__ == "__main__": common_args.run_cli(cli_parse, cli_action)