Chapter 3: Sampling Distributions and Random Streams

How Utils loads the fitted distribution files and samples from them with named, independently seeded random streams.

Written for: Intermediate Coder Confident Coder Decision-makers Domain experts General public

Chapter 2 ended with each replication getting its own seed from the master seed. This chapter picks that seed up inside the model and follows it into Utils.

What it is

Utils (in air_ambulance_des/utils.py) is the model’s sampling layer. Its constructor reads every fitted input the model needs from distribution_data/ (plus upper_allowable_time_bounds.csv from actual_data/) into pandas frames and dictionaries. After setup_seeds() it exposes selection methods (ampds_code_selection, care_category_selection, activity_time, inc_per_day_samples, sample_ad_hoc_reason and friends) that each return one sampled value.

The contract callers rely on:

• Input: a SeedSequence (the master_seed argument) and the files in distribution_data/, which Distribution Fitting produces.
• Output: one value per call, drawn from a stream named for the purpose, so the same seed always gives the same sequence for that purpose.
• Invariant: each sampling purpose has its own numpy Generator. Drawing more or fewer values for one purpose must not move the numbers any other purpose sees.
• Callers: DES_HEMS builds one Utils per run (Calls and Patient Journeys). HEMSAvailability receives that same object and calls sample_ad_hoc_reason during allocation (Resource Allocation).

The use case this chapter keeps returning to: a patient’s mobile time is sampled with activity_time("helicopter", "time_mobile"), and you want it to be reproducible and unaffected by someone adding a new sampled attribute elsewhere.

Utils also carries unrelated odds and ends: folder constants, school-holiday date calculators, and module-level helpers such as format_diff and get_rota_month_strings. They are out of scope here except where they bite (see Section 4).

How it works

Construction and seeding

DES_HEMS.__init__ turns the run’s integer seed into a SeedSequence, spawns one child for Utils, and calls setup_seeds() straight away.

air_ambulance_des/des_hems.py:61-71

        self.random_seed = random_seed
        self.random_seed_sequence = SeedSequence(self.random_seed)

        self.print_debug_messages = print_debug_messages

        self.utils = Utils(
            master_seed=self.random_seed_sequence.spawn(1)[0],
            print_debug_messages=self.print_debug_messages,
        )

        self.utils.setup_seeds()

So the seed chain from Chapter 2 is: master seed, then per-run seed (integer), then a SeedSequence per model, then one child per Utils, then one stream per purpose (Figure 1).

flowchart TD
    A["random_seed (int, per run)"] --> B["SeedSequence(random_seed)"]
    B -->|"spawn(1)[0]"| C["Utils.master_seed_sequence"]
    C -->|"spawn(22) in setup_seeds"| D["22 named Generators in Utils.rngs"]
    C -->|"spawn(n) in build_seeded_distributions"| E["one Generator per activity-time fit"]
    E --> F["SeededDistribution"]

Figure 1: Seed hierarchy inside one run

setup_seeds() does the spawning:

air_ambulance_des/utils.py:184-221

    def setup_seeds(self):
        #########################################################
        # Control generation of required random number streams  #
        #########################################################

        # Spawn substreams for major simulation modules
        module_keys = [
            "activity_times",
            "ampds_code_selection",
            "care_category_selection",
            "callsign_group_selection",
            "vehicle_type_selection",
            "hems_result_by_callsign_group_and_vehicle_type_selection",
            "hems_result_by_care_category_and_helicopter_benefit_selection",
            "hems_results_by_patient_outcome_and_quarter_and_vehicle_type_and_callsign_group_selection",
            "pt_outcome_selection",
            "sex_selection",
            "age_sampling",
            "calls_per_day",
            "calls_per_hour",
            "predetermine_call_arrival",
            "call_iat",
            "helicopter_benefit_from_reg",
            "hems_case",
            "ad_hoc_reason_selection",
            "know_heli_benefit",
            "helicopter_benefit_from_cc",
            "helicopter_benefit_from_ec",
            "know_cc_ec_benefit",
        ]
        # Efficiently spawn substreams
        spawned = self.master_seed_sequence.spawn(len(module_keys))

        # Map substreams and RNGs to keys
        self.seed_substreams = dict(zip(module_keys, spawned))
        self.rngs = {key: default_rng(ss) for key, ss in self.seed_substreams.items()}

        self.build_seeded_distributions(self.master_seed_sequence)

