Feature Example: How Many Replications Are Enough?

Warning

This example has not yet been updated to reflect all of the new features and recommendations from vidigi 2.0.0, which have further simplified the process of adding vidigi to your model and accessing and modifying your animation.

For now, all of the code below will still work, but check out the getting started page for a full guide to the recommended way to use vidigi 2.0.0.

A single run of a stochastic simulation tells you what happened once. Averaging several replications gets you closer to the true steady-state answer, but “several” is not a number until you ask how much the estimate is still moving around as you add more runs. Too few replications and a reported mean is barely better than a guess; too many and you have spent compute on precision nobody needed.

vidigi.analysis.replication_precision() and vidigi.plots.plot_replication_analysis() (also available as TrialLogger.get_replication_precision() / .plot_replication_analysis()) answer this with a confidence-interval-based diagnostic: recompute the interval after each replication, watch its width relative to the mean shrink as more runs are added, and read off the point after which it stays acceptably tight.

This is the companion to feat_warm_up.ipynb, which covers where to cut the warm-up period - this notebook covers how many replications to run once that cut is decided. It reuses the same single-queue clinic model and the same 20-replication trial, so the two notebooks’ numbers are directly comparable.

References

• Hoad, K., Robinson, S., & Davies, R. (2010). Automated selection of the number of replications for a discrete-event simulation. Journal of the Operational Research Society, 61(11), 1632-1644. doi.org/10.1057/jors.2009.121. The source for the confidence-interval-based replication procedure implemented here: run replications until the interval’s relative half-width stays under a chosen threshold.
• Law, A. M. Simulation Modeling and Analysis (McGraw-Hill) - see the “Output Data Analysis” chapter in any edition for the underlying confidence-interval-for-a-mean theory (vidigi.analysis.mean_confidence_interval) this diagnostic is built on. Publisher/author page.
• Rossetti, M. D. Simulation Modeling and Arena - a freely-readable online textbook with its own worked replication-count example: rossetti.github.io/RossettiArenaBook (see the statistical output analysis chapters).

import random

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"

Model setup

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

    arrival_rate = 8

    sim_duration = 3000
    number_of_runs = 20

class Patient:
    """Class defining details for a patient entity"""

    def __init__(self, p_id):
        self.id = p_id

# Identical to feat_warm_up.ipynb's model - same clinic, same
# parameters, so the two notebooks' numbers are directly comparable.
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"
        )

    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() as req:
            treatment_resource = yield req
            self.logger.log_resource_use_start(
                entity_id=patient.id,
                event="treatment_begins",
                resource_id=treatment_resource.id,
                unique_resource_id=treatment_resource.unique_id,
            )
            yield self.env.timeout(self.treat_dist.sample())
            self.logger.log_resource_use_end(
                entity_id=patient.id,
                event="treatment_complete",
                resource_id=treatment_resource.id,
                unique_resource_id=treatment_resource.unique_id,
            )
        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)

class Trial:
    def __init__(self):
        self.all_event_logs = []
        self.run_trial()

    def run_trial(self):
        for run in range(1, g.number_of_runs + 1):
            random.seed(run)
            my_model = Model(run)
            my_model.run()
            self.all_event_logs.append(my_model.logger)

clinic_trial = Trial()
trial_logs = TrialLogger(clinic_trial.all_event_logs)
trial_logs.summary()

{'number_of_runs': 20, 'label': None, 'scenario_attached': False}

Precision as replications accumulate

plot_replication_analysis() needs an event pair, just like plot_metric_bar - here, waiting time from treatment_wait_begins to treatment_begins. For each replication count k = 1..20 (in the order the runs were generated), it recomputes the mean and its 95% confidence interval using only the first k replications, then plots both the cumulative mean (with its CI band) and the interval’s relative half-width - deviation = half_width / mean - underneath.

fig = trial_logs.plot_replication_analysis("treatment_wait_begins", "treatment_begins")
fig.update_layout(width=900, height=650)
fig.show()

