Comprehensive Global Sensitivity Analysis Tutorial¶

# Import required libraries
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from mcpost import gsa_pipeline, gsa_for_target
from mcpost import plot_sensitivity_metrics

# Quieten sklearn's per-restart optimiser chatter. MCPost's own warnings --
# surrogate collapse, Sobol range violations, length scales on a bound -- are
# NOT suppressed: those are the ones you want to see.
import warnings
from sklearn.exceptions import ConvergenceWarning
warnings.filterwarnings("ignore", category=ConvergenceWarning)

# Set random seed for reproducibility
np.random.seed(42)

print("MCPost comprehensive GSA tutorial loaded!")

# Define the Ishigami function
def ishigami_function(X):
    """
    Ishigami function: f(x1,x2,x3) = sin(x1) + 7*sin(x2)^2 + 0.1*x3^4*sin(x1)
    
    Parameters should be in range [-π, π]
    """
    x1, x2, x3 = X[:, 0], X[:, 1], X[:, 2]
    return np.sin(x1) + 7 * np.sin(x2)**2 + 0.1 * x3**4 * np.sin(x1)

# Generate samples
n_samples = 1200
X_ishigami = np.random.uniform(-np.pi, np.pi, (n_samples, 3))

# Evaluate function
y_ishigami = ishigami_function(X_ishigami)
Y_ishigami = y_ishigami.reshape(-1, 1)

print(f"Generated {n_samples} Ishigami function samples")
print("(A GP is O(n^3) to fit; n=1200 is plenty for a 3-parameter demo.)")
print(f"Output range: [{y_ishigami.min():.3f}, {y_ishigami.max():.3f}]")
print(f"Output std: {y_ishigami.std():.3f}")

# Run comprehensive GSA on Ishigami function
param_names = ["x1", "x2", "x3"]
feature_names = ["ishigami"]

print("Running comprehensive GSA on Ishigami function...")
results_ishigami = gsa_pipeline(
    X_ishigami, Y_ishigami,
    param_names=param_names,
    feature_names=feature_names,
    scaler="minmax",
    enable_sobol=True,
    enable_gp=True,
    enable_perm=True,
    make_pdp=True,
    topk_pdp=3,
    N_sobol=4096,  # Higher for better Sobol accuracy
)

# Display results
ishigami_table = results_ishigami["results"]["ishigami"]["table"]
print("\nIshigami Function - GSA Results:")
print(ishigami_table.round(4))

# Compare with the CLOSED-FORM Sobol indices.
#
# For f = sin(x1) + a*sin(x2)^2 + b*x3^4*sin(x1) with a = 7, b = 0.1 and
# x ~ U(-pi, pi)^3:
#     V1  = 0.5 * (1 + b*pi^4/5)^2       V2 = a^2/8       V3 = 0
#     V13 = 8 * b^2 * pi^8 / 225
#     V   = a^2/8 + b*pi^4/5 + b^2*pi^8/18 + 0.5
# The table columns are S1 (first order) and ST (total effect).
a, b = 7.0, 0.1
V1 = 0.5 * (1 + b * np.pi**4 / 5) ** 2
V2 = a**2 / 8
V13 = 8 * b**2 * np.pi**8 / 225
V = a**2 / 8 + b * np.pi**4 / 5 + b**2 * np.pi**8 / 18 + 0.5
theoretical = {
    "x1": {"S1": V1 / V,  "ST": (V1 + V13) / V},
    "x2": {"S1": V2 / V,  "ST": V2 / V},
    "x3": {"S1": 0.0,     "ST": V13 / V},
}

print("\nComparison with the analytic Sobol indices:")
header = f"{'Param':<6}{'S1 (GP)':>10}{'S1 (exact)':>12}{'err':>9}"
header += f"{'ST (GP)':>10}{'ST (exact)':>12}{'err':>9}"
print(header)
print("-" * len(header))
for param in param_names:
    s1_c = ishigami_table.loc[param, "S1"]
    st_c = ishigami_table.loc[param, "ST"]
    s1_t = theoretical[param]["S1"]
    st_t = theoretical[param]["ST"]
    print(f"{param:<6}{s1_c:>10.4f}{s1_t:>12.4f}{abs(s1_c - s1_t):>9.4f}"
          f"{st_c:>10.4f}{st_t:>12.4f}{abs(st_c - st_t):>9.4f}")