SeedSequence.spawn is stateful: it hands out children numbered from an internal counter. That is why the second spawn inside build_seeded_distributions (called on the last line above) gives the activity-time generators streams that differ from the 22 named ones, without any extra bookkeeping.

Weighted table lookups

Most selection methods follow one pattern: filter a probability table by the conditioning variables, then draw one row with proportion as weights, passing the named generator as random_state.

air_ambulance_des/utils.py:289-301

    def ampds_code_selection(self, hour: int) -> int:
        """
        This function will allocate and return an AMPDS card category
        based on the hour of day
        """

        df = self.hour_by_ampds_df[self.hour_by_ampds_df["hour"] == hour]

        return pd.Series.sample(
            df["ampds_card"],
            weights=df["proportion"],
            random_state=self.rngs["ampds_code_selection"],
        ).iloc[0]

care_category_selection, pt_outcome_selection, callsign_group_selection, vehicle_type_selection_qtr and the HEMS result sampler Utils.hems_results_by_patient_outcome_and_time_of_day_and_quarter_and_vehicle_type_and_callsign_group_probs are the same shape with different filters and different stream keys. Two methods deviate: sex_selection compares rngs["sex_selection"].uniform(0, 1) with a single probability, and sample_ad_hoc_reason uses Generator.choice with normalised probabilities.

Continuous distributions and activity times

The files activity_time_distributions.txt, age_distributions.txt and the two inc_per_day_* files hold Python literals, read with ast.literal_eval in the constructor. Activity times are special: build_seeded_distributions replaces self.activity_time_distr with a dictionary of SeededDistribution objects, keyed by (vehicle_type, time_type).

air_ambulance_des/utils.py:890-955

    def build_seeded_distributions(self, seed_seq):
        """
        Build a dictionary of seeded distributions keyed by (vehicle_type, time_type)

        Returns
        -------
        dict of (vehicle_type, time_type) -> SeededDistribution

        e.g. {('helicopter', 'time_allocation'): <utils.SeededDistribution object at 0x000001521627C750>,
        ('car', 'time_allocation'): <utils.SeededDistribution object at 0x000001521627D410>,
        ('helicopter', 'time_mobile'): <utils.SeededDistribution object at 0x0000015216267E90>}
        """
        n = len(self.activity_time_distr)
        rngs = [default_rng(s) for s in seed_seq.spawn(n)]

        dist_map = {
            "poisson": poisson,
            "bernoulli": bernoulli,
            "triang": triang,
            "erlang": erlang,
            "weibull_min": weibull_min,
            "expon_weib": exponweib,
            "betabinom": betabinom,
            "pearson3": pearson3,
            "cauchy": cauchy,
            "chi2": chi2,
            "expon": expon,
            "exponential": expon,  # alias for consistency
            "exponpow": exponpow,
            "gamma": gamma,
            "lognorm": lognorm,
            "norm": norm,
            "normal": norm,  # alias for consistency
            "powerlaw": powerlaw,
            "rayleigh": rayleigh,
            "uniform": uniform,
        }

        seeded_distributions = {}

        # print(activity_time_distr)
        # print(len(self.activity_time_distr))
        i = 0
        for entry, rng in zip(self.activity_time_distr, rngs):
            i += 1
            vt = entry["vehicle_type"]
            tt = entry["time_type"]
            best_fit = entry["best_fit"]

            # dist_name = best_fit['dist'].lower()
            dist_name = list(best_fit.keys())[0]
            dist_cls = dist_map.get(dist_name)
            self.debug(
                f"{i}/{len(self.activity_time_distr)} for {vt} {tt} set up {dist_name} {dist_cls} with params: {best_fit[dist_name]}"
            )

            if dist_cls is None:
                raise ValueError(f"Unsupported distribution type: {dist_name}")

            params = best_fit[dist_name]

            seeded_distributions[(vt, tt)] = SeededDistribution(dist_cls, rng, **params)

        self.activity_time_distr = seeded_distributions

        return seeded_distributions