The top panel’s confidence band is wide with few replications and narrows as more are added - the visual signature of 1/sqrt(n) convergence. The bottom panel turns that into a single number to threshold on: the dashed line marks deviation_threshold (5% by default), and the title reports the smallest replication count after which deviation stays at or under it for every replication count from there to the end - not just the first one that happens to dip below by chance, which a noisy early curve can do long before the interval has actually settled.

Here, the title reports that deviation never stays below 5% within the 20 replications this trial has. That is a genuine, useful answer, not a failed one: it means waiting time at this system’s ~78% utilisation is variable enough between replications that 20 runs is not yet enough to pin the mean down to within 5% - a finding consistent with feat_warm_up.ipynb’s own observation that this queue is autocorrelated enough for 20 replications to be “little enough”. The honest response is to run more replications, not to treat whichever number happened to look lowest as good enough.

Trying more replications

20 replications wasn’t enough for this metric at this system’s utilisation - the natural next question is how many would be. Rather than guess, run more and check again: same model, same event pair, only g.number_of_runs changes.

One presentation wrinkle at this many points: plot_replication_analysis’s default marker size (tuned for the handful of replications a typical trial runs) overlaps into a thick, hard-to-read smear once there are hundreds of them on one axis. marker_size=/line_width= turn both down.

g.number_of_runs = 400

more_trial = Trial()
trial_logs_400 = TrialLogger(more_trial.all_event_logs)

fig = trial_logs_400.plot_replication_analysis(
    "treatment_wait_begins",
    "treatment_begins",
    marker_size=3,
    line_width=1.5,
)
fig.update_layout(width=900, height=650)
fig.show()

g.number_of_runs = 20  # restore for the rest of this notebook

At 400 replications, deviation does settle: the title reports a recommended replication count of 311, comfortably inside the 400 run here, and confirms it stays below 5% for every count from there to the end - not just the last one. Getting there took roughly sixteen times the 20 replications this notebook otherwise runs, which is itself the point: this queue’s ~78% utilisation makes waiting time variable enough between replications that a quick batch this small was never going to be enough on its own, consistent with feat_warm_up.ipynb’s own observation that this queue is autocorrelated enough for 20 replications to be “little enough”.

Two caveats on that 311, worth stating plainly rather than leaving implicit. First, it is a property of this random-number stream on this metric, not a fixed fact about the system - a different random_number_set, a different event pair, or a different run ordering would very likely settle at a different count, sometimes by a lot. Second, like the 20-replication analysis above it, it is still computed on the full 3000-time-unit run with no warm_up= trimming - the next section explains why that biases the mean the interval is drawn around. A tight interval at n=311 is not, on its own, a reason to skip that step; treat 311 as a demonstration that convergence is achievable for this metric, not as the replication count to actually run.

Everything from here on deliberately goes back to the original 20-replication trial_logs - still under-replicated by the standard just established. That’s intentional: the rest of this notebook is about how to read and use the diagnostic itself - the table it returns, what stays_below_threshold means, the caveats on trusting its number - not about re-establishing “enough” replications every time a new question is asked.

Checking the claims above against the real implementation

As with the warm-up notebook, the cell below prints the actual replication_precision() 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.

import inspect

from vidigi.analysis import replication_precision

print(inspect.getsource(replication_precision))

