"""Sobol sensitivity analysis for model parameters (Priority 7.5)."""
from typing import Callable, Dict, List, Tuple
import numpy as np
from SALib.analyze import sobol as sobol_analyze
from SALib.sample import sobol
[docs]
class SobolSensitivity:
"""
Perform Sobol sensitivity analysis on model parameters.
Determines which parameters most influence a given output metric.
"""
[docs]
def __init__(
self,
problem: Dict,
n_samples: int = 1024,
calc_second_order: bool = True,
n_processors: int = 1,
):
"""
Parameters
----------
problem : dict
SALib problem definition with 'num_vars', 'names', 'bounds'.
n_samples : int
Number of samples (must be power of 2).
calc_second_order : bool
Whether to compute second-order interactions.
n_processors : int
Number of CPUs for parallel evaluation.
"""
self.problem = problem
self.n_samples = n_samples
self.calc_second_order = calc_second_order
self.n_processors = n_processors
# Generate parameter samples
self.param_values = sobol.sample(problem, n_samples, calc_second_order=calc_second_order)
[docs]
def evaluate_model(
self, model_func: Callable[[np.ndarray], float], batch_size: int = 100
) -> np.ndarray:
"""
Evaluate model for all parameter samples.
Parameters
----------
model_func : function
Takes a parameter vector (1D array) and returns a scalar output.
batch_size : int
Number of samples to evaluate in each batch (for memory).
Returns
-------
Y : np.ndarray
Model outputs for each sample.
"""
n = len(self.param_values)
Y = np.zeros(n)
for i in range(0, n, batch_size):
batch = self.param_values[i : i + batch_size]
batch_results = []
for params in batch:
try:
result = model_func(params)
batch_results.append(result)
except Exception as e:
print(f"Error with params {params}: {e}")
batch_results.append(np.nan)
Y[i : i + batch_size] = batch_results
return Y
[docs]
def analyze(self, Y: np.ndarray) -> Dict:
"""
Compute Sobol indices.
Returns
-------
dict
Contains 'S1' (first-order), 'ST' (total-order), and optionally 'S2' (second-order).
"""
Si = sobol_analyze.analyze(
self.problem, Y, calc_second_order=self.calc_second_order, print_to_console=False
)
return Si
[docs]
def get_top_parameters(
self, Si: Dict, threshold: float = 0.05, order: str = "total"
) -> List[Tuple[str, float]]:
"""
Return parameters with sensitivity above threshold, sorted by importance.
Parameters
----------
Si : dict
Result from analyze().
threshold : float
Minimum sensitivity value to include.
order : str
'first' for S1, 'total' for ST.
"""
if order == "first":
values = Si["S1"]
else:
values = Si["ST"]
names = self.problem["names"]
important = [
(names[i], values[i])
for i in range(len(names))
if not np.isnan(values[i]) and values[i] > threshold
]
important.sort(key=lambda x: x[1], reverse=True)
return important
[docs]
def create_default_problem() -> Dict:
"""
Create a default problem definition for MERCURIAL parameters.
Returns
-------
dict
Problem specification with parameter names and bounds.
"""
return {
"num_vars": 12,
"names": [
"kappa_bias", # Top-down biasing strength
"gamma0", # Base decay rate
"S_crit", # Critical entropy threshold
"k_align", # Alignment Boltzmann constant
"tau_align", # Alignment decay time
"alpha_form", # Impression formation rate
"beta_steepness", # Modality threshold steepness
"coupling_base", # Cross-level baseline coupling
"decay_length", # LADDER adjacency decay length
"noise_amplitude", # Stochastic noise strength
"prior_strength", # Bayesian prior strength
"learning_rate", # Perception learning rate
],
"bounds": [
[0.01, 0.5], # kappa_bias
[0.001, 0.1], # gamma0
[5.0, 20.0], # S_crit
[0.5, 5.0], # k_align
[0.1, 10.0], # tau_align
[0.01, 0.5], # alpha_form
[1.0, 10.0], # beta_steepness
[0.01, 0.5], # coupling_base
[0.5, 2.0], # decay_length
[0.001, 0.1], # noise_amplitude
[0.1, 2.0], # prior_strength
[0.01, 0.5], # learning_rate
],
}
[docs]
def example_model(params: np.ndarray) -> float:
"""
Deterministic synthetic model for sensitivity testing (no random noise).
Output is a smooth function of parameters.
"""
(
kappa_bias,
gamma0,
S_crit,
k_align,
tau_align,
alpha_form,
beta_steepness,
coupling_base,
decay_length,
noise_amplitude,
prior_strength,
learning_rate,
) = params
# Deterministic output – no random noise
output = (
0.5 * kappa_bias
+ 0.3 * np.exp(-gamma0)
+ 0.2 * (S_crit / 10.0)
+ 0.1 * k_align
+ 0.05 * tau_align
+ 0.1 * alpha_form
+ 0.05 * beta_steepness
+ 0.1 * coupling_base
+ 0.05 * decay_length
+ 0.02 * noise_amplitude
+ 0.1 * prior_strength
+ 0.03 * learning_rate
)
return output