Comparison, selection and sensitivity¶

comparison ¶

Scenario comparison against a baseline.

Comparison ¶

Comparison(rows: Sequence[ComparisonRow], confidence: float, adjust: AdjustMethod = 'none')

Differences scenario - baseline for each metric, with confidence intervals.

The intervals quantify simulation sampling error only. A significant difference means the model, as specified, responds to the change; it says nothing about whether the real system would (that depends on model validity). With adjust="none" no multiple-comparison correction is applied; see :func:compare for the alternatives.

compare ¶

compare(result: ExperimentResult, baseline: str = 'baseline', scenarios: Iterable[str] | None = None, metrics: Iterable[str] | None = None, *, confidence: float = 0.95, adjust: AdjustMethod = 'none') -> Comparison

Compare scenarios of one experiment with its baseline scenario.

When the experiment used common random numbers and both scenarios have the same replication seeds, a paired-t interval is used (replication i of each scenario saw the same random streams); otherwise Welch's t.

adjust corrects for testing many (metric, scenario) pairs at once: "bonferroni" widens every interval to level 1 - alpha/m (simultaneous coverage) and adjusts p-values; "holm" (family-wise error) and "bh" (false discovery rate) adjust p-values only, and significant then means p_adjusted < 1 - confidence.

compare_samples ¶

compare_samples(baseline: Mapping[str, Sequence[float]], scenario: Mapping[str, Sequence[float]], *, paired: bool = False, confidence: float = 0.95, adjust: AdjustMethod = 'none', baseline_name: str = 'baseline', scenario_name: str = 'scenario') -> Comparison

Compare two {metric: values} samples (e.g. collected outside an Experiment).

selection ¶

Ranking and selection: pick the best scenario with a statistical guarantee.

:func:select_best implements the fully sequential procedure KN of Kim and Nelson (2001). It runs every scenario n0 times, then adds one replication at a time to the scenarios still in contention and eliminates any that is clearly worse than another. With probability at least confidence the scenario selected is the true best, or within indifference of it (the smallest difference you care about).

Replications use common random numbers, which KN exploits: it works with the variance of paired differences, which CRN makes small. Compared with a fixed number of replications per scenario, clearly inferior scenarios stop early and the budget goes to the close contenders.

select_best ¶

select_best(model: Model, scenarios: Sequence[Scenario] | Mapping[str, Mapping[str, Any]], metric: str, *, indifference: float, minimize: bool = True, confidence: float = 0.95, n0: int = 10, max_replications: int = 500, seed: int = 0) -> SelectionResult

Select the best of scenarios on metric with procedure KN.

indifference (delta) is in the metric's units: differences smaller than this are treated as ties, and any scenario within delta of the best counts as a correct selection. Smaller delta means more replications. If max_replications per scenario is reached first, the scenario with the best mean among those remaining is returned with converged=False.

sensitivity ¶

Sensitivity analysis: which inputs move which outputs, and by how much.

Methods implemented (all simple and well understood):

:func:one_at_a_time - change one parameter at a time around a base point, with common random numbers, and report the change in each output plus a paired confidence interval and an elasticity.
:func:finite_difference - local derivative d output / d parameter by central (or forward) differences with common random numbers.
:func:correlation_sensitivity - global screening from Monte Carlo samples: Pearson, Spearman rank correlation and standardised regression coefficients (SRC, with the regression R^2 so you can tell whether a linear summary is adequate).
:func:sobol_indices - variance-based global sensitivity: first-order and total-effect Sobol indices (Saltelli 2010 / Jansen estimators) with bootstrap confidence intervals. Captures interactions and non-linear effects at a cost of N * (d + 2) model evaluations.
:func:morris_screening - Morris elementary effects: a cheap global screen (r * (d + 1) evaluations) that separates negligible inputs from important ones (mu_star) and flags non-linearity or interactions (sigma). Use it to prune many inputs before a Sobol analysis.

Correlation measures are cheap screening tools that miss interactions and non-monotonic effects; use Sobol indices when that matters.

SensitivityResult ¶

SensitivityResult(rows: Sequence[SensitivityRow], info: Mapping[str, Any] | None = None)

ranking ¶

ranking(output: str | None = None, method: str | None = None) -> list[SensitivityRow]

Rows sorted by absolute effect size, largest first.

one_at_a_time ¶

one_at_a_time(model: Model, changes: Mapping[str, Sequence[Any]], outputs: Sequence[str], *, base: Mapping[str, Any] | None = None, replications: int = 10, seed: int = 0, confidence: float = 0.95) -> SensitivityResult

Vary each parameter in changes alone; report output change vs the base point.