def replication_precision(
    values,
    *,
    ci_level: float = 0.95,
    deviation_threshold: float = 0.05,
) -> pd.DataFrame:
    """
    Running confidence-interval precision as replications accumulate.

    For k = 1..n, in the order `values` is given, computes the confidence
    interval using only the first k replications, plus its *relative*
    half-width (`deviation`) - the standard diagnostic (Hoad, Robinson &
    Davies, 2010) for deciding how many replications a study needs: run more
    until the interval is tight enough relative to the mean, not just "run
    some fixed number and hope."

    Parameters
    ----------
    values : array-like
        Per-replication values, **already in run order** - typically
        `replication_means(...)["value"]`. This function does not sort them;
        a cumulative diagnostic is only meaningful walked through in the
        order replications were actually generated.
    ci_level : float, default=0.95
        Confidence level for each cumulative interval - see
        `mean_confidence_interval`.
    deviation_threshold : float, default=0.05
        Relative half-width threshold used for `stays_below_threshold` (see
        *Returns*). `0.05` means the CI half-width must be within 5% of the
        cumulative mean.

    Returns
    -------
    pandas.DataFrame
        One row per k = 1..n, columns:

        - ``n_replications`` : int - k.
        - ``cumulative_mean`` : float - mean of `values[:k]`.
        - ``half_width``, ``lower``, ``upper`` : float - the confidence
          interval from the first k values. `NaN` at k=1 - a spread needs at
          least 2 points, matching `mean_confidence_interval`.
        - ``deviation`` : float - `half_width / abs(cumulative_mean)`, the
          relative precision. Always non-negative, so a metric with a
          negative mean (e.g. a before/after difference) is not read as
          trivially "precise" by a negative ratio. `NaN` wherever
          `half_width` is `NaN`, or if `cumulative_mean` is `0`.
        - ``stays_below_threshold`` : bool - True at row k if `deviation` is
          defined and no greater than `deviation_threshold` at k *and every
          later row, up to n*. Deliberately "stays below", not "first drops
          below": a noisy early curve can dip under the threshold once by
          chance and rise again, which would be a spurious recommendation.
          The smallest `n_replications` with `stays_below_threshold=True` is
          the recommended minimum replication count; always `False` at k=1.
          This is a property of the batch of `n` replications actually
          supplied, not a guarantee that deviation stays low forever - a run
          flagged `True` from a 20-replication batch could fail to qualify
          once replications 21+ are added and re-checked.

    Raises
    ------
    ValueError
        If `values` is empty.
    ImportError
        If `scipy` is not installed and `n >= 2` - see
        `mean_confidence_interval`.

    See Also
    --------
    mean_confidence_interval : The single-k confidence interval this is built from.
    replication_means : Produces the per-replication values this function consumes.
    vidigi.plots.plot_replication_analysis : Plots this table.

    Notes
    -----
    Recomputing a confidence interval at every k and reading off the first
    crossing of `deviation_threshold` risks "early convergence" - a spurious
    dip from a run of similar-by-chance values, not a genuinely settled
    interval. Hoad, Robinson & Davies (2010) name this problem directly and
    address it with a "look ahead": once precision first crosses the
    threshold, their algorithm runs a further, fixed number of replications
    (their `kLimit`, for which they recommend a default of 5) and checks
    precision stays crossed before accepting the result - shown empirically
    to remove the coverage failures a naive first-crossing rule produced in
    their own tests. `stays_below_threshold` here is a simpler stand-in for
    that same idea - checking all the way to the end of the supplied batch,
    not a fixed look-ahead window - rather than a literal implementation of
    `kLimit`, and, like the source paper's own procedure, validated only
    empirically rather than derived as a formal statistical correction for
    the underlying repeated-test problem. Treat the reported
    `stays_below_threshold`/recommended count as a starting point for
    judgement, not a fully-automated stopping rule - matching this
    package's deliberate choice not to fully automate
    `welch_moving_average`'s warm-up selection either.
    """
    series = pd.Series(values).reset_index(drop=True)
    n = len(series)
    if n == 0:
        raise ValueError("`values` must contain at least one replication.")

    rows = []
    for k in range(1, n + 1):
        window = series.iloc[:k]
        mean = window.mean()
        if k < 2:
            half_width = lower = upper = float("nan")
        else:
            ci = mean_confidence_interval(window, ci_level=ci_level)
            half_width, lower, upper = ci.half_width, ci.lower, ci.upper
        deviation = half_width / abs(mean) if mean != 0 else float("nan")
        rows.append(
            {
                "n_replications": k,
                "cumulative_mean": mean,
                "half_width": half_width,
                "lower": lower,
                "upper": upper,
                "deviation": deviation,
            }
        )

    result = pd.DataFrame(rows)
    suffix_max = result["deviation"][::-1].cummax()[::-1]
    result["stays_below_threshold"] = result["deviation"].notna() & (
        suffix_max <= deviation_threshold
    )
    return result

