import numpy as np
import plotly.io as pio
import simpy
from sim_tools.distributions import Exponential, Lognormal
from vidigi.logging import EventLogger, TrialLogger
from vidigi.resources import VidigiStore
pio.renderers.default = "notebook"Feature Example: Choosing a Warm-up Length with Welch’s Procedure
A simulation that starts empty is not yet representative of the steady-state system you usually want to study - the first entities see an artificially short queue, because nobody was there before them. Discarding this warm-up period is standard practice, but “how long is it?” is not something you can just guess: too short, and the bias is still there; too long, and you throw away good data for nothing.
Welch’s (1983) moving-average procedure gives a visual way to answer this: run several replications, average them at each point in time, smooth the result, and look for where the curve stops drifting. vidigi implements this as vidigi.analysis.welch_moving_average() and vidigi.plots.plot_warm_up_diagnostic() (also available as TrialLogger.plot_warm_up_diagnostic()).
This notebook reuses the single-queue clinic model from feat_event_logger.ipynb’s “Simple Example” section - see that notebook for the full walkthrough of EventLogger/VidigiStore. The only thing that changes here is the arrival rate, turned up to make the warm-up transient visible.
References
- Welch, P. D. (1983). The statistical analysis of simulation results. In S. S. Lavenberg (Ed.), The Computer Performance Modeling Handbook (pp. 268-328). Academic Press. The original source for the procedure implemented here - a book chapter with no DOI or stable free-standing URL, being from before either existed.
- Law, A. M. Simulation Modeling and Analysis (McGraw-Hill). The secondary source this implementation’s docstrings cite for the exact edge-of-window formula - see the “Output Data Analysis” chapter in any edition. Publisher/author page, for anyone wanting to check the formula against the original working.
- Rossetti, M. D. Simulation Modeling and Arena - a freely-readable online textbook covering the same procedure with a runnable worked example: rossetti.github.io/RossettiArenaBook (see the statistical output analysis chapters).
- Robinson, S. (2004). Simulation: The Practice of Model Development and Use. Wiley. Cited by the DES RAP book (below) for the general “look for where the curve stabilises” principle behind time series inspection.
- Heather, A., Monks, T., Harper, A., Alidoost, F., Challen, R., Slater, T., & Mustafee, N. (2026). Reproducible analytical pipelines for healthcare discrete-event simulation: An open guide and worked examples. NIHR Open Research, 6, 68. doi.org/10.3310/nihropenres.14296.1. The “DES RAP book” - its warm-up length page demonstrates the cumulative-mean technique
method="cumulative"(below) is modelled on.
Model setup
Show the global parameter class code
class g:
"""
Create a scenario to parameterise the simulation model
Parameters:
-----------
random_number_set: int
Set to control the initial seeds of each stream of pseudo
random numbers used in the model.
n_cubicles: int
The number of treatment cubicles
treat_mean, treat_var: float
Mean and variance of the treatment duration distribution (Lognormal)
arrival_rate: float
Mean of the exponential inter-arrival time distribution
sim_duration: int
The number of time units the simulation will run for
number_of_runs: int
The number of replications
"""
random_number_set = 42
n_cubicles = 4
treat_mean = 25
treat_var = 5
# NOT the same value as feat_event_logger.ipynb's original (arrival_rate=5)
# deliberately: with that arrival rate, this system's utilisation is
# rho = (1/5) / (4/40) = 2.0 - overloaded, so the queue never reaches a
# steady state to diagnose at all, it just grows without bound. 8 gives
# rho = (1/8) / (4/25) = 0.78 - busy enough for a real transient, stable
# enough to actually settle. Pushing this further towards 1 makes the
# transient slower and considerably noisier - a real property of
# near-capacity queues (their "relaxation time" grows sharply as
# utilisation approaches 100%), not a bug in the diagnostic. Try it.
arrival_rate = 8
sim_duration = 3000
number_of_runs = 20Show the patient class code
class Patient:
"""Class defining details for a patient entity"""
def __init__(self, p_id):
self.id = p_idShow the model code
# Class representing our model of the clinic - identical in structure to
# feat_event_logger.ipynb's "Simple Example", only the g.arrival_rate above differs.
class Model:
def __init__(self, run_number):
self.env = simpy.Environment()
self.run_number = run_number
self.logger = EventLogger(env=self.env, run_number=self.run_number)
self.patient_counter = 0
self.init_distributions()
self.init_resources()
def init_distributions(self):
self.patient_inter_arrival_dist = Exponential(
mean=g.arrival_rate, random_seed=self.run_number * g.random_number_set
)
self.treat_dist = Lognormal(
mean=g.treat_mean,
stdev=g.treat_var,
random_seed=self.run_number * g.random_number_set,
)
def init_resources(self):
self.treatment_cubicles = VidigiStore(
self.env,
num_resources=g.n_cubicles,
label="treatment_cubicle",
logger=self.logger,
)
def generator_patient_arrivals(self):
while True:
self.patient_counter += 1
p = Patient(self.patient_counter)
self.env.process(self.attend_clinic(p))
yield self.env.timeout(self.patient_inter_arrival_dist.sample())
def attend_clinic(self, patient):
self.logger.log_arrival(entity_id=patient.id)
self.logger.log_queue(entity_id=patient.id, event="treatment_wait_begins")
with self.treatment_cubicles.request(
entity_id=patient.id,
start_event="treatment_begins",
end_event="treatment_complete",
) as req:
yield req
yield self.env.timeout(self.treat_dist.sample())
self.logger.log_departure(entity_id=patient.id)
def run(self):
self.env.process(self.generator_patient_arrivals())
self.env.run(until=g.sim_duration)event_logs = []
for run in range(1, g.number_of_runs + 1):
my_model = Model(run)
my_model.run()
event_logs.append(my_model.logger)
trial_logs = TrialLogger(event_logs)
trial_logs.summary(){'number_of_runs': 20, 'label': None, 'scenario_attached': False}Reading the warm-up length off a plot
plot_warm_up_diagnostic() needs to know which series to diagnose. series="queue" looks at the length of a named queue at regular snapshots, ensemble-averaged across every replication and then smoothed. windows= overlays a few different smoothing widths - if they broadly agree on where the curve flattens, that agreement is the signal Welch’s procedure is looking for.
fig = trial_logs.plot_warm_up_diagnostic(
series="queue",
event="treatment_wait_begins",
every_x_time_units=25,
windows=(5, 10, 20),
)
fig.update_layout(title="Welch's procedure: treatment queue length")
fig.show()The raw ensemble mean (the thin dotted line) is too noisy to read by eye - that noise is exactly what the smoothing is for. All three window= curves agree on the sharp part: the queue climbs from empty to around 1 within the first few hundred time units. What they don’t agree on so cleanly is what happens after - the smoothed curves keep wandering by roughly +/-0.2-0.3 (a quarter or so of the level itself) well past t=1000, rather than settling onto a flat line. That’s a realistic outcome, not a failure of the method: at ~78% utilisation this queue is autocorrelated enough, and 20 replications is little enough, that Welch’s procedure here narrows the question to “somewhere in the first few hundred time units” rather than answering it with a single confident number. There is deliberately no automatic answer for exactly this reason: Welch’s procedure is a visual one, and a mechanical flatness threshold would report false confidence on a series that wanders like this one does. We’ll pick warm_up=500 below and then check how much that specific choice actually matters.
Checking the claims above against the real implementation
Rather than take the description above on trust, the cell below prints the actual welch_moving_average() source - pulled live from the installed vidigi package via inspect.getsource, not pasted in and liable to drift out of sync with the real code. _ensemble_mean is the private helper it (and plot_warm_up_diagnostic) shares for the replication-averaging step.
Show the welch_moving_average source, read live from the installed package
import inspect
from vidigi.analysis import _ensemble_mean, welch_moving_average
print(inspect.getsource(_ensemble_mean))
print(inspect.getsource(welch_moving_average))def _ensemble_mean(series_by_run: Sequence[Sequence[float]]):
"""Truncate every run to the shortest length and average across runs.
Shared by `welch_moving_average` and `plot_warm_up_diagnostic`, so the
latter can draw the raw ensemble-average trace and feed the *same*
length-matched runs onward to `welch_moving_average` for each window -
truncating twice would otherwise warn twice for the same unequal-length
input.
Returns
-------
tuple
`(m, ensemble_mean, trimmed)`: the common length, the length-`m`
array of per-index means, and the list of length-`m` per-run arrays
it was computed from.
"""
lengths = [len(run) for run in series_by_run]
m = min(lengths)
if max(lengths) != m:
warnings.warn(
f"Runs of unequal length were passed (lengths {lengths}); every "
f"run is truncated to the shortest, {m}.",
UserWarning,
stacklevel=3,
)
trimmed = [np.asarray(run[:m], dtype=float) for run in series_by_run]
ensemble_mean = np.mean(np.array(trimmed), axis=0)
return m, ensemble_mean, trimmed
def welch_moving_average(
series_by_run: Sequence[Sequence[float]],
window: int | None = None,
*,
method: WarmUpMethod = "welch",
) -> np.ndarray:
"""
Smooth an ensemble-averaged series to visually locate a warm-up length.
Underlies `plot_warm_up_diagnostic`. Takes one series per replication -
e.g. queue length at each snapshot, resource occupancy at each snapshot,
or per-entity durations in arrival order - and reduces the replications
to a single smoothed curve whose early instability (or lack of it) is the
visual signal a modeller reads off to choose `warm_up=`.
Parameters
----------
series_by_run : sequence of sequence of float
One series per replication, each already in the order its warm-up
length will be judged against (snapshot time order, or entity
arrival order). Runs of unequal length are truncated to the
shortest, with a warning - the ensemble average at index i needs a
value from every run at that index, so a longer run contributes no
extra information past the shortest run's end.
window : int, optional
Half-width of the moving-average window. Required, and must be a
positive integer less than the (truncated) series length, when
`method="welch"`. Ignored when `method="cumulative"` or
`method="none"`, neither of which has a window.
method : {"welch", "cumulative", "none"}, default="welch"
- `"welch"`: Welch's (1983) moving-average procedure, as described
in Law, *Simulation Modeling and Analysis*. The ensemble average is
smoothed with a symmetric window of full width `2 * window + 1`;
for the first `window` points, where a full-width window would run
off the start of the series, it shrinks to the narrower symmetric
window that does exist (`2i - 1` points, i.e. only the `i - 1`
points actually available on each side - `i` here is Welch's own
1-indexed position in the series, `i = 1` being the first point,
not this function's 0-indexed output array) rather than returning
`NaN` there - this is the detail a `pandas.rolling(center=True)`
call gets wrong. A window this wide also has no full-width
neighbourhood in the last `window` points of the series, so those
are not returned at all: the output has length
`len(series) - window`, not `len(series)` - a larger `window`
gives a smoother curve at the cost of a shorter one.
- `"cumulative"`: the plain expanding mean of the ensemble average,
`mean(ensemble[:i + 1])` for every `i` - simpler and needs no
window, but each new point only shifts the running mean by `1/i`,
so a biased early transient decays out of it far more slowly than
out of Welch's fixed-width window. The risk is not that this curve
looks rougher - it is *smoother*, in a way that can make it look
settled long after (or hide a real shift long before) the series
has actually stabilised. This is the "cumulative mean" time series
inspection technique demonstrated in the DES RAP book (Heather et
al., 2026 -
https://pythonhealthdatascience.github.io/des_rap_book/pages/guide/output_analysis/length_warmup.html),
which plots both the pooled cumulative mean and each replication's
own cumulative mean (matching `show_runs=True` below) and cites
Robinson, S. (2004), *Simulation: The Practice of Model Development
and Use* (Wiley), for the general "look for where the curve
stabilises" principle. Returned at the series' full length.
- `"none"`: the ensemble average itself, completely unsmoothed - no
window, no running mean, nothing. A further step beyond the
cumulative-mean inspection above - here dropping the running-mean
smoothing entirely rather than only the fixed window - not itself
the specific technique the DES RAP book demonstrates. Noisier than
either other method by construction, since nothing here reduces
the within-replication variance the ensemble average didn't
already remove; useful specifically when you want to see that
noise rather than have it smoothed away. Returned at the series'
full length, identical to what `_ensemble_mean` computes.
Returns
-------
numpy.ndarray
The smoothed series. See `method` above for its length.
Raises
------
ValueError
If `series_by_run` is empty; if `method` is not `"welch"`,
`"cumulative"` or `"none"`; or if `method="welch"` and `window` is
`None`, not a positive integer, or `>=` the (truncated) series
length, leaving nothing to return.
Notes
-----
Deliberately has no automatic warm-up-length selector. Welch's procedure
is explicitly a visual one - a flatness threshold returns a confident
number that is wrong on any series with slow drift, silently discarding
the wrong amount of data with no signal anything went awry. Overlay
several `window` values (see `plot_warm_up_diagnostic`) and read off the
point where they agree, by eye.
See Also
--------
plot_warm_up_diagnostic : Builds the figure this function's output is drawn into.
"""
if len(series_by_run) == 0:
raise ValueError("`series_by_run` must contain at least one run.")
if method not in ("welch", "cumulative", "none"):
raise ValueError(
f"`method` must be 'welch', 'cumulative' or 'none'; got {method!r}."
)
m, ensemble_mean, _ = _ensemble_mean(series_by_run)
if method == "none":
return ensemble_mean
if method == "cumulative":
return np.cumsum(ensemble_mean) / np.arange(1, m + 1)
if window is None or window < 1:
raise ValueError(
f"`window` must be a positive integer when method='welch'; got {window!r}."
)
if window >= m:
raise ValueError(
f"`window` ({window}) leaves nothing to return from a series of "
f"length {m} - Welch's procedure only defines a moving average up "
f"to `len(series) - window` points. Pass a smaller `window`."
)
output_length = m - window
smoothed = np.empty(output_length)
for i in range(output_length):
if i < window:
smoothed[i] = ensemble_mean[: 2 * i + 1].mean()
else:
smoothed[i] = ensemble_mean[i - window : i + window + 1].mean()
return smoothed
method="cumulative": a simpler, but riskier, alternative
method="cumulative" draws the plain running mean of the ensemble average instead - no window to choose, but every new point only nudges it by 1/i. That is not simply “rougher” - if anything the curve below looks smoother than the Welch curves above, which is exactly the problem: a slow-decaying early bias can leave it looking settled long before the queue has actually reached its steady-state behaviour. This is the “time series inspection” technique the DES RAP book (Heather et al., 2026) demonstrates: plot the cumulative mean - both pooled and per replication - and look for where it settles.
fig = trial_logs.plot_warm_up_diagnostic(
series="queue",
event="treatment_wait_begins",
every_x_time_units=25,
method="cumulative",
)
fig.update_layout(title="Cumulative mean: treatment queue length")
fig.show()method="none": no smoothing at all
A third option skips smoothing entirely and plots the ensemble average exactly as computed - one step further than the cumulative-mean view above, dropping the running mean too rather than just Welch’s fixed window. It is by construction the noisiest of the three views here - nothing about it reduces the within-replication variance the ensemble average didn’t already remove - which is also exactly its appeal: nothing about the underlying data is hidden from you by a choice of window or averaging method.
fig = trial_logs.plot_warm_up_diagnostic(
series="queue",
event="treatment_wait_begins",
every_x_time_units=25,
method="none",
)
fig.update_layout(title='method="none": raw ensemble-averaged queue length')
fig.show()The DES RAP book’s own figures go one step further and plot every individual replication’s own trace alongside the pooled line, rather than just the pooled line on its own - though there, each replication’s trace is itself a cumulative mean (matching method="cumulative" above), not the fully raw series show_runs=True draws here. show_runs=True reproduces the same idea of overlaying every run: every run drawn faintly in the background, under one shared legend entry (“individual runs” - a legend entry per run would swamp the windows=/method entries that are the actual point of the rest of this plot, once there are twenty of them). It works with any method, not just "none", since it’s answering a different question - how much do the replications actually disagree with each other? - from whichever smoothing method is chosen.
fig = trial_logs.plot_warm_up_diagnostic(
series="queue",
event="treatment_wait_begins",
every_x_time_units=25,
method="none",
show_runs=True,
)
fig.update_layout(title='method="none", show_runs=True: every replication, faintly')
fig.show()Occupancy converges faster than the queue does
series="occupancy" diagnoses how many treatment cubicles are busy, rather than how many patients are waiting. It is worth comparing the two: with four cubicles and patients arriving faster than they’re treated on average, the cubicles fill up almost as soon as the simulation starts, while the queue behind them takes much longer to settle into its steady-state distribution. A warm-up length chosen from occupancy alone would badly understate what the queue actually needs.
fig = trial_logs.plot_warm_up_diagnostic(
series="occupancy",
event="treatment_begins",
every_x_time_units=25,
windows=(5, 10, 20),
)
fig.update_layout(title="Welch's procedure: cubicles occupied")
fig.show()/home/runner/work/vidigi/vidigi/src/vidigi/analysis.py:2165: UserWarning: 66 resource_use row(s) were still open at the end of the window - censored at 2999.8678663039827 (unclosed='censor', the default). Pass unclosed='drop' to exclude them instead.
intervals = resource_use_intervals(Diagnosing a duration instead of a time series
series="duration" diagnoses per-entity durations - here, how long each patient waited for a cubicle - in arrival order rather than simulated time, since there is no shared time grid to average durations against. The x-axis is “the n-th patient to arrive”, not a time.
fig = trial_logs.plot_warm_up_diagnostic(
series="duration",
first_event="treatment_wait_begins",
second_event="treatment_begins",
windows=(20, 50),
)
fig.update_layout(title="Welch's procedure: individual waiting times, in arrival order")
fig.show()/home/runner/work/vidigi/vidigi/src/vidigi/logging.py:2726: UserWarning: Runs of unequal length were passed (lengths [380, 364, 374, 387, 350, 380, 378, 357, 375, 393, 354, 362, 365, 373, 374, 381, 400, 396, 367, 365]); every run is truncated to the shortest, 350.
return _plot_warm_up_diagnostic(Applying the reading
We’ll use warm_up=500 - past the sharp initial rise, though (per the previous section) not obviously past all of the drifting. Every warm-up-aware function in vidigi.analysis (and the matching TrialLogger methods) takes it as warm_up=. Given the ambiguity above, it’s worth checking how much this specific choice actually matters before trusting it - not just comparing “with” against “without”.
candidate_warm_ups = [0, 300, 500, 700]
for candidate in candidate_warm_ups:
mean_wait = trial_logs.get_event_duration_stat(
"treatment_wait_begins", "treatment_begins", what="mean", warm_up=candidate
)
print(f"warm_up={candidate:>4}: mean wait = {mean_wait:.2f}")
warm_up = 500warm_up= 0: mean wait = 8.32
warm_up= 300: mean wait = 8.55
warm_up= 500: mean wait = 8.41
warm_up= 700: mean wait = 8.43warm_up= on get_event_duration_stat excludes a pairing by when it started - a patient who began waiting before warm_up is dropped even if they were still waiting when the cutoff passed, since a wait time is one observation, not something that can be partially inside the window. That’s a genuinely different rule from resource_use-based warm_up= (below), which censors a bout rather than excluding it outright - see event_durations’s docstring for the full reasoning.
The same warm_up= is available on plot_duration_distribution and plot_metric_bar. Unlike get_event_duration_stat above (which pools every entity from every run into one number), across="runs" computes the mean within each run first and then averages those - a different, and for this comparison more useful, statistic, because it comes with a confidence interval:
fig_before = trial_logs.plot_metric_bar(
[
{
"label": "Including startup",
"first_event": "treatment_wait_begins",
"second_event": "treatment_begins",
}
],
across="runs",
error_bars="ci",
title="Mean wait, before discarding warm-up",
)
fig_before.show()
mean_before, ci_before = fig_before.data[0].y[0], fig_before.data[0].error_y.array[0]
print(f"Mean wait, including startup: {mean_before:.2f} +/- {ci_before:.2f} (95% CI)")Mean wait, including startup: 8.24 +/- 1.44 (95% CI)fig_after = trial_logs.plot_metric_bar(
[
{
"label": "After warm-up",
"first_event": "treatment_wait_begins",
"second_event": "treatment_begins",
}
],
across="runs",
error_bars="ci",
warm_up=warm_up,
title="Mean wait, after discarding warm-up",
)
fig_after.show()
mean_after, ci_after = fig_after.data[0].y[0], fig_after.data[0].error_y.array[0]
print(
f"Mean wait, after warm_up={warm_up}: {mean_after:.2f} +/- {ci_after:.2f} (95% CI)"
)Mean wait, after warm_up=500: 8.32 +/- 1.59 (95% CI)A word of caution
The two confidence intervals overlap almost completely - discarding warm-up here moved the mean by less than the noise between replications. That is a genuine result, not a failure of the method: at this system’s ~78% utilisation the startup bias is real but small, and discarding data to remove it also cost precision (the CI got wider, not narrower - fewer observations went into each run’s mean). Whether that trade is worth making depends on how large the bias actually is relative to the noise, which is exactly what comparing the two intervals - not just the two point estimates - tells you. On a more heavily-loaded system (or a longer warm-up transient) the bias would be larger relative to the noise, and the case for discarding it stronger.
This notebook stops at the visual Welch/cumulative-mean approach vidigi implements. Related questions it doesn’t cover - how autocorrelated these snapshots are, whether 20 replications is enough to trust a Welch plot, how to choose window= systematically, and alternative approaches such as batch means or automated rules like MSER - are covered in Law’s and Rossetti’s books (see References above).
resource_utilisation and queue_size_over_time already had their own warm_up= (added earlier in 2.0.0) - trimming the snapshot window itself, rather than filtering individual observations, since presence at a snapshot has to be worked out from the whole log, not just the rows after the cutoff:
trial_logs.get_resource_utilisation(
by="resource", resource_col_name="unique_resource_id", warm_up=warm_up
).head()/home/runner/work/vidigi/vidigi/src/vidigi/analysis.py:1945: UserWarning: 66 resource_use row(s) were still open at the end of the window - censored at 2999.8678663039827 (unclosed='censor', the default). Pass unclosed='drop' to exclude them instead.
intervals = resource_use_intervals(| run_number | resource_id | busy_time | mean_in_use | capacity | utilisation | |
|---|---|---|---|---|---|---|
| 0 | 1 | treatment_cubicle_1 | 1960.646548 | 0.784300 | 1.0 | 0.784300 |
| 1 | 1 | treatment_cubicle_2 | 1983.439444 | 0.793418 | 1.0 | 0.793418 |
| 2 | 1 | treatment_cubicle_3 | 1982.435127 | 0.793016 | 1.0 | 0.793016 |
| 3 | 1 | treatment_cubicle_4 | 1918.352165 | 0.767381 | 1.0 | 0.767381 |
| 4 | 2 | treatment_cubicle_1 | 1958.160425 | 0.783306 | 1.0 | 0.783306 |
trial_logs.plot_queue_size(
["treatment_wait_begins"],
every_x_time_units=25,
limit_duration=g.sim_duration,
warm_up=warm_up,
)welch_moving_average() directly
plot_warm_up_diagnostic() is a thin wrapper: it builds a per-run series from the event log, then hands it to vidigi.analysis.welch_moving_average(), which does the actual ensemble-averaging and smoothing. Calling it directly is useful if you already have a series from somewhere else - a KPI computed outside vidigi, say - and just want the smoothing.
from vidigi.analysis import welch_moving_average
# Eight short, noisy "runs" of the same underlying series
rng = np.random.default_rng(42)
runs = [rng.normal(loc=np.linspace(0, 10, 40), scale=1.5) for _ in range(8)]
smoothed = welch_moving_average(runs, window=5, method="welch")
print(
f"{len(runs[0])} points in -> {len(smoothed)} points out (window=5 costs 5 off the end)"
)
print(np.round(smoothed, 2))40 points in -> 35 points out (window=5 costs 5 off the end)
[-0.1 0.17 0.41 0.65 0.97 1.21 1.46 1.78 1.94 2.09 2.35 2.69
2.95 3.21 3.46 3.78 3.98 4.24 4.48 4.81 5.18 5.49 5.74 6.
6.29 6.55 6.75 7.01 7.32 7.59 7.88 8.07 8.27 8.53 8.72]That covers choosing where to cut the warm-up period. The related question of how many replications are enough to trust the result - covered by confidence-interval-based replication analysis - is covered in feat_replication_analysis.ipynb, which reuses this same clinic model and 20-replication trial.
A third, related question - applying an already-chosen warm-up length to an animation, where naively filtering the event log by time silently corrupts the result rather than just excluding early data - is covered separately in feat_animation_warm_up.ipynb.