print("\nNote x3: its first-order effect is exactly zero, but its total effect")
print("is not -- x3 acts only through its interaction with x1. A surrogate that")
print("confused S1 with ST could not reproduce both columns.")

# Visualize sensitivity results
# plot_sensitivity_metrics(table, metrics=None, title=..., cmap=...) -- it
# sizes the figure itself, so there is no `figsize` argument. Pass `metrics` to
# restrict the heatmap to comparable columns; mixing S1 (a fraction) with MI
# (unbounded) on one colour scale makes both harder to read.
fig, ax = plot_sensitivity_metrics(
    ishigami_table,
    metrics=["MI", "dCor", "PermMean", "S1", "ST"],
    title="Ishigami Function - Comprehensive Sensitivity Analysis",
)
plt.tight_layout()
plt.show()

print("\nInterpretation:")
print("- x2 shows highest sensitivity across all metrics (expected)")
print("- x1 shows moderate main effect but high total effect (interactions with x3)")
print("- x3 shows low main effect but contributes through interactions")
print("- Different metrics provide complementary insights")

# Compare different scaling options
#
# gsa_for_target returns a (table, extras) PAIR. `extras` carries the fit
# diagnostics that tell you whether the numbers in `table` mean anything.
scaling_options = ["minmax", "standard", None]
scaling_results = {}

print("Comparing scaling options on Ishigami function:")
print("=" * 62)

for scaler in scaling_options:
    table, extras = gsa_for_target(
        X_ishigami, y_ishigami,
        param_names=param_names,
        scaler=scaler,
        enable_perm=True,
        enable_gp=True,
        enable_sobol=True,
        N_sobol=2048,
    )
    scaling_results[str(scaler)] = (table, extras)

    print(f"\nscaler = {scaler!r}")
    print(f"  MI ranking : {table.sort_values('MI', ascending=False).index.tolist()}")
    print(f"  S1 ranking : {table.sort_values('S1', ascending=False).index.tolist()}")
    print(f"  sum(S1)    : {table['S1'].sum():.4f}   max ST: {table['ST'].max():.4f}")
    print(f"  surrogate collapsed : {extras['surrogate_collapsed']}")
    print(f"  log marginal lik.   : {extras['log_marginal_likelihood']:.1f}")
    print(f"  length scales on a bound : {extras['ard_at_bounds']}")

print("\nWhat to look at here:")
print("- 'surrogate collapsed' is the guard that matters. A collapsed GP predicts")
print("  a constant, so every Sobol index comes back 0 -- which reads as")
print("  'no parameter matters' and is indistinguishable from a real result.")
print("- log marginal likelihood compares fits on the SAME data. A fit that is")
print("  hundreds of nats worse than another is stuck, not merely different.")
print("- Rankings should agree across scalers. If they do not, treat the GP-based")
print("  columns (ARD_LS, S1, ST) with suspicion before the screening ones.")

# Compare different kernel types for GP
kernel_options = ["rbf", "matern32", "matern52"]
kernel_results = {}

print("Comparing GP kernel options:")
print("=" * 46)

for kernel in kernel_options:
    table, extras = gsa_for_target(
        X_ishigami, y_ishigami,
        param_names=param_names,
        scaler="minmax",
        kernel_kind=kernel,
        enable_gp=True,
        enable_perm=False,   # Skip for speed
        enable_sobol=False,  # Skip for speed
    )
    kernel_results[kernel] = (table, extras)

    ard_ranking = table.sort_values("1/ARD_LS", ascending=False).index.tolist()
    print(f"\nkernel = {kernel}")
    print(f"  ARD ranking : {ard_ranking}")
    print(f"  1/ARD_LS    : {table['1/ARD_LS'].round(3).to_dict()}")
    print(f"  log marginal likelihood : {extras['log_marginal_likelihood']:.1f}")
    print(f"  length scales on a bound: {extras['ard_at_bounds']}")