The table behind the plot

get_replication_precision() (or the free function vidigi.analysis.replication_precision, given replication_means(...)["value"]) returns the same per-k numbers the plot draws, for anyone who wants the table rather than the chart - to log it, threshold on it in code, or feed it into a report.

precision = trial_logs.get_replication_precision(
    "treatment_wait_begins", "treatment_begins"
)
precision

	n_replications	cumulative_mean	half_width	lower	upper	deviation	stays_below_threshold
0	1	9.896811	NaN	NaN	NaN	NaN	False
1	2	8.025390	23.778659	-15.753269	31.804050	2.962929	False
2	3	8.817961	5.765513	3.052448	14.583473	0.653837	False
3	4	9.792569	4.325842	5.466727	14.118411	0.441747	False
4	5	8.767955	4.079034	4.688920	12.846989	0.465221	False
5	6	9.038577	3.161068	5.877509	12.199645	0.349731	False
6	7	8.735946	2.648682	6.087264	11.384628	0.303193	False
7	8	8.359814	2.388462	5.971352	10.748276	0.285708	False
8	9	8.110806	2.132952	5.977853	10.243758	0.262977	False
9	10	8.895537	2.579493	6.316045	11.475030	0.289976	False
10	11	8.558402	2.417805	6.140597	10.976207	0.282507	False
11	12	8.297400	2.254657	6.042743	10.552057	0.271731	False
12	13	8.130054	2.085215	6.044839	10.215269	0.256482	False
13	14	8.197671	1.919757	6.277913	10.117428	0.234183	False
14	15	8.204198	1.774362	6.429837	9.978560	0.216275	False
15	16	8.087540	1.668079	6.419461	9.755619	0.206253	False
16	17	8.267246	1.604289	6.662956	9.871535	0.194054	False
17	18	8.458315	1.558382	6.899934	10.016697	0.184243	False
18	19	8.269444	1.520556	6.748888	9.790000	0.183876	False
19	20	8.236491	1.438757	6.797734	9.675248	0.174681	False

stays_below_threshold is False at k=1 (a single replication has no spread to estimate a confidence interval from at all) and again wherever a later replication count’s deviation rises back above the threshold - only a row where every subsequent row also qualifies is flagged True. The smallest n_replications with stays_below_threshold True is the number the plot’s title reports as “recommended”.

Reading the recommendation honestly

That recommendation is a statement about precision, not about correctness. Three things it does not protect against:

• Warm-up bias. replication_precision() was computed above on the full 3000-time-unit run, startup transient included - a tight interval around a biased mean is still wrong, just confidently so. Combine this with feat_warm_up.ipynb’s warm_up= first, then check replication precision on the post-warm-up statistic, as the cell below does.
• Autocorrelation between replications. The confidence interval here assumes each replication’s mean is an independent draw - true as long as every run uses its own, independent random number stream (as Model.init_distributions above does via random_seed=self.run_number * g.random_number_set), but not automatically true of every simulation setup.
• The recommendation is bounded by the batch you ran, not open-ended. stays_below_threshold only checks “below the threshold from here to the last replication you supplied” - a run flagged as converged from a 20-replication batch is not guaranteed to still qualify once replications 21+ are added and re-checked. Hoad, Robinson & Davies (2010) address this same “early convergence” risk with a “look ahead”: once precision first crosses the threshold, their algorithm runs a further, fixed number of replications (their kLimit, for which they recommend a default of 5) and checks it stays crossed before accepting the result - shown empirically in their own tests to fix the coverage failures a naive first-crossing rule produced. stays_below_threshold here is a simpler stand-in for that same idea - checking to the end of the batch rather than a fixed look-ahead window - not a literal implementation of kLimit, and, like the source paper’s own procedure, validated only empirically rather than derived as a formal statistical correction for the underlying repeated-test problem. Treat the reported recommendation as a starting point for judgement, the same way feat_warm_up.ipynb treats a Welch-diagnostic reading - not as a number the tool has decided for you.