value is the mean change in the output; detail includes the elasticity (dy/y) / (dx/x) when both base values are non-zero numbers.

finite_difference ¶

finite_difference(model: Model, parameters: Sequence[str], outputs: Sequence[str], *, base: Mapping[str, Any] | None = None, relative_step: float = 0.05, replications: int = 10, seed: int = 0, scheme: Literal['central', 'forward'] = 'central', confidence: float = 0.95) -> SensitivityResult

Local derivatives by finite differences with common random numbers.

Integer parameters use a step of at least 1. The CI comes from the per-replication paired differences, so it reflects simulation noise only; a step that is too small amplifies noise, too large adds bias.

correlation_sensitivity ¶

correlation_sensitivity(mc: MonteCarloResult, outputs: Sequence[str] | None = None, *, methods: Sequence[Literal['pearson', 'spearman', 'src']] = ('pearson', 'spearman', 'src')) -> SensitivityResult

Global screening from Monte Carlo samples of uncertain inputs.

sobol_indices ¶

sobol_indices(model: Callable[..., Any] | Model, parameters: Mapping[str, DistributionLike | Mapping[str, Any]], n: int = 1024, *, outputs: Sequence[str] | None = None, seed: int = 0, fixed: Mapping[str, Any] | None = None, vectorized: bool = False, confidence: float = 0.95, n_bootstrap: int = 500) -> SensitivityResult

First-order and total-effect Sobol indices with bootstrap confidence intervals.

Uses the Saltelli (2010) design: two independent sample matrices A and B of n rows plus, for each of the d parameters, A with that column taken from B. First-order indices use the Saltelli (2010) estimator, total effects the Jansen (1999) estimator. Inputs are assumed independent. For simulation models every matrix row j uses the same replication seed (common random numbers), which keeps simulation noise from swamping the index estimates; noise still adds some bias, so use enough replication length or average several runs per row if outputs are very noisy.

value is the index estimate (method = sobol-first or sobol-total). Estimates can fall slightly outside [0, 1] with small n; the bootstrap interval shows how precise they are.

morris_screening ¶

morris_screening(model: Callable[..., Any] | Model, ranges: Mapping[str, tuple[float, float]], r: int = 20, *, levels: int = 4, outputs: Sequence[str] | None = None, seed: int = 0, fixed: Mapping[str, Any] | None = None, confidence: float = 0.95, n_bootstrap: int = 500) -> SensitivityResult

Morris elementary-effects screening over ranges ({name: (low, high)}).

Builds r random one-at-a-time trajectories on a levels-level grid (Morris 1991), each costing d + 1 evaluations. Effects are measured on inputs scaled to [0, 1], so they are comparable across parameters and in output units. Per parameter and output it reports:

morris-mu_star - mean absolute effect (Campolongo et al. 2007), the importance measure, with a bootstrap CI over trajectories;
morris-mu - mean signed effect (direction; can cancel out);
morris-sigma - standard deviation of effects: large relative to mu_star means non-linear effects or interactions.

Integer model parameters are rounded. For simulation models every point of a trajectory uses the same replication seed (common random numbers).

warmup ¶

Warm-up (initial transient) detection.