Each fitted entry is looked up in a name-to-scipy-class map (dist_map). Unknown names raise ValueError; exponential and normal are aliases. The generator for entry i is the i-th spawned child, so the order of entries in activity_time_distributions.txt decides which entry gets which stream.

SeededDistribution is a thin pair of a frozen scipy distribution and its generator:

air_ambulance_des/utils.py:1053-1059

class SeededDistribution:
    def __init__(self, dist, rng, **kwargs):
        self.dist = dist(**kwargs)
        self.rng = rng

    def sample(self, size=None):
        return self.dist.rvs(size=size, random_state=self.rng)

activity_time then draws from it, rejecting values outside the bounds for that time_type:

air_ambulance_des/utils.py:580-607

    def activity_time(self, vehicle_type: str, time_type: str) -> float:
        """
        This function will return a dictionary containing
        the distribution and parameters for the distribution
        that match the provided HEMS vehicle type and time type

        """

        dist = self.activity_time_distr.get((vehicle_type, time_type))
        # print(dist)

        if dist is None:
            raise ValueError(f"No distribution found for ({vehicle_type}, {time_type})")

        try:
            min_time, max_time = self.min_max_cache[time_type]
            # print(f"Min time {time_type}: {min_time}")
            # print(f"Max time {time_type}: {max_time}")
        except KeyError:
            raise ValueError(f"Min/max bounds not found for time_type='{time_type}'")

        sampled_time = -1
        while not (min_time <= sampled_time <= max_time):
            sampled_time = dist.sample()

        # print(sampled_time)

        return sampled_time

The bounds come from min_max_cache, built in the constructor from upper_allowable_time_bounds.csv:

air_ambulance_des/utils.py:177-182

        # Turn the min-max activity times data into a format that supports easier/faster
        # lookups
        self.min_max_cache = {
            row["time"]: (row["min_value_mins"], row["max_value_mins"])
            for _, row in self.min_max_values_df.iterrows()
        }

Note the lookup is by time_type only. The bounds are shared across vehicle types.

Ages and daily call counts use the other path, sample_from_distribution, which takes a one-entry dict such as {"gamma": {...}} and the stream to use:

air_ambulance_des/utils.py:673-698

    def sample_from_distribution(self, distr: dict, rng: np.random.Generator) -> float:
        """
        Sample a single float value from a seeded scipy distribution.

        Parameters
        ----------
        distr : dict
            A dictionary with one key (the distribution name) and value as the parameters.
        rng : np.random.Generator
            A seeded RNG from the simulation's RNG stream pool.

        Returns
        -------
        float
            A positive sampled value from the specified distribution.
        """
        if len(distr) != 1:
            raise ValueError("Expected one distribution name in distr dictionary")

        dist_name, params = list(distr.items())[0]
        sci_distr = getattr(scipy.stats, dist_name)

        while True:
            sampled_value = np.floor(sci_distr.rvs(random_state=rng, **params))
            if sampled_value > 0:
                return sampled_value

It floors the draw and loops until the result is positive. age_sampling and inc_per_day wrap it in further loops (a maximum age; per-season or per-quarter min and max counts).

The model itself calls inc_per_day_samples, not inc_per_day. It resamples from empirical daily counts stored in inc_per_day_samples.json:

air_ambulance_des/utils.py:652-671

    def inc_per_day_samples(
        self, quarter: int, quarter_or_season: str = "season"
    ) -> float:
        """
        This function will return a dictionary containing
        the distribution and parameters for the distribution
        that match the provided HEMS vehicle type and time type

        """

        season = "summer" if quarter in [2, 3] else "winter"

        if quarter_or_season == "season":
            return self.rngs["calls_per_day"].choice(
                self.incident_per_day_samples[season]
            )
        else:
            return self.rngs["calls_per_day"].choice(
                self.incident_per_day_samples[f"Q{quarter}"]
            )