print("\nNote: Different kernels encode different smoothness assumptions")
print("- RBF: Infinitely smooth (most common)")
print("- Matern 3/2: Less smooth, more flexible")
print("- Matern 5/2: Moderately smooth")
print("\nA length scale resting ON a bound is where the optimiser stopped, not")
print("where the likelihood peaked, so the 1/ARD_LS entry for that parameter is")
print("a property of the bound rather than a measurement. Compare kernels by")
print("log marginal likelihood, not by which one gives the ranking you expected.")

# Create a function with known characteristics for interpretation
def interpretation_function(X):
    """
    Function designed to demonstrate different sensitivity patterns:
    - x1: Linear effect
    - x2: Quadratic effect  
    - x3: Interaction with x1
    - x4: Threshold effect
    - x5: Noise parameter (minimal effect)
    """
    x1, x2, x3, x4, x5 = X[:, 0], X[:, 1], X[:, 2], X[:, 3], X[:, 4]
    
    linear_term = 2 * x1
    quadratic_term = x2**2
    interaction_term = x1 * x3
    threshold_term = np.where(x4 > 0, x4**2, 0)
    noise_term = 0.01 * x5
    
    return linear_term + quadratic_term + interaction_term + threshold_term + noise_term

# Generate samples
n_samples = 1200
X_interp = np.random.uniform(-1, 1, (n_samples, 5))
y_interp = interpretation_function(X_interp)

param_names_interp = ["x1_linear", "x2_quadratic", "x3_interaction", "x4_threshold", "x5_noise"]

# Work out what to expect, rather than guessing. On U(-1,1) the terms are
# independent, so the variances just add:
v_x1  = 4 * (1 / 3)                          # Var(2*x1)
v_x2  = 1 / 5 - (1 / 3) ** 2                 # Var(x2^2)
v_x13 = (1 / 3) * (1 / 3)                    # Var(x1*x3)  -- an INTERACTION
v_x4  = 0.5 * (1 / 9) - (0.5 * (1 / 5)) ** 2 # Var(threshold term)
v_x5  = (0.01 ** 2) / 3                      # Var(0.01*x5)
V = v_x1 + v_x2 + v_x13 + v_x4 + v_x5

print("Interpretation Function Analysis:")
print(f"Generated {n_samples} samples")
print("\nAnalytic first-order indices (S1):")
for name, v in zip(param_names_interp, [v_x1, v_x2, 0.0, v_x4, v_x5]):
    print(f"  {name:<16}{v / V:.4f}")
print(f"\nx3_interaction has NO main effect at all -- x1*x3 averages to zero over")
print(f"x1 -- so its S1 is 0 while its total effect is {v_x13 / V:.4f}.")
print("Expected S1 ranking: x1 >> x2 > x4 > x5, with x3 at zero.")
print("Note this contradicts the intuition that a 'quadratic effect' must")
print("dominate: 2*x1 has 15x the variance of x2^2 here.")

# Run GSA on interpretation function
results_interp, extras_interp = gsa_for_target(
    X_interp, y_interp,
    param_names=param_names_interp,
    scaler="minmax",
    enable_perm=True,
    enable_gp=True,
    enable_sobol=True,
    N_sobol=2048,
)

print("\nInterpretation Function - GSA Results:")
print(results_interp[["MI", "dCor", "PermMean", "S1", "ST"]].round(4))
print(f"\nsurrogate collapsed: {extras_interp['surrogate_collapsed']}   "
      f"sum(S1) = {results_interp['S1'].sum():.4f}")

# Analyze rankings
metrics = ["MI", "dCor", "PermMean", "S1", "ST"]
print("\nRankings by different metrics:")
for metric in metrics:
    ranking = results_interp.sort_values(metric, ascending=False).index.tolist()
    print(f"{metric:<10}: {ranking}")

exact_s1 = dict(zip(param_names_interp, [v_x1 / V, v_x2 / V, 0.0, v_x4 / V, v_x5 / V]))
print("\nS1 against the analytic values:")
print(f"{'Param':<18}{'S1 (GP)':>10}{'S1 (exact)':>12}{'abs err':>10}")
print("-" * 50)
for name in param_names_interp:
    got = results_interp.loc[name, "S1"]
    print(f"{name:<18}{got:>10.4f}{exact_s1[name]:>12.4f}{abs(got - exact_s1[name]):>10.4f}")