warm_up = 500  # see feat_warm_up.ipynb for how this was chosen

precision_after_warm_up = trial_logs.get_replication_precision(
    "treatment_wait_begins", "treatment_begins", warm_up=warm_up
)
precision_after_warm_up

	n_replications	cumulative_mean	half_width	lower	upper	deviation	stays_below_threshold
0	1	7.811345	NaN	NaN	NaN	NaN	False
1	2	7.366359	5.654083	1.712276	13.020442	0.767555	False
2	3	8.595153	5.401394	3.193758	13.996547	0.628423	False
3	4	9.888387	4.991902	4.896485	14.880288	0.504825	False
4	5	8.717439	4.685018	4.032420	13.402457	0.537431	False
5	6	9.035387	3.634751	5.400635	12.670138	0.402280	False
6	7	8.765966	2.997529	5.768437	11.763495	0.341951	False
7	8	8.457946	2.612232	5.845714	11.070178	0.308849	False
8	9	8.185957	2.332568	5.853390	10.518525	0.284947	False
9	10	9.149343	2.989686	6.159657	12.139029	0.326765	False
10	11	8.803429	2.772877	6.030553	11.576306	0.314977	False
11	12	8.547113	2.563282	5.983831	11.110395	0.299900	False
12	13	8.310873	2.390203	5.920670	10.701077	0.287600	False
13	14	8.340041	2.195070	6.144971	10.535111	0.263197	False
14	15	8.326217	2.028977	6.297240	10.355194	0.243685	False
15	16	8.236357	1.895832	6.340525	10.132189	0.230178	False
16	17	8.340860	1.784980	6.555880	10.125840	0.214004	False
17	18	8.521457	1.717679	6.803778	10.239136	0.201571	False
18	19	8.295266	1.686258	6.609008	9.981524	0.203280	False
19	20	8.320213	1.594566	6.725646	9.914779	0.191650	False

Excluding the biased startup observations shifts the cumulative means (compare the cumulative_mean column above with the one two cells up) - a reminder that “how many replications” and “how much warm-up” are both answers a single number can only be trusted alongside, not instead of. deviation still never stays below the 5% threshold here either: warm-up and replication count are separate problems, and fixing one does not fix the other.

deviation_threshold= and ci_level=

Both are just numbers passed through to mean_confidence_interval and the “stays below” check - tightening deviation_threshold (wanting a more precise estimate) or widening ci_level (wanting more confidence in that estimate) both push the recommended replication count later. The cell above already showed 20 replications isn’t enough to reach the default 5% at 95% confidence for this metric; asking for 2% at 99% confidence only widens that gap, as the deviation curve below confirms.

# Tighter precision (2%) at a higher confidence level (99%) than the default call above.
fig = trial_logs.plot_replication_analysis(
    "treatment_wait_begins",
    "treatment_begins",
    deviation_threshold=0.02,
    ci_level=0.99,
)
fig.update_layout(width=900, height=650)
fig.show()

show_deviation=False

If the deviation panel isn’t needed - the recommendation in the title is often enough on its own - show_deviation=False drops it and plots only the cumulative mean and its CI band.

fig = trial_logs.plot_replication_analysis(
    "treatment_wait_begins",
    "treatment_begins",
    show_deviation=False,
)
fig.update_layout(width=900, height=450)
fig.show()

Closing note

This and feat_warm_up.ipynb together cover the two standard questions asked before trusting a stochastic simulation’s output: how much of each run to discard, and how many runs to run. vidigi implements the confidence-interval precision check Hoad, Robinson & Davies describe, but stops short of their paper’s own further step of automating the stopping decision - matching the deliberate choice already made for Welch’s procedure in the companion notebook: a stays_below_threshold column and a title to read, not a number the tool decides for you and runs with.

A third, related question - whether a metric drifts within a run’s steady operation, depending on when the entity that produced it arrived - is covered separately in feat_metric_vs_arrival_time.ipynb.