Chapter 31: Ensemble Scenarios and Sensitivity Analysis¶
Part VII — Scaling, Migration & Reference
Draft chapter
This chapter is an early draft of the ensemble-and-sensitivity-analysis workflow, not a settled final version — expect it to evolve. The ensemble example itself, built on Chapter 23's PredatorPrey model, was tested end to end, including the Sobol sampling and analysis pipeline; what's likely to change is scope and framing as this topic develops alongside the rest of the book.
Learning Objectives¶
By the end of this chapter you will be able to:
- Explain why a single deterministic model run is rarely enough to trust a result
- Design and run a parameter ensemble
- Compute global sensitivity indices with
SALib, and read what they mean - Interpret an ensemble's spread as a statement about uncertainty, not a single "correct" answer
# Standard imports
import contextlib
import io
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
Every model this book has run so far — Chapter 23's PredatorPrey, Chapter 24's FireModel, Chapter 26's eventual LUCC executors — has run exactly once, with one fixed set of parameters, producing one deterministic answer. That answer is only as trustworthy as the parameter values feeding it, and in a real project those values are rarely known exactly — a growth rate estimated from noisy field data, a mortality rate borrowed from a different study system entirely. This chapter is about turning "here's what the model says with these parameters" into "here's what the model says across the plausible range of parameters, and here's which of those parameters actually matters."
Introduction to Ensemble Modeling¶
An ensemble is simply many runs of the same model, each with a different parameter combination, treated as one collective result rather than many separate ones. The reason to build one at all is that a single run answers the wrong question — "what does the model predict" is less useful, on its own, than "how much does the model's prediction change as the inputs vary within their plausible range." A tight ensemble (every run lands close together) says the model's conclusion is robust to what you don't know precisely about the parameters; a wide, divergent ensemble says the opposite — and that distinction is invisible from any single run, however carefully chosen.
Designing an Ensemble: Parameter Sampling¶
The naive approach — loop over every combination of a handful of values per parameter — scales catastrophically: 3 parameters at 10 values each is already 1,000 runs, and every parameter added multiplies the total again. Quasi-random sampling (Sobol sequences, used below) instead spreads a fixed sampling budget evenly across the whole parameter space, covering it far more uniformly than either a full grid or genuinely random points would for the same number of runs — the property Sobol analysis specifically needs to work correctly, not an incidental convenience.
SALib's sobol.sample() needs a problem definition — parameter names and their bounds — and returns exactly the array of parameter combinations a Sobol sensitivity analysis requires:
from SALib.sample import sobol as sobol_sample
problem = {
"num_vars": 3,
"names": ["prey_growth", "prey_death_pred", "pred_death"],
"bounds": [
[0.05, 0.15], # prey_growth
[0.0005, 0.002], # prey_death_pred
[0.01, 0.05], # pred_death
],
}
param_values = sobol_sample.sample(problem, N=32)
print("Parameter combinations generated:", param_values.shape[0])
print("First 3 rows:")
print(param_values[:3])
N=32 here is a base sample size, not the run count directly — Sobol's specific sampling scheme needs extra evaluation points to separate each parameter's individual effect from its interactions with the others, so the actual number of model runs comes out several times larger than N alone. That overhead is real, and it's exactly what the Global Sensitivity Analysis section below is paying for.
Running an Ensemble¶
Reusing Chapter 23's PredatorPrey unchanged, one run per sampled parameter combination — only the driving loop is new, not the model:
from dissmodel.core import Environment
from dissmodel_sysdyn.models import PredatorPrey
def run_predator_prey(prey_growth, prey_death_pred, pred_death):
with contextlib.redirect_stdout(io.StringIO()): # silence Environment's per-run log
env = Environment(start_time=0, end_time=100)
pp = PredatorPrey(
predator=40.0, prey=1000.0,
prey_growth=prey_growth, prey_death_pred=prey_death_pred,
pred_death=pred_death, pred_growth_kills=0.00002,
)
env.run()
return pp.prey # final prey population — this ensemble's output of interest
final_prey = np.array([run_predator_prey(*row) for row in param_values])
print(f"Ran {len(final_prey)} simulations")
print(f"Final prey population: min={final_prey.min():.1f}, median={np.median(final_prey):.1f}, max={final_prey.max():.1f}")
Interpreting Ensemble Output¶
A histogram of the ensemble's output is the first thing worth looking at — not the mean alone, which can hide a distribution that's actually bimodal or wildly skewed:
fig, ax = plt.subplots(figsize=(8, 4))
ax.hist(final_prey, bins=30, color="seagreen", alpha=0.7)
ax.axvline(np.median(final_prey), color="black", linestyle="--", label="median")
ax.set_xlabel("Final prey population (tick 100)")
ax.set_ylabel("Number of ensemble runs")
ax.set_title(f"Ensemble spread — {len(final_prey)} runs")
ax.legend()
plt.show()
A wide, heavy-tailed spread here means the model's tick-100 prey population depends strongly on exactly which parameter values within the sampled ranges turn out to be correct — the ecosystem sometimes essentially collapses (values near zero) and sometimes overshoots into the thousands, and no single run captures that range. Reporting only "the model predicts X prey" without also reporting this spread would misrepresent how much the model actually knows.
Global Sensitivity Analysis: Sobol Indices¶
A wide ensemble spread raises the obvious next question: spread caused by which parameter? SALib's Sobol analysis decomposes the ensemble's total output variance into a share attributable to each parameter individually — a first-order index, S1 — using nothing but the same param_values and final_prey arrays already in hand:
from SALib.analyze import sobol as sobol_analyze
Si = sobol_analyze.analyze(problem, final_prey, print_to_console=False)
sensitivity = pd.DataFrame({
"parameter": problem["names"],
"S1": Si["S1"],
"S1_conf": Si["S1_conf"],
})
sensitivity.sort_values("S1", ascending=False)
fig, ax = plt.subplots(figsize=(6, 4))
ax.barh(sensitivity["parameter"], sensitivity["S1"], xerr=sensitivity["S1_conf"], color="steelblue")
ax.set_xlabel("First-order Sobol index (S1)")
ax.set_title("Which parameter drives final prey population?")
plt.show()
An S1 of 0 would mean that parameter, varied alone, explains none of the ensemble's output variance; an S1 near 1 would mean it explains nearly all of it. S1 values across all parameters summing to noticeably less than 1 signals that interactions between parameters — not captured by any single S1 — also matter, which is precisely the kind of effect a naive one-parameter-at-a-time sweep would have missed entirely, and exactly why the Sobol sampling scheme in Designing an Ensemble exists in the first place.
Interpreting GSA for Policy¶
The parameter with the highest S1 is where uncertainty reduction pays off the most: narrowing that one parameter's plausible range — better field data, a more targeted study — shrinks the ensemble's spread more than narrowing any other single parameter would, for the same amount of research effort spent. The parameters with the lowest S1 are, symmetrically, the ones where a modeler can reasonably use a rough estimate without materially changing the model's conclusions — refining them further is effort spent where the model already doesn't care much.
This is the sentence a GSA result is actually for: not "parameter X is important" as an abstract fact, but "further precision on parameter X's estimate would change this model's usefulness more than precision on any other single parameter would" — a direct, actionable answer to "where should the next unit of research effort go," grounded in the model itself rather than intuition.
Exercises¶
- A different output. Rerun the ensemble from Running an Ensemble, but track
pp.predator(final predator population) instead ofpp.prey. Rerun the Sobol analysis on this new output — does the same parameter still dominateS1, or does a different one take the lead? - Increase the sample size. Rerun Designing an Ensemble with
N=64instead of32. How many total simulations does that produce, and do theS1values change meaningfully from theN=32result — or have they already converged? - A fourth parameter. Add
pred_growth_killstoproblem'snames/bounds(try[0.00001, 0.00004]), rerun the whole pipeline, and report where it ranks among the fourS1values. - From index to recommendation. Using your own ensemble's
S1results, write two sentences in the style of Interpreting GSA for Policy — one about where refining parameter estimates would help most, one about where it wouldn't.
# Your code here
Summary¶
Key concepts introduced¶
- Why a single deterministic run under-reports what a model actually knows: it hides how much the answer would change under different, equally plausible parameter values
- Quasi-random (Sobol) sampling, spreading a fixed run budget evenly across parameter space — a requirement of Sobol analysis, not an incidental choice
- Building an ensemble by reusing an existing model (Chapter 23's
PredatorPrey) completely unchanged, varying only the driving loop around it - Reading an ensemble's output distribution — spread, not just a mean — as the honest summary of a result
- First-order Sobol indices (
S1), decomposing output variance by parameter, and what a low total across allS1values implies about parameter interactions - Turning a sensitivity ranking into an actionable research-effort recommendation, not just a ranked list of numbers
This closes the paradigm-and-analysis arc of this book. Chapter 32 turns outward — a concept-by-concept guide for migrating an existing TerraME or LUCCME model into the ecosystem this book has spent fifteen chapters building.
Further Reading¶
- Saltelli, A. et al. (2008). Global Sensitivity Analysis: The Primer. Wiley — the standard reference behind Sobol indices and the sampling scheme this chapter uses
- SALib documentation: https://salib.readthedocs.io/
- Sobol, I. M. (2001). "Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates." Mathematics and Computers in Simulation, 55(1-3), 271-280
- Saltelli, A. et al. (2010). "Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index." Computer Physics Communications, 181(2), 259-270