quarter_or_season picks the JSON key ("summer"/"winter" or "Q1" to "Q4"). DES_HEMS sets daily_calls_by_quarter_or_season = "quarter" (des_hems.py:76), so the quarter keys are the live path; the season default here is not what the model uses.

Ad hoc unavailability

sample_ad_hoc_reason(hour, quarter, registration) maps the hour to a six-hour bin, filters ad_hoc_unavailability.csv by registration, bin and quarter, and falls back in stages:

1. no rows for that registration: use all registrations for that bin and quarter;
2. still none: return "available";
3. any NaN probability: impute the registration’s mean across bins for that reason, else the bin’s mean across registrations, else return "available";
4. probabilities still missing or summing to zero: return "available".

Only then does it normalise and draw:

air_ambulance_des/utils.py:1044-1050

        norm_probs = filled_subset["probability"] / total_prob

        sampled_reason = self.rngs["ad_hoc_reason_selection"].choice(
            filled_subset["reason"].tolist(), p=norm_probs.tolist()
        )

        return sampled_reason

Every degraded path returns "available", so missing data silently makes aircraft more available. The constructor is consistent with this: if ad_hoc_unavailability.csv is absent, the FileNotFoundError is caught (utils.py:131-138) and self.ad_hoc_probs is simply never set. The failure then surfaces as an AttributeError at the first call, not at construction.

Extending this

Add a new sampled quantity (a new probability table).

1. Produce the CSV in DistributionFitUtils (Distribution Fitting), writing to distribution_data/.
2. Load it in Utils.__init__ next to the others.
3. Add a key to module_keys in setup_seeds at the end of the list. See the first sharp edge below for why.
4. Write the selection method following ampds_code_selection: filter, then pd.Series.sample(..., random_state=self.rngs["<key>"]).
5. Call it from DES_HEMS (Calls and Patient Journeys).

Add a new activity time or vehicle type. Add the fit in the activity-time distributions file via the fitting pipeline, and add a row for the time_type in upper_allowable_time_bounds.csv. activity_time raises ValueError for a missing distribution or missing bounds. A new scipy family must be added to dist_map in build_seeded_distributions.

Change which file backs a sampler (for example per-season versus per-quarter call counts): the switch is daily_calls_by_quarter_or_season on DES_HEMS. The JSON must have the matching keys.

Conventions the code expects

• One stream key per sampling purpose, and the method reads self.rngs[<key>]. Keys are strings, so a typo is a KeyError at the first sample.
• Probability tables have a proportion column. Method-specific filter columns are named in the existing methods.
• Tests that exercise sampling construct Utils(master_seed=SeedSequence(42)) and call setup_seeds() explicitly (tests/test_resource_allocation.py:34-35).

Tests that cover this area: tests/test_reproducibility.py, tests/test_against_reality.py, tests/test_backtest.py and tests/test_resource_allocation.py. The repo’s tests/README.md documents running them from the repo root:

Example usage (see tests/README.md)

pytest tests/test_reproducibility.py
pytest tests/test_resource_allocation.py
pytest -m backtest

If a change to sampling is meant to alter results, tests/README.md describes pytest -m backtest --update-backtest-baseline for refreshing the stored baseline after you have checked the differences.

Sharp edges