Steady-state estimates are biased by the start-up period of a run (an empty system has short queues). :func:mser implements MSER (White 1997; MSER-5 with batches of 5, Franklin & White 2008): delete the first d observations where d minimises the marginal standard error of the remaining mean. :func:suggest_warmup applies it to a model: it averages a recorded time series across replications (reducing noise, as in Welch's method), and reports the truncation point as a simulation time.

MSER is a heuristic. If the suggested truncation is close to the maximum allowed fraction, the run is probably too short to reach steady state - the result flags this rather than returning a misleading number.

MserResult `dataclass` ¶

MserResult(truncation: int, statistic: float, n: int, batch_size: int, reliable: bool, note: str = '')

truncation `instance-attribute` ¶

truncation: int

Number of original observations to delete.

mser ¶

mser(values: Sequence[float] | FloatArray, batch_size: int = 5, max_fraction: float = 0.5) -> MserResult

MSER-batch_size truncation point for one output sequence (in observation order).

suggest_warmup ¶

suggest_warmup(model: Model, params: Mapping[str, Any] | None = None, *, series: str | None = None, replications: int = 5, seed: int = 0, duration: float | None = None, bins: int = 200, batch_size: int = 5) -> WarmupAdvice

Suggest a warm-up for model from a recorded time series.

series is a key of SimulationResult.series (default: the first resource queue_length). Each replication is run without warm-up, the series is time-averaged into bins equal bins, averaged across replications, and MSER-batch_size picks the truncation bin.

input_uncertainty ¶

How much of the uncertainty comes from the input data?

Simulation confidence intervals usually cover simulation noise only. But input distributions are fitted to finite data, and a different sample of data would give different parameters. :func:input_uncertainty estimates that extra uncertainty by bootstrapping: resample each data set, refit its distribution, run the model, repeat. The spread of the results across bootstrap fits, after removing simulation noise, is the input uncertainty (Barton and Schruben 2001; Barton 2012).

interarrival times observed at a desk, used to set the arrival rate:¶

input_uncertainty(mmc, {"arrival_rate": InputData(gaps, "exponential", lambda d: 1 / d.mean)},¶

"resource.server.wait.mean")¶

InputData `dataclass` ¶

InputData(values: Sequence[float], family: str = 'exponential', transform: Callable[[Distribution[Any]], Any] | None = None)

Observed data for one model parameter.

family is the distribution family to fit ("exponential", "gamma", ...). transform turns the fitted distribution into the parameter value (for example lambda d: 1 / d.mean for a rate); by default the distribution itself is passed, for distribution parameters.

InputUncertaintyResult `dataclass` ¶

InputUncertaintyResult(metric: str, point_estimate: float, simulation_half_width: float, ci_low: float, ci_high: float, input_variance: float, simulation_variance: float, bootstrap_means: list[float], replications: int, confidence: float)

input_share `property` ¶

input_share: float

Share of the estimator's variance caused by finite input data.

input_uncertainty ¶

input_uncertainty(model: Model, inputs: Mapping[str, InputData], metric: str, *, n_bootstrap: int = 30, replications: int = 5, seed: int = 0, fixed: Mapping[str, Any] | None = None, confidence: float = 0.95) -> InputUncertaintyResult

Bootstrap the input data to see how much it moves metric.

For each of n_bootstrap rounds every data set in inputs is resampled with replacement and refitted; the model then runs replications times with common seeds. The total variance of the round means is split into simulation noise (mean within-round variance / replications) and input variance (the rest). The returned interval is the bootstrap percentile interval, which includes both.

rare_events ¶

Rare events: probabilities of extreme outcomes, with honest intervals.

"How likely is it that the hospital overflows?" is a question about a small probability. :func:rare_event_probability runs replications, counts how often metric exceeds threshold and reports an exact (Clopper-Pearson) interval plus the number of replications needed for a target relative error. When the event is too rare to observe often, it also fits a generalised Pareto distribution to the upper tail (peaks over threshold, extreme value theory) and extrapolates the probability, with a bootstrap interval. Extrapolation assumes the tail keeps its shape beyond the data, so treat it as an informed estimate rather than a measurement.

gpd_tail_probability ¶

gpd_tail_probability(values: Sequence[float] | ndarray[Any, Any], threshold: float, *, tail_quantile: float = 0.9, n_bootstrap: int = 200, confidence: float = 0.95, seed: int = 0) -> TailEstimate

P(X > threshold) from a generalised Pareto fit to the values above a high quantile.

clopper_pearson ¶

clopper_pearson(hits: int, n: int, confidence: float = 0.95) -> tuple[float, float]

Exact binomial confidence interval.

rare_event_probability ¶

rare_event_probability(model: Model, metric: str, threshold: float, *, replications: int = 200, seed: int = 0, params: Mapping[str, Any] | None = None, tail: bool = True, confidence: float = 0.95) -> RareEventResult

Estimate P(metric > threshold) from replications independent runs.

Comparison, selection and sensitivity¶

comparison ¶

Comparison ¶

compare ¶

compare_samples ¶

selection ¶

select_best ¶

sensitivity ¶

SensitivityResult ¶

ranking ¶

one_at_a_time ¶

finite_difference ¶

correlation_sensitivity ¶

sobol_indices ¶

morris_screening ¶

warmup ¶

MserResult dataclass ¶

truncation instance-attribute ¶

mser ¶

suggest_warmup ¶

input_uncertainty ¶

interarrival times observed at a desk, used to set the arrival rate:¶

input_uncertainty(mmc, {"arrival_rate": InputData(gaps, "exponential", lambda d: 1 / d.mean)},¶

"resource.server.wait.mean")¶

InputData dataclass ¶

InputUncertaintyResult dataclass ¶

input_share property ¶

input_uncertainty ¶

rare_events ¶

gpd_tail_probability ¶

clopper_pearson ¶

rare_event_probability ¶

MserResult `dataclass` ¶

truncation `instance-attribute` ¶

InputData `dataclass` ¶

InputUncertaintyResult `dataclass` ¶

input_share `property` ¶