flowchart LR
A[Cleaned job extract CSV] --> D[DistributionFitUtils]
B[actual_data/upper_allowable_time_bounds.csv] --> D
D --> E[distribution_data/]
D --> F[historical_data/]
D --> G[historical_data/calculated/]
E --> H[Utils: sampling in a run]
F --> I[HistoricResults]
G --> I
D -.runs the model.-> J[parallelProcessJoblib]
J -.-> G
Chapter 8: Distribution Fitting
How DistributionFitUtils turns a historical job extract into the distribution and historical-data files the rest of the system reads.
In Trial Results and Historic Comparison we saw HistoricResults loading aggregated files from historical_data/, and in Sampling Distributions and Random Streams we saw Utils loading fitted inputs from distribution_data/. Both chapters treated those files as a given. This chapter is about where they come from.
Why fit distributions at all?
The simulation does not replay history. When it needs an activity time, a patient age or the number of calls on a given day, it draws a random value from something that describes the real data. Those descriptions have to be built from somewhere.
Our central use case is this: a new extract of cleaned job data arrives (say, another year of Devon Air Ambulance jobs), and you want the model and the results pages to reflect it. Doing that by hand would mean working out dozens of probability tables, picking a suitable statistical distribution for each phase of a job, and aggregating the observed totals that HistoricResults compares against.
DistributionFitUtils automates the lot. It is an offline pipeline: you run it occasionally from the command line, not during a simulation run and not from the Streamlit app. Its output is a set of files that the rest of the system reads.
Distribution fitting pipeline: the offline sequence of steps that turns a job-level CSV into the files in distribution_data/ and historical_data/. See the glossary.
Inputs and outputs
Figure 1 shows the whole picture. One CSV goes in, three main families of output come out, and each is consumed by a part of the system you have already met.
distribution_data/holds probability tables, fitted distributions and daily call-count samples.Utilsloads these when a model is built.historical_data/holds aggregated observed totals, such as jobs per month or missed calls by hour.HistoricResultsreads these.historical_data/calculated/holds the output of a simulated historical baseline: the model run with the historic rotas and servicing, so that you can compare like with like.actual_data/is also written:school_holidays.csv(if requested), and the rota, lookup and service files are overwritten by the baseline step, so back up any scenario you care about first.
Simulated historical baseline: a set of model runs using historic rotas and servicing rather than the user’s scenario, giving the simulation’s own picture of “how things were”. See the glossary.
Running the pipeline
The entry point is the small script run_dist_fit_utils.py at the repository root.
run_dist_fit_utils.py:1-10
if __name__ == "__main__":
from air_ambulance_des.distribution_fit_utils import DistributionFitUtils
test = DistributionFitUtils(
"external_data/clean_daa_import_missing_2023_2024.csv",
calculate_school_holidays=True,
)
# test = DistributionFitUtils('external_data/clean_daa_import.csv')
test.import_and_wrangle()
test.run_sim_on_historical_params()The constructor takes the path to the extract. import_and_wrangle() does the fitting and writing; run_sim_on_historical_params() then produces the baseline. Passing calculate_school_holidays=True also regenerates the school holiday dates used by the servicing logic in HEMS Resources and Servicing.
Two practical points follow from how the class is written:
- File paths such as
"distribution_data/hour_by_ampds_card_probs.csv"are relative strings, so the script has to be run from the repository root. distribution_fit_utils.pyuses bare imports (from utils import Utils,from des_parallel_process import ...), as some other modules in the package do. We describe them as they are; just be aware that the module expects the package directory to be importable in that way.
Key concept: a best-fit distribution
Most of the continuous inputs (activity times, ages, daily call counts) are stored as a best-fit distribution: the name and parameters of whichever candidate from a fixed list matched the data most closely. The choice is delegated to the Fitter library, wrapped by getBestFit.
air_ambulance_des/distribution_fit_utils.py:142-163
def getBestFit(self, q_times, distr=get_common_distributions(), show_summary=False):
"""
Convenience function for Fitter.
Returns model and parameters that is considered
the 'best fit'.
TODO: Determine how Fitter works this out
"""
if q_times.size > 0:
if len(distr) > 0:
f = Fitter(q_times, timeout=60, distributions=distr)
else:
f = Fitter(q_times, timeout=60)
f.fit()
if show_summary == True:
f.summary()
return f.get_best()
else:
return {}Fitter.fit() tries each candidate distribution against the sample and get_best() returns a dictionary such as {"lognorm": {...parameters...}}. If there is no data at all (q_times.size is zero) it returns an empty dictionary instead of failing, which matters because some combinations of vehicle type and HEMS result have no jobs.
Which candidates are tried? The constructor hard-codes a list, self.sim_tools_distr_plus (poisson, gamma, lognorm, weibull_min and so on). A comment in the constructor notes that this was fixed deliberately, rather than using get_common_distributions(), to make random number generation setup more robust. That fits the sampling side: Sampling Distributions and Random Streams builds samplers from these names, so every name written here has to be one the sampler understands.
What import_and_wrangle does
import_and_wrangle is the conductor. It first reads the CSV and adds helper columns (date_only, hour, day_of_week, month, quarter, first_day_of_month) derived from inc_date. It also reads actual_data/upper_allowable_time_bounds.csv, which holds the minimum and maximum plausible value for each activity time. These bounds replaced an earlier approach based on the upper quartile, and have the advantage that they can be edited from the front end.
Then it clears out the old results and calls one method per output.
air_ambulance_des/distribution_fit_utils.py:212-252
self.removeExistingResults(Utils.HISTORICAL_FOLDER)
self.removeExistingResults(Utils.DISTRIBUTION_FOLDER)
# get proportions of AMPDS card by hour of day
self.hour_by_ampds_card_probs()
# Determine 'best' distributions for time-based data
self.activity_time_distributions()
# Calculate probability patient will be female based on AMPDS card
self.sex_by_ampds_card_probs()
# Determine 'best' distributions for age ranges straitifed by AMPDS card
self.age_distributions()
# Alternative approach to IA times. Start with probabilty of call at given hour stratified by quarter
self.hourly_arrival_by_qtr_probs()
# Calculates the mean and standard deviation of the number of incidents per day stratified by quarter
self.incidents_per_day()
self.incidents_per_day_samples()
# Calculate probability of enhanced or critical care being required based on AMPDS card
self.enhanced_or_critical_care_by_ampds_card_probs()
# Calculate HEMS result
self.hems_results_by_patient_outcome_and_time_of_day_and_quarter_and_vehicle_type_and_callsign_group_probs()
# Calculate probability of callsign being allocated to a job based on AMPDS card and hour of day
# self.callsign_group_by_ampds_card_and_hour_probs()
# self.callsign_group_by_ampds_card_probs()
# self.callsign_group_by_care_category()
self.callsign_group()
# Calculate probability of a particular vehicle type based on callsign group and month of year
# self.vehicle_type_by_month_probs()
self.vehicle_type_by_quarter_probs()
# self.vehicle_type_probs() # Similar to previous but without monthly stratification since ad hoc unavailability should account for this.
# Calculate the patient outcome (conveyed, deceased, unknown)
self.patient_outcome_by_care_category_and_quarter_probs()The two removeExistingResults calls delete every file with an extension at the top level of historical_data/ and distribution_data/ (subfolders, including historical_data/calculated/, are left alone), so an interrupted run leaves you with missing files rather than a mixture of old and new. After that, each call writes its own file. Some alternatives are commented out, and the original source labels others as archived code; we do not cover those.
Below this excerpt, the method calls the school-holiday calculation (if requested), then the twelve historical_* methods that write to historical_data/, and finally ad_hoc_unavailability(), wrapped in a try that prints a message if its input file is missing.
Probability tables
The simplest outputs are tables of proportions. hourly_arrival_by_qtr_probs is typical: count calls per hour within each quarter, divide by the quarter’s total, and write a CSV.
air_ambulance_des/distribution_fit_utils.py:863-873
ia_df = self.df[["quarter", "hour"]].dropna()
hourly_counts = (
ia_df.groupby(["hour", "quarter"]).size().reset_index(name="count")
)
total_counts = hourly_counts.groupby(["quarter"])["count"].transform("sum")
hourly_counts["proportion"] = round(hourly_counts["count"] / total_counts, 4)
hourly_counts.sort_values(by=["quarter", "hour"]).to_csv(
"distribution_data/hourly_arrival_by_qtr_probs.csv", mode="w+"
)This is the hourly arrival profile that Calls and Patient Journeys uses in calls_per_hour() to spread a day’s calls over the hours. hour_by_ampds_card_probs() follows the same shape, grouping by hour and ampds_card, and many of the other *_probs methods do likewise with different grouping columns (care category, vehicle type, callsign group, quarter).
Fitting activity times
activity_time_distributions is where getBestFit earns its keep. Each HEMS result has its own list of phases that apply (a “Stand Down” has no time on scene, for instance), held in self.times_to_fit. The method loops over every HEMS result, every phase in its list and every vehicle type.
air_ambulance_des/distribution_fit_utils.py:370-421
for row in self.times_to_fit:
# print(row)
for ttf in row["times_to_fit"]:
for vt in vehicle_type:
# print(f"HEMS result is {row['hems_result']} times_to_fit is {ttf} and vehicle type is {vt}")
# This line might not be required if data quality is determined when importing the data
max_time = self.min_max_values_df[
self.min_max_values_df["time"] == ttf
].max_value_mins.iloc[0]
min_time = self.min_max_values_df[
self.min_max_values_df["time"] == ttf
].min_value_mins.iloc[0]
# print(f"Max time is {max_time} and Min time is {min_time}")
if ttf == "time_on_scene":
# There is virtually no data for HEMS_result other than patient conveyed
# which is causing issues with fitting. For time on scene, will
# use a simplified fitting ignoring hems_result as a category
fit_times = self.df[
(self.df.vehicle_type == vt)
& (self.df[ttf] >= min_time)
& (self.df[ttf] <= max_time)
][ttf]
else:
fit_times = self.df[
(self.df.vehicle_type == vt)
& (self.df[ttf] >= min_time)
& (self.df[ttf] <= max_time)
& (self.df.hems_result == row["hems_result"])
][ttf]
# print(fit_times[:10])
best_fit = self.getBestFit(
fit_times, distr=self.sim_tools_distr_plus
)
# print(best_fit)
return_dict = {
"vehicle_type": vt,
"time_type": ttf,
"best_fit": best_fit,
"hems_result": row["hems_result"],
"n": len(fit_times),
}
# print(return_dict)
final_distr.append(return_dict)
with open(
"distribution_data/activity_time_distributions.txt", "w+"
) as convert_file:
convert_file.write(json.dumps(final_distr))For each combination, it selects the matching times, drops anything outside the bounds read earlier, fits, and records the winner along with n, the number of observations used. The special case for time_on_scene ignores the HEMS result, because the comment in the code says there is very little data for results other than “patient conveyed”, which caused fitting problems. The list of dictionaries is written as JSON to distribution_data/activity_time_distributions.txt, which is the file Utils wraps in SeededDistribution objects. Note that Utils keys the loaded fits only by (vehicle type, time type), so when several HEMS results give an entry for the same pair, only the last one in the file is used when sampling (see Sampling Distributions and Random Streams).
age_distributions() works in the same spirit, with one fit per AMPDS card.
Daily call counts
How many calls arrive each day is handled twice. incidents_per_day() fits a distribution by season and by quarter, giving inc_per_day_distributions.txt and inc_per_day_qtr_distributions.txt. Alongside that, incidents_per_day_samples stores the empirical counts directly.
air_ambulance_des/distribution_fit_utils.py:612-645
# Year weight map
if weight_map is None:
weight_map = {2023: 1.0, 2024: 1.1}
# Compute weights
daily_counts["weight"] = daily_counts["year"].map(weight_map).fillna(1.0)
# Storage
empirical_samples = {}
# Season-based
for season in daily_counts["season"].dropna().unique():
filtered = daily_counts[daily_counts["season"] == season].copy()
repeated = filtered.loc[
filtered.index.repeat(
(filtered["weight"] * scale_factor).round().astype(int)
)
]["jobs_per_day"].tolist()
empirical_samples[season] = repeated
# Quarter-based
for quarter in daily_counts["quarter"].dropna().unique():
filtered = daily_counts[daily_counts["quarter"] == quarter].copy()
repeated = filtered.loc[
filtered.index.repeat(
(filtered["weight"] * scale_factor).round().astype(int)
)
]["jobs_per_day"].tolist()
empirical_samples[f"Q{int(quarter)}"] = repeated
with open("distribution_data/inc_per_day_samples.json", "w") as f:
json.dump(empirical_samples, f)Each observed daily count is repeated in proportion to a weight for its year (default {2023: 1.0, 2024: 1.1}, scaled by scale_factor), so that recent, busier years count for more. The result is a dictionary of plain lists keyed by "winter", "summer", "Q1" to "Q4", saved to inc_per_day_samples.json. Drawing from such a list at random reproduces the observed spread of busy and quiet days, with no assumption that it follows a textbook curve.
Note that incidents_per_day() uses a different year-weight map from the one above, so check both methods before changing how recent years are treated.
Historical aggregates
The historical_* methods do no fitting. They simply summarise what actually happened, in the shape that HistoricResults expects.
air_ambulance_des/distribution_fit_utils.py:1548-1572
monthly_df = self.df[
["inc_date", "first_day_of_month", "hems_result", "vehicle_type"]
].drop_duplicates(subset="inc_date", keep="first")
is_stand_down = monthly_df["hems_result"].str.contains("Stand Down")
monthly_df["stand_down_car"] = (
(monthly_df["vehicle_type"] == "car") & is_stand_down
).astype(int)
monthly_df["stand_down_helicopter"] = (
(monthly_df["vehicle_type"] == "helicopter") & is_stand_down
).astype(int)
monthly_totals_df = (
monthly_df.groupby("first_day_of_month")
.agg(
stand_down_car=("stand_down_car", "sum"),
stand_down_helicopter=("stand_down_helicopter", "sum"),
total_jobs=("vehicle_type", "size"),
)
.reset_index()
)
monthly_totals_df.rename(columns={"first_day_of_month": "month"}).to_csv(
"historical_data/historical_jobs_per_month.csv", mode="w+", index=False
)historical_monthly_totals drops duplicate inc_date rows first, since several resources can attend one job, then counts jobs and stand-downs per month and writes historical_jobs_per_month.csv. The other methods in the family cut the data by callsign, day of week, hour of day, activity durations, utilisation, missed jobs and care category. upper_allowable_time_bounds also lives in this part of the class, but the pipeline currently reads the bounds from actual_data/ instead of calling it (the call in import_and_wrangle is commented out).
The simulated historical baseline
Observed totals describe the past; to judge a scenario fairly you also want the model’s version of the past. run_sim_on_historical_params() provides it. First it overwrites the rota, callsign lookup and servicing files in actual_data/ with the historic versions from tests/rotas_historic/.
air_ambulance_des/distribution_fit_utils.py:2691-2730
# Ensure all rotas are using default values
rota = pd.read_csv("tests/rotas_historic/HISTORIC_HEMS_ROTA.csv")
rota.to_csv("actual_data/HEMS_ROTA.csv", index=False)
callsign_reg_lookup = pd.read_csv(
"tests/rotas_historic/HISTORIC_callsign_registration_lookup.csv"
)
callsign_reg_lookup.to_csv(
"actual_data/callsign_registration_lookup.csv", index=False
)
service_history = pd.read_csv(
"tests/rotas_historic/HISTORIC_service_history.csv"
)
service_history.to_csv("actual_data/service_history.csv", index=False)
service_sched = pd.read_csv(
"tests/rotas_historic/HISTORIC_service_schedules_by_model.csv"
)
service_sched.to_csv("actual_data/service_schedules_by_model.csv", index=False)
print("Generating simulation results...")
removeExistingResults()
total_runs = 30
sim_years = 2
sim_duration = 60 * 24 * 7 * 52 * sim_years
parallelProcessJoblib(
total_runs=total_runs,
sim_duration=sim_duration,
warm_up_time=0,
sim_start_date=datetime.strptime(
"2023-01-01 05:00:00", "%Y-%m-%d %H:%M:%S"
),
amb_data=False,
print_debug_messages=True,
)
collateRunResults()Then it runs 30 replications of 104 weeks (728 days), with no warm-up and debug output switched on, using parallelProcessJoblib and collateRunResults exactly as described in Simulation Runs and Replications, and the collated data/run_results.csv is summarised into the historical_data/calculated/SIM_hist_params_* files.
WarningKnown issue: the missed-jobs summary cannot currently run
The method ends by calling _job_outcome_calculation.get_missed_call_df(...) to write the two SIM_hist_params_missed_jobs_care_cat_* files. But the import that provides that name is commented out at the top of the module, under a TODO: Important fix for rerunning with new input data:
air_ambulance_des/distribution_fit_utils.py:20-21
# TODO: Important fix for rerunning with new input data
# import visualisation._job_outcome_calculation as _job_outcome_calculationNo other line defines _job_outcome_calculation, so those calls would raise a NameError. The surrounding try only catches FileNotFoundError, so the error would propagate. The earlier steps (the two “suboptimal” summaries) have already been written by then. The SIM_hist_params_missed_jobs_* files currently in historical_data/calculated/ were therefore produced by an earlier version of the code. Treat this as something to fix before relying on a fresh run.
run_sim_on_historical_params is not part of the generated API reference, so it appears here as plain code only.
Under the hood: one pass through the pipeline
- You construct
DistributionFitUtilswith the extract path. Nothing is read yet; the constructor stores the path and the lists of times and candidate distributions. import_and_wrangle()loads the CSV, adds date and time columns, and reads the allowable time bounds.- It deletes every file with an extension at the top level of
historical_data/anddistribution_data/(subfolders, includinghistorical_data/calculated/, are left alone). - Each fitting method filters the DataFrame, then either counts and divides (probability tables) or calls
getBestFit(distributions), and writes one file. - The
historical_*methods write observed aggregates. run_sim_on_historical_params()swaps in the historic rota files and runs the model to build the baseline.
sequenceDiagram
participant U as run_dist_fit_utils.py
participant D as DistributionFitUtils
participant F as Fitter
participant FS as Output folders
participant P as parallelProcessJoblib
U->>D: import_and_wrangle()
D->>FS: removeExistingResults
D->>D: activity_time_distributions()
D->>F: getBestFit(times)
F-->>D: best fit and parameters
D->>FS: write distribution_data files
D->>FS: write historical_data files
U->>D: run_sim_on_historical_params()
D->>P: run 30 replications
P-->>D: run_results.csv
D->>FS: write calculated baseline files
The design is deliberately plain: one long class, one method per output file, and the folder layout acting as the contract with the rest of the system. The cost of that simplicity is that the pipeline does little checking. The source itself notes a TODO to validate the imported CSV for missing values and wrong column names, so a malformed extract usually fails deep inside a method.
Conclusion
DistributionFitUtils is the start of the data’s journey through this project. It reads a cleaned job extract, fits distributions with Fitter, tabulates probabilities and call-count samples for Utils, writes observed aggregates for HistoricResults, and runs a baseline simulation on historic rotas for comparison.
That completes the tour: from the Streamlit parameters in chapter 1, through seeded runs, patient journeys, resources and allocation, to the results and the fitting that underpins them. To look up a term, see the glossary. For signatures and parameters of every documented class and function, head to the API reference.