• Inserting a stream key shifts everyone after it. Children are numbered by position in spawn(len(module_keys)), and the activity-time generators are spawned afterwards from the same sequence. Adding a key anywhere moves the activity-time streams, and adding one mid-list also moves the later named streams. Results change even for a fixed seed. This follows from how SeedSequence.spawn works; the repo does not say it in a comment (inferred).
• setup_seeds() is not idempotent. build_seeded_distributions overwrites self.activity_time_distr (a list of dicts) with a dict. A second call would iterate the dict’s tuple keys and fail on entry["vehicle_type"].
• The default master_seed=SeedSequence(42) is evaluated once, at function definition. Every Utils() that uses the default shares one SeedSequence object and its spawn counter. Two such objects do not get identical streams.
• Bare Utils() objects have no streams. HEMS objects (Chapter 5) each build a bare Utils() for their time-range helper (class_hems.py:33). It has no rngs because setup_seeds() is never called, so never sample from it.
• Import-time file reads. HEMS_ROTA, HEMS_ROTA_DEFAULT and SERVICING_SCHEDULES_BY_MODEL are read when the module is imported (utils.py:49-56), not when Utils is constructed. The app overwrites those CSVs (Streamlit App and Scenario Parameters), so a process that imported utils earlier holds the old class-level frames. Utils.__init__ itself reads distribution_data/ fresh each time.
• Rejection loops have no cap. activity_time, sample_from_distribution, age_sampling and inc_per_day loop until a draw lands in range. Bounds that the fitted distribution almost never reaches will hang a run rather than raise.
• Stale-looking TODO: RANDOM SEED SETTING comments sit above sex_selection and age_sampling (utils.py:535, :552), yet both already use named streams. Whatever they originally flagged, the current code is seeded.
• Ad hoc fallbacks are silent. See above: every degraded path returns "available", and sample_ad_hoc_reason is annotated -> bool while returning a string.
• Hour ranges are inclusive in one place. hems_results_by_patient_outcome_and_time_of_day_and_quarter_and_vehicle_type_and_callsign_group_probs treats hours 7 to 18 inclusive as "day", so a call at 18:xx is day; sample_ad_hoc_reason uses different, six-hour bins.
• Two methods lack self. Utils.current_time() and Utils.biased_mean(series, bias) are defined without self, so they only work when called on the class, not an instance. Commit 7f5798f notes biased_mean “may need watching”.

Design notes

Documented.

• Random number streams for activity times were introduced in commit 2af661f, and commit 9a05b28 moved seed handling to numpy seed sequences “more robustly”. Together these are the origin of the SeedSequence plus named-generator design.
• Commit 2eb0ffd fixed different results in tests where a seed was set, which is the behaviour tests/test_reproducibility.py now guards.
• Commit 7f5798f added the option to stratify calls per day by season or quarter, which is the quarter_or_season argument.
• Commit ceb157e switched vehicle type sampling to per quarter, which is why vehicle_type_selection_qtr is live and the non-quarter version is commented out.
• Commit 5ad2dec moved utils.py into the air_ambulance_des package; file history before that move was not followed.

Inferred.

• The named-stream split presumably exists so that changing one part of the model (say, adding a sampled field) does not reshuffle unrelated draws, which makes scenario comparisons less noisy. The code is consistent with that but no commit message says so.
• Pairing each scipy distribution with its own generator in SeededDistribution suggests the aim of keeping every (vehicle_type, time_type) fit independent of how many times the others are drawn.
• The large block of commented-out methods and # ARCHIVED CODE suggests the sampling tables were reworked more than once; the live versions are the quarter-stratified ones. This is a reading of the code only.

Limitations and open questions

NoteWhy the model makes simplifying assumptions

Every simulation model rests on assumptions and simplifications, and it has to. No model can be a perfect copy of reality. Adding ever more detail does not fix that, and can make things worse: each extra detail is one more thing to get wrong, so an over-complicated model can drift further from reality, not closer. Judge the items below, and any change you are tempted to make, against whether the assumption is reasonable for the questions asked of the model. Closing every gap with reality is not the aim, and extra fidelity can add more error than it removes.

• Issue #183 (open as of 2026-10-19): add tests to quantify how much the random seed matters.
• Issue #241 (open as of 2026-10-19): a car is sometimes sent when the ad hoc unavailability spreadsheet gives no obvious reason, which touches the sampling fallbacks above.
• The two TODO: RANDOM SEED SETTING markers at air_ambulance_des/utils.py:535 and :552.

Check the issue tracker for the current state.

Conclusion

With sampling and seeding in hand, Calls and Patient Journeys shows where DES_HEMS calls these methods as each call becomes a patient.