print("\nInterpretation insights:")
print("- x1_linear dominates: 2*x1 carries ~84% of the variance on its own.")
print("- x3_interaction is the instructive one. Its S1 is ~0 (analytically")
print("  exactly 0) while its ST is clearly positive -- it acts ONLY through")
print("  the x1*x3 term. Reading S1 alone would call it irrelevant.")
print("- A small NEGATIVE S1, as x3 shows here, is estimator noise around zero,")
print("  not a negative sensitivity. Report it as 0.")
print("- The screening metrics (MI, dCor, PermMean) rank x3 second because they")
print("  see the dependence without decomposing it; Sobol S1 does decompose it.")
print("  Disagreement between the two families is informative, not a bug.")
print("- x5_noise ranks last on every metric, as it should.")

# Challenge 1: Correlated inputs
print("Challenge 1: Handling Correlated Inputs")
print("=" * 40)

# Generate correlated inputs
n_samples = 1200
mean = [0, 0, 0]
cov = [[1.0, 0.7, 0.3],    # x1 strongly correlated with x2
       [0.7, 1.0, 0.1],    # x2 moderately correlated with x3
       [0.3, 0.1, 1.0]]    # x3 weakly correlated with others

X_corr = np.random.multivariate_normal(mean, cov, n_samples)

# Simple linear model
def correlated_model(X):
    return 2*X[:, 0] + X[:, 1] + 0.5*X[:, 2]

y_corr = correlated_model(X_corr)

# Run GSA
results_corr, extras_corr = gsa_for_target(
    X_corr, y_corr,
    param_names=["x1", "x2", "x3"],
    scaler='standard',  # Better for correlated normal variables
    enable_perm=True,
    enable_sobol=True
)

print("\nCorrelated inputs - GSA results:")
print(results_corr[['MI', 'dCor', 'PermMean', 'S1', 'ST']].round(4))

# Check correlation matrix
corr_matrix = np.corrcoef(X_corr.T)
print("\nInput correlation matrix:")
print(pd.DataFrame(corr_matrix, index=["x1", "x2", "x3"], columns=["x1", "x2", "x3"]).round(3))

print("\nNote: With correlated inputs, interpretation becomes more complex.")
print("Permutation importance may be more reliable than Sobol indices.")

# Challenge 2: Computational efficiency for large problems
print("\nChallenge 2: Computational Efficiency")
print("=" * 35)

# Simulate a larger problem
n_params_large = 20
n_samples_large = 800

print(f"Large problem: {n_params_large} parameters, {n_samples_large} samples")

# Generate data
X_large = np.random.uniform(-1, 1, (n_samples_large, n_params_large))

# Simple additive model with decreasing importance
def large_model(X):
    weights = np.exp(-np.arange(X.shape[1]) * 0.3)  # Exponentially decreasing importance
    return np.sum(X * weights, axis=1)

y_large = large_model(X_large)
param_names_large = [f"param_{i+1}" for i in range(n_params_large)]

# Fast GSA configuration (skip expensive methods)
import time
start_time = time.time()

results_large, extras_large = gsa_for_target(
    X_large, y_large,
    param_names=param_names_large,
    scaler='minmax',
    enable_perm=False,    # Skip permutation importance (expensive)
    enable_gp=True,       # Keep GP (relatively fast)
    enable_sobol=False,   # Skip Sobol (expensive for many parameters)
)

elapsed_time = time.time() - start_time

print(f"\nFast GSA completed in {elapsed_time:.2f} seconds")
print("\nTop 10 most important parameters:")
top_params = results_large.sort_values('MI', ascending=False).head(10)
print(top_params[['MI', 'dCor']].round(4))

print(f"surrogate collapsed: {extras_large['surrogate_collapsed']}")

print("\nComputational tips:")
print("- Use MI and dCor for quick screening")
print("- Enable GP for moderate-sized problems (< 50 parameters)")
print("- Use Sobol indices only for final analysis of top parameters")
print("- Consider chunked processing for very large datasets")
print("- GP fitting is O(n^3) in the number of SAMPLES, so halving n is worth")
print("  far more than dropping parameters. n=800 here rather than 5000 for")
print("  exactly that reason -- the default now runs the optimiser 5 times.")