from __future__ import annotations
import matplotlib.pyplot as plt
import numpy as np
from example_models import get_linear_chain_2v
from mxlpy import scan, sensitivity
Global sensitivity analysis¶
Before committing to an expensive analysis — parameter fitting, profile-likelihood identifiability, or a full variance-based Sobol decomposition — it pays to ask a cheaper question first: which parameters actually move the output at all?
The Morris method answers exactly that. It is an elementary-effects screening method: it perturbs one parameter at a time along a handful of random trajectories through parameter space and measures how much the output changes. With only n_trajectories * (k + 1) model evaluations (k = number of parameters) it sorts parameters into three buckets:
| Reading | Interpretation |
|---|---|
low mu_star |
unimportant — safe to fix to a constant |
high mu_star, low sigma |
important and roughly linear |
high mu_star, high sigma |
important and nonlinear or interacting |
where mu_star is the mean absolute elementary effect and sigma its standard deviation across trajectories.
The output function¶
Sensitivity is always measured with respect to some scalar output of the model. Rather than ship a fixed menu of reductions, sensitivity.morris asks you for an output function:
output(model, samples: pd.DataFrame) -> np.ndarray
It receives the model and the full matrix of parameter samples (one row per sample, one column per parameter) and returns one scalar per row. The natural way to write it is to lean on mxlpy.scan, which already evaluates the whole sample matrix in parallel, caches results, and — importantly — returns NaN for any sample whose integration failed. That NaN flows straight through to the sensitivity indices, so failed runs surface honestly instead of being silently masked.
Any reduction scan can express is available for free:
- steady-state value of a species:
scan.steady_state(model, to_scan=samples).variables[name] - peak value over a time course:
scan.time_course(model, to_scan=samples, time_points=...).get_agg_per_run("max")[name] - area under the curve: pass
np.trapezoidas the aggregator toget_agg_per_run
Here we screen the steady-state concentration of x in a simple linear chain $\varnothing \xrightarrow{v_1} x \xrightarrow{v_2} y \xrightarrow{v_3} \varnothing$ against its three rate constants.
def steady_state_x(model, samples):
return scan.steady_state(model, to_scan=samples).variables["x"].to_numpy()
result = sensitivity.morris(
get_linear_chain_2v(),
output=steady_state_x,
param_bounds={"k1": (0.5, 2.0), "k2": (0.5, 2.0), "k_out": (0.5, 2.0)},
n_trajectories=5,
seed=0,
)
result.sort_values("mu_star", ascending=False)
0%| | 0/20 [00:00<?, ?it/s]
5%|▌ | 1/20 [00:05<01:41, 5.32s/it]
20%|██ | 4/20 [00:05<00:16, 1.04s/it]
65%|██████▌ | 13/20 [00:05<00:01, 4.13it/s]
100%|██████████| 20/20 [00:05<00:00, 3.58it/s]
| mu | mu_star | sigma | mu_star_conf | |
|---|---|---|---|---|
| k2 | -2.100000e+00 | 2.100000e+00 | 1.341641e+00 | 1.009701e+00 |
| k1 | 1.900000e+00 | 1.900000e+00 | 1.024695e+00 | 8.569927e-01 |
| k_out | -1.477828e-10 | 1.477953e-10 | 3.304699e-10 | 2.676228e-10 |
The result is a labelled pandas.DataFrame indexed by parameter name with columns mu, mu_star, sigma, and mu_star_conf (a bootstrap confidence interval on mu_star).
The reading is immediate: at steady state $x = k_1 / k_2$, so k1 and k2 dominate while k_out — which only drains the downstream species y — has an elementary effect of essentially zero. A screening run like this lets you drop k_out from a subsequent fit without losing anything.
Visualising the screen¶
Two plots make the screen easy to read: a ranking of mu_star (with its confidence interval), and a mu_star-vs-sigma scatter that separates linear from nonlinear influence.
ranked = result.sort_values("mu_star")
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(9, 3.5))
ax1.barh(ranked.index, ranked["mu_star"], xerr=ranked["mu_star_conf"])
ax1.set(xlabel=r"$\mu^*$ (mean absolute effect)", title="Parameter ranking")
ax2.scatter(result["mu_star"], result["sigma"])
for name, row in result.iterrows():
ax2.annotate(name, (row["mu_star"], row["sigma"]))
ax2.set(xlabel=r"$\mu^*$", ylabel=r"$\sigma$", title="Linear vs. nonlinear")
fig.tight_layout()
plt.show()
Cost and reproducibility¶
Morris is cheap: with n_trajectories=20 and k=3 parameters the screen above used 20 * (3 + 1) = 80 model evaluations — a tiny fraction of what a variance-based method would need. That is the point: screen first with Morris, then quantify the survivors with a more expensive method.
Pass a seed for reproducible sampling and bootstrap confidence intervals, and tune num_levels (default 4) to control the granularity of the sampling grid.
n_params = 3
n_trajectories = 20
print(f"Evaluations: {n_trajectories * (n_params + 1)}")
Evaluations: 80
From screening to quantification: Sobol indices¶
Morris tells you which parameters matter; it does not tell you how much each one explains, nor cleanly separate a parameter's own effect from the effect of its interactions. The Sobol method does both. It is a variance-based method: it decomposes the total variance of the output into the share attributable to each parameter.
It returns two indices per parameter:
| Index | Meaning |
|---|---|
S1 (first-order) |
fraction of output variance explained by that parameter alone |
ST (total-order) |
fraction explained by that parameter including all of its interactions |
The gap ST - S1 is the interaction signal: if it is large, the parameter matters mostly through other parameters rather than on its own.
sobol takes the same output and param_bounds contract as morris, so the reduction you already wrote carries over unchanged. The one new knob is n_samples (the Saltelli base sample size), which must be a power of two — passing anything else raises a ValueError. The total cost is n_samples * (k + 2) model evaluations, which is why you screen with Morris first and reserve Sobol for the survivors.
result_sobol = sensitivity.sobol(
get_linear_chain_2v(),
output=steady_state_x,
param_bounds={"k1": (0.5, 2.0), "k2": (0.5, 2.0), "k_out": (0.5, 2.0)},
n_samples=512,
seed=0,
)
result_sobol.sort_values("ST", ascending=False)
0%| | 0/2560 [00:00<?, ?it/s]
0%| | 1/2560 [00:05<3:48:40, 5.36s/it]
0%| | 2/2560 [00:05<1:39:30, 2.33s/it]
0%| | 4/2560 [00:05<38:48, 1.10it/s]
3%|▎ | 66/2560 [00:05<01:22, 30.07it/s]
4%|▍ | 109/2560 [00:05<00:44, 54.91it/s]
6%|▋ | 163/2560 [00:06<00:25, 94.32it/s]
8%|▊ | 215/2560 [00:06<00:16, 138.47it/s]
10%|█ | 265/2560 [00:06<00:12, 184.00it/s]
13%|█▎ | 323/2560 [00:06<00:09, 244.77it/s]
15%|█▍ | 373/2560 [00:06<00:07, 289.74it/s]
17%|█▋ | 427/2560 [00:06<00:06, 339.83it/s]
19%|█▊ | 477/2560 [00:06<00:05, 375.15it/s]
21%|██ | 529/2560 [00:06<00:04, 407.43it/s]
23%|██▎ | 583/2560 [00:06<00:04, 440.34it/s]
25%|██▍ | 635/2560 [00:06<00:04, 453.64it/s]
27%|██▋ | 691/2560 [00:07<00:03, 481.24it/s]
29%|██▉ | 744/2560 [00:07<00:03, 494.89it/s]
31%|███ | 797/2560 [00:07<00:03, 497.51it/s]
33%|███▎ | 849/2560 [00:07<00:03, 501.10it/s]
35%|███▌ | 904/2560 [00:07<00:03, 514.03it/s]
37%|███▋ | 957/2560 [00:07<00:03, 516.51it/s]
39%|███▉ | 1010/2560 [00:07<00:02, 518.41it/s]
42%|████▏ | 1063/2560 [00:07<00:02, 520.71it/s]
44%|████▎ | 1119/2560 [00:07<00:02, 531.91it/s]
46%|████▌ | 1173/2560 [00:07<00:02, 529.56it/s]
48%|████▊ | 1227/2560 [00:08<00:02, 518.99it/s]
50%|█████ | 1282/2560 [00:08<00:02, 527.57it/s]
52%|█████▏ | 1336/2560 [00:08<00:02, 529.13it/s]
54%|█████▍ | 1391/2560 [00:08<00:02, 530.89it/s]
56%|█████▋ | 1445/2560 [00:08<00:02, 528.89it/s]
59%|█████▊ | 1498/2560 [00:08<00:02, 524.55it/s]
61%|██████ | 1551/2560 [00:08<00:01, 517.44it/s]
63%|██████▎ | 1603/2560 [00:08<00:01, 517.75it/s]
65%|██████▍ | 1656/2560 [00:08<00:01, 521.21it/s]
67%|██████▋ | 1709/2560 [00:08<00:01, 521.16it/s]
69%|██████▉ | 1762/2560 [00:09<00:01, 515.92it/s]
71%|███████ | 1816/2560 [00:09<00:01, 516.40it/s]
73%|███████▎ | 1870/2560 [00:09<00:01, 517.25it/s]
75%|███████▌ | 1923/2560 [00:09<00:01, 520.00it/s]
77%|███████▋ | 1976/2560 [00:09<00:01, 515.77it/s]
79%|███████▉ | 2030/2560 [00:09<00:01, 517.86it/s]
81%|████████▏ | 2084/2560 [00:09<00:00, 520.36it/s]
84%|████████▎ | 2138/2560 [00:09<00:00, 520.48it/s]
86%|████████▌ | 2191/2560 [00:09<00:00, 521.78it/s]
88%|████████▊ | 2244/2560 [00:10<00:00, 516.52it/s]
90%|████████▉ | 2296/2560 [00:10<00:00, 496.23it/s]
92%|█████████▏| 2352/2560 [00:10<00:00, 512.57it/s]
94%|█████████▍| 2404/2560 [00:10<00:00, 514.09it/s]
96%|█████████▌| 2456/2560 [00:10<00:00, 511.92it/s]
98%|█████████▊| 2508/2560 [00:10<00:00, 509.46it/s]
100%|██████████| 2560/2560 [00:10<00:00, 240.55it/s]
| S1 | ST | S1_conf | ST_conf | |
|---|---|---|---|---|
| k2 | 5.427225e-01 | 6.108146e-01 | 9.621455e-02 | 8.662555e-02 |
| k1 | 3.829592e-01 | 4.491552e-01 | 6.922274e-02 | 6.469512e-02 |
| k_out | -5.538662e-11 | 1.247718e-19 | 6.100698e-11 | 8.766044e-20 |
The result is a labelled pandas.DataFrame indexed by parameter name with columns S1, ST, and their bootstrap confidence intervals S1_conf, ST_conf.
The story matches the Morris screen but is now quantitative: with $x = k_1 / k_2$ at steady state, k1 and k2 carry essentially all of the variance while k_out contributes nothing. Here S1 and ST are close for every parameter, telling us the influence is largely additive over this box rather than interaction-driven.
ranked = result_sobol.sort_values("ST")
fig, ax = plt.subplots(figsize=(6, 3.5))
y = np.arange(len(ranked))
ax.barh(
y - 0.2, ranked["S1"], height=0.4, xerr=ranked["S1_conf"], label="S1 (first-order)"
)
ax.barh(
y + 0.2, ranked["ST"], height=0.4, xerr=ranked["ST_conf"], label="ST (total-order)"
)
ax.set_yticks(y, ranked.index)
ax.set(xlabel="Sobol index", title="Variance decomposition")
ax.legend()
fig.tight_layout()
plt.show()
Screen, then quantify. Morris is the cheap first pass that discards parameters you can safely fix; Sobol is the more expensive second pass that puts numbers on the ones that survive. Used together they turn a high-dimensional parameter space into a short, ranked list of the parameters worth fitting or measuring.