Research workflow

SimulSI is designed so that a simulation study can be repeated by someone else - or by you in six months - and give the same numbers.

A reproducible study, step by step

1. Version the model. Put the model in a git repository and set Model(version="..."). Bump the version whenever logic changes.
2. Validate before you trust. Run validate_model (or simulsi validate) for each scenario. Then check the model against something known: a closed-form case (as the M/M/c tests do with Erlang C), historical data, or a hand-worked trace (trace=True).
3. Choose warm-up and run length. For steady-state questions, run suggest_warmup(model, params) (or simulsi warmup model.py): it averages a time series over replications and applies MSER-5. Check the averaged series too; if the result is flagged unreliable, the run is too short to reach steady state. For terminating systems (a working day, a week of operations) the run length is the real horizon and no warm-up is needed.
4. Pilot, then size the experiment. Run ~10 replications and ask result.replication_advice(metric, relative_precision=0.05) how many you need. Re-check after running them.
5. Compare with common random numbers. Keep common_random_numbers=True so scenario differences are estimated with paired intervals.
6. Record everything. result.save(dir) writes the records together with seeds, resolved parameters, model version, git commit (and whether the tree was dirty), timestamp, runtime and environment.

from simulsi import Experiment, ExperimentResult, Scenario
from simulsi.models import mmc

model = mmc.with_options(duration=2_000, warmup=200)
pilot = Experiment(model, Scenario("baseline", {"arrival_rate": 0.9}), replications=10, seed=2026).run()
advice = pilot.replication_advice("resource.server.wait.mean", relative_precision=0.10)
print(advice.relative_half_width, advice.required_n)

study = Experiment(model, Scenario("baseline", {"arrival_rate": 0.9}),
                   replications=max(10, advice.required_n or 10), seed=2026).run()
study.save("results/mm1-study")
again = ExperimentResult.load("results/mm1-study")
md = again.metadata
print(md.experiment_id, md.model_name, md.model_version, md.seed, md.git, md.environment["python"])
rec = again.records[0]
replay = model.simulate(rec.parameters, seed=rec.seed)       # any replication, re-run alone
assert replay.metrics["resource.server.wait.mean"] == rec.metrics["resource.server.wait.mean"]

What the statistics do and do not tell you

• Replications are independent. Different replications use independent seeds, so t-intervals across replications are valid (with the usual normal-approximation caveats for very skewed metrics with few replications - use more replications or bootstrap_ci).
• Observations within one run are not independent. Do not compute a t-interval over all customers' waiting times in one run; use replication means, or batch_means for a single long run.
• Intervals cover sampling error only. They measure how precisely the model's expected output has been estimated, not how well the model matches reality. Validation is a separate, domain-specific task.
• Comparisons are not causal claims about the world. A significant difference means the model responds to the change.
• Many comparisons inflate false positives. With 20 metrics at 95%, expect about one spurious "significant" difference. Pre-register the metrics that matter, or adjust: compare(..., adjust="holm") (simulsi analyze --adjust holm).
• Optimising over noise flatters the winner. Re-evaluate the selected design with fresh seeds before reporting it (optimize(..., indifference=...) does this with select_best).

Uncertainty from the input data

Simulation intervals measure simulation noise. If the inputs were fitted to a small data set, the parameters themselves are uncertain, and that can matter more. input_uncertainty bootstraps the data, refits, re-runs and splits the variance into simulation noise and input uncertainty:

import numpy as np

from simulsi.analysis import InputData, input_uncertainty
from simulsi.models import mmc

gaps = np.random.default_rng(1).exponential(1 / 0.7, size=40)   # 40 observed interarrival times
res = input_uncertainty(
    mmc.with_options(duration=2_000, warmup=200),
    {"arrival_rate": InputData(gaps, "exponential", lambda d: 1 / d.mean)},
    "resource.server.wait.mean",
    n_bootstrap=10, replications=3,
)
print(res.format())   # with 40 data points, most of the uncertainty is in the input

When input_share is high, collecting more data narrows the answer far more than running more replications.

Rare events

rare_event_probability(model, metric, threshold) estimates the chance of an extreme outcome with an exact (Clopper-Pearson) interval and the number of replications needed for a 10% relative error. With enough runs it also fits a generalised Pareto distribution to the upper tail and extrapolates beyond the observed data (gpd_tail_probability), with a bootstrap interval:

from simulsi.analysis import rare_event_probability

risk = rare_event_probability(mmc.with_options(duration=2_000), "resource.server.wait.max", 40.0,
                              params={"arrival_rate": 0.8}, replications=60)
print(risk.format())

Tail extrapolation assumes the tail keeps the shape seen in the data; check it against a longer run before relying on it.

Design of experiments

• grid(...) / Experiment.grid(...) - full factorial designs.
• one_at_a_time(...) scenarios and the analysis.one_at_a_time method.
• monte_carlo(..., sampling="lhs") - Latin hypercube designs over uncertain inputs, followed by correlation_sensitivity (Spearman, SRC with R^2).
• finite_difference(...) - local gradients with common random numbers.
• morris_screening(...) - cheap global screening of many inputs, then sobol_indices(...) on the ones that matter.
• select_best(...) - choose among alternatives with a probability of correct selection guarantee (procedure KN).
• fit_surrogate(...) / optimize(...) - Latin hypercube designs, Gaussian process metamodels, Bayesian optimisation and the cross-entropy method.
• pareto_search(...) - trade-off fronts between several objectives.

Reporting checklist

• model version and git commit, SimulSI version (result.metadata)
• run length, warm-up and how they were chosen
• number of replications and why it is enough
• confidence level and interval method (paired-t / Welch / bootstrap)
• validation performed and its limits
• all scenario definitions (result.metadata.scenarios)

Reading list

• A. M. Law, Simulation Modeling and Analysis, 5th ed., McGraw-Hill, 2015 - output analysis, replication sizing (the approximation used by required_replications), common random numbers.
• J. Banks et al., Discrete-Event System Simulation, 5th ed., Pearson, 2010.
• B. L. Nelson, Foundations and Methods of Stochastic Simulation, 2nd ed., Springer, 2021.
• A. Saltelli et al., Global Sensitivity Analysis: The Primer, Wiley, 2008.
• M. D. Morris, "Factorial sampling plans for preliminary computational experiments", Technometrics 33(2), 1991; F. Campolongo, J. Cariboni and A. Saltelli, "An effective screening design for sensitivity analysis of large models", Environmental Modelling & Software 22(10), 2007.
• S.-H. Kim and B. L. Nelson, "A fully sequential procedure for indifference-zone selection in simulation", ACM TOMACS 11(3), 2001.
• R. R. Barton and L. W. Schruben, "Resampling methods for input modeling", Proceedings of the Winter Simulation Conference, 2001; R. Y. Rubinstein and D. P. Kroese, The Cross-Entropy Method, Springer, 2004; S. Coles, An Introduction to Statistical Modeling of Extreme Values, Springer, 2001.
• C. E. Rasmussen and C. K. I. Williams, Gaussian Processes for Machine Learning, MIT Press, 2006; D. R. Jones, M. Schonlau and W. J. Welch, "Efficient global optimization of expensive black-box functions", Journal of Global Optimization 13, 1998.