Simulation for what-if questions

“What happens to ED waits if we get a third of discharges out before noon?” A spreadsheet can’t answer that. Discharging at 10:00 instead of 16:00 doesn’t change admissions, midnight bed counts or length of stay in days, so every average in the spreadsheet stays the same. What changes is timing, and timing only shows up in a model that has a clock in it.

This article covers when a simulation is worth building, how to build a small one in plain Python and NumPy, and the three habits that make its answers trustworthy: a warm-up period, many replications, and validation against things you already know. The examples come from my bed occupancy simulation, which models 120 acute medical beds on synthetic data.

When a simulation beats a spreadsheet

A spreadsheet is the right tool when the question is about totals or averages. Beds needed on average = admissions per day × mean length of stay. That is Little’s law, and it holds whatever the distributions look like (queueing basics for capacity planning covers it).

Reach for a simulation when any of these is true:

• The question is “how often” or “how long”, not “how many”. Days with patients waiting for a bed, or 12-hour waits in the ED, depend on variation around the average, and a spreadsheet has no variation in it.
• Timing matters. Morning discharges, weekend discharge rates and clinic templates change when capacity frees up, not how much there is.
• Things interact. A patient waiting for a bed delays the next one. Queues feed back on themselves, and the effect grows fast near capacity.
• The textbook formulas don’t fit. Erlang C assumes a constant arrival rate and exponential service times. Admissions follow a daily and weekly cycle, and discharges bunch in the afternoon.

If none of these apply, don’t build a simulation. It costs more to build, explain and maintain than a spreadsheet, and a wrong simulation looks just as convincing as a right one.

Building one in plain Python and NumPy

You don’t need a simulation package for most operational questions. The bed model has three parts, inputs, an engine and measurement, and the first two come to under 300 lines of Python.

Inputs as named constants. Every planted pattern (arrivals by hour, the Monday peak, length of stay by patient group, discharge times, the weekend effect) sits in one file, assumptions.py. When a manager asks “what did you assume about weekends?”, the answer is one line of code.

Arrivals, vectorised. Emergency admissions arrive at random, at a rate that varies by hour and weekday. Draw a Poisson count for every hour of the run, then spread each hour’s arrivals uniformly within it:

counts = arrivals_rng.poisson(rates)
arrival = np.repeat(np.arange(len(rates)), counts) + arrivals_rng.random(counts.sum())
arrival = np.sort(arrival)

The engine. A first-come-first-served queue with many identical servers needs no event calendar. Process arrivals in order and keep a heap of the times at which each bed next becomes free. Each patient takes the earliest free bed, or waits for it:

def first_come_first_served(
    arrival: Sequence[float], servers: int, leave: Callable[[int, float], float]
) -> tuple[np.ndarray, np.ndarray]:
    """A queue with identical servers (beds) that starts empty.

    Arrivals are served in order. Each takes the server that frees up first, waiting if none is
    free; `leave(i, start)` says when arrival i releases its server.
    """
    free_at = [0.0] * servers  # a heap of the times each server next becomes free
    starts, ends = [], []
    for i, arrived in enumerate(arrival):
        start = max(arrived, free_at[0])
        end = leave(i, start)
        heapq.heapreplace(free_at, end)
        starts.append(start)
        ends.append(end)
    return np.array(starts), np.array(ends)

leave holds the domain logic: the patient’s ready day, a delay if that falls at a weekend, and a discharge time of day. Keeping it separate from the queue means the same engine can run a textbook check (below) with a one-line leave. This loop handles about a million patients in a few seconds. The whole project, 1,200 simulated years plus the checks below, runs in under two minutes on a laptop.

If your problem has priorities, patients who give up and leave, or capacity that changes through the day, you need a proper event list or a library such as SimPy. Start with the simplest engine that answers the question.

Measurement. Record each patient’s arrival, bed and discharge times, then count occupancy at every hour with a difference array: +1 when someone enters, −1 when they leave, and a cumulative sum.

def counts_at_each_hour(enter: np.ndarray, exit_: np.ndarray, hours: int) -> np.ndarray:
    """How many people are between `enter` and `exit_` at each whole hour (exit not included)."""

    def by_hour(times: np.ndarray) -> np.ndarray:
        index = np.minimum(np.ceil(times).astype(int), hours)
        return np.bincount(index, minlength=hours + 1)[:hours]

    return np.cumsum(by_hour(enter) - by_hour(exit_))

Warm-up

A simulation that starts with an empty ward spends its first weeks filling up, and those weeks are not what you want to measure. Either start from a realistic state or throw the start away. Throwing it away is simpler and harder to get wrong.

To choose how much to discard, run many replications from empty, average them day by day, and look for the point where the average stops climbing. This is the idea behind Welch’s graphical method.

Data table

Day of the run	Mean of replications	Mean of replications (low)	Mean of replications (high)
1	8%	5%	10%
2	24%	19%	29%
3	35%	29%	42%
4	44%	36%	51%
5	50%	42%	57%
6	56%	49%	64%
7	64%	56%	72%
8	68%	60%	78%
9	67%	59%	76%
10	70%	61%	78%
11	71%	63%	81%
12	72%	63%	82%
13	76%	68%	85%
14	82%	73%	91%
15	84%	74%	93%
16	80%	69%	90%
17	81%	70%	92%
18	82%	71%	93%
19	82%	71%	94%
20	84%	74%	96%
21	90%	79%	100%
22	90%	80%	98%
23	84%	74%	94%
24	84%	73%	94%
25	84%	73%	95%
26	84%	75%	95%
27	86%	77%	97%
28	91%	82%	100%
29	91%	83%	99%
30	87%	76%	98%
31	87%	77%	97%
32	87%	77%	97%
33	87%	77%	97%
34	88%	78%	98%
35	93%	82%	100%
36	92%	83%	99%
37	87%	77%	98%
38	87%	77%	97%
39	87%	77%	97%
40	87%	78%	98%
41	90%	80%	99%
42	94%	85%	100%
43	93%	84%	100%
44	89%	79%	98%
45	89%	79%	99%
46	89%	80%	98%
47	88%	79%	98%
48	90%	80%	100%
49	94%	86%	100%
50	94%	85%	100%
51	89%	79%	99%
52	89%	78%	100%
53	89%	78%	99%
54	88%	77%	99%
55	90%	79%	99%
56	94%	84%	100%
57	94%	86%	100%
58	89%	79%	99%
59	89%	80%	99%
60	88%	78%	99%
61	88%	78%	99%
62	90%	80%	100%
63	94%	85%	100%
64	94%	86%	100%
65	90%	80%	100%
66	90%	79%	99%
67	89%	78%	99%
68	89%	78%	99%
69	90%	79%	100%
70	94%	85%	100%
71	94%	86%	100%
72	90%	80%	99%
73	89%	80%	99%
74	89%	80%	98%
75	89%	80%	99%
76	91%	81%	100%
77	94%	86%	100%
78	94%	86%	100%
79	89%	79%	100%
80	89%	79%	99%
81	88%	78%	98%
82	88%	78%	99%
83	90%	79%	100%
84	94%	86%	100%
85	94%	86%	100%
86	90%	81%	99%
87	90%	79%	100%
88	89%	78%	100%
89	89%	80%	99%
90	90%	81%	100%
91	94%	86%	100%
92	94%	88%	100%
93	90%	81%	100%
94	90%	80%	99%
95	90%	80%	99%
96	90%	80%	99%
97	91%	81%	100%
98	95%	86%	100%
99	94%	87%	100%
100	91%	80%	100%
101	90%	81%	99%
102	90%	82%	99%
103	90%	80%	99%
104	92%	82%	100%
105	95%	86%	100%
106	94%	86%	100%
107	90%	81%	100%
108	91%	81%	100%
109	90%	81%	99%
110	90%	80%	99%
111	91%	81%	100%
112	95%	87%	100%
113	94%	86%	100%
114	91%	82%	100%
115	91%	80%	100%
116	90%	80%	99%
117	90%	79%	99%
118	91%	81%	100%
119	94%	86%	100%
120	94%	85%	100%
121	90%	80%	100%
122	89%	79%	100%
123	89%	79%	99%
124	89%	78%	99%
125	90%	80%	100%
126	94%	85%	100%
127	94%	86%	100%
128	89%	79%	100%
129	89%	78%	99%
130	89%	78%	99%
131	89%	80%	99%
132	91%	80%	100%
133	95%	86%	100%
134	94%	86%	100%
135	90%	80%	99%
136	89%	80%	99%
137	89%	80%	98%
138	89%	79%	99%
139	91%	82%	100%
140	95%	86%	100%

Synthetic data. Source: projects/bed-occupancy-simulation. Current practice. Nothing before the end of warm-up is scored.

Mean bed occupancy was 40.0% in the first week, 90.3% in week 7 and 90.4% in week 8, against 90.9% over the scored weeks. After week 7, the weekly averages move within their ordinary noise. I discard eight weeks. Err on the long side: extra warm-up costs run time, too little biases every number.

Replications and confidence intervals

One run of a simulation is an anecdote. The same inputs with a different seed give a different year, and in this model the middle 80% of simulated years had between 88 and 156 days with patients waiting for a bed (the 10th and 90th percentiles). So run many independent replications, record one number per measure per replication, and summarise across them.

The estimate settles as replications are added:

Data table

Replications	Running mean	Running mean (low)	Running mean (high)
3	92	72	112
4	94	82	105
5	100	80	121
6	111	80	141
7	106	78	133
8	109	85	133
9	107	85	129
10	109	89	128
11	113	93	132
12	114	97	132
13	116	99	133
14	117	102	132
15	118	103	132
16	118	104	131
17	118	106	131
18	117	105	129
19	117	106	129
20	119	108	130
21	118	108	129
22	117	106	127
23	116	105	126
24	115	105	125
25	116	106	126
26	116	106	125
27	116	106	125
28	114	105	123
29	116	106	126
30	116	107	126
31	116	106	125
32	118	108	127
33	119	109	129
34	120	110	129
35	121	111	131
36	121	112	131
37	122	113	131
38	123	113	132
39	122	113	131
40	122	113	131
41	122	114	131
42	123	114	131
43	122	114	131
44	122	113	130
45	122	114	130
46	122	114	130
47	123	115	131
48	124	116	131
49	123	115	131
50	123	115	130
51	123	115	131
52	123	115	130
53	123	115	130
54	123	115	130
55	122	114	129
56	121	114	129
57	122	114	129
58	122	115	129
59	121	114	128
60	120	113	128
61	120	113	127
62	120	113	127
63	120	113	126
64	119	112	126
65	119	112	126
66	119	113	126
67	119	113	126
68	119	113	126
69	119	112	125
70	119	113	126
71	119	113	126
72	119	113	125
73	119	113	125
74	120	113	126
75	120	114	126
76	119	113	125
77	119	113	125
78	119	113	125
79	119	113	124
80	118	113	124
81	118	112	124
82	118	113	124
83	118	112	124
84	118	112	124
85	118	113	124
86	118	112	123
87	118	112	123
88	117	112	123
89	117	112	122
90	117	112	122
91	117	112	122
92	117	112	123
93	117	112	123
94	117	112	122
95	118	112	123
96	118	113	123
97	118	113	123
98	118	113	123
99	118	113	123
100	118	113	123
101	118	113	123
102	118	113	123
103	118	113	123
104	118	113	123
105	118	113	123
106	118	113	123
107	118	113	123
108	118	113	123
109	118	113	123
110	118	113	123
111	118	113	123
112	118	113	122
113	118	113	123
114	118	114	123
115	118	114	123
116	118	114	123
117	118	114	123
118	119	114	123
119	119	114	123
120	118	114	123
121	118	114	123
122	119	114	123
123	119	114	123
124	119	114	123
125	119	115	124
126	119	115	124
127	119	115	123
128	119	115	124
129	119	115	124
130	119	115	124
131	119	115	124
132	119	115	124
133	120	115	124
134	120	115	124
135	120	115	124
136	120	116	124
137	120	116	124
138	120	116	125
139	120	116	125
140	121	116	125
141	121	117	125
142	121	117	126
143	121	117	125
144	121	117	126
145	122	117	126
146	122	118	126
147	122	118	126
148	122	118	126
149	122	118	126
150	122	118	126
151	122	118	126
152	122	118	126
153	122	118	126
154	122	118	126
155	122	118	126
156	122	118	126
157	122	118	126
158	122	117	126
159	121	117	126
160	121	117	126
161	121	117	125
162	121	117	125
163	121	117	125
164	121	117	125
165	121	117	125
166	121	117	125
167	121	117	125
168	121	117	125
169	121	117	125
170	120	116	124
171	121	117	125
172	121	117	125
173	121	117	125
174	120	116	124
175	120	116	124
176	120	116	124
177	120	116	124
178	120	116	124
179	120	116	124
180	120	116	124
181	120	116	124
182	120	116	124
183	120	117	124
184	121	117	125
185	121	117	124
186	120	117	124
187	121	117	124
188	120	117	124
189	121	117	124
190	120	117	124
191	120	117	124
192	120	117	124
193	120	117	124
194	120	117	124
195	120	117	124
196	120	117	124
197	121	117	124
198	120	117	124
199	120	117	124
200	121	117	124

Synthetic data. Source: projects/bed-occupancy-simulation. Current practice. The interval is a t interval for the mean.

Five replications would have put the answer 20 days below the 200-replication estimate. The interval is a t interval, which matters with few replications, and statsmodels computes it in one line:

def interval(values: np.ndarray) -> tuple[float, float]:
    """95% t confidence interval for the mean across replications."""
    low, high = DescrStatsW(values).tconfint_mean()
    return float(low), float(high)

Keep two kinds of spread apart, because they answer different questions:

• The 10th to 90th percentile across replications says how different one year can look from the next. This is what the service will live through, and it doesn’t shrink with more replications.
• The confidence interval for the mean says how precisely the long-run average is known. It does shrink, in proportion to one over the square root of the number of replications: four times as many replications roughly halve it.

Common random numbers

When you compare scenarios, give them the same random numbers. In the bed model, replication k has the same patients arriving at the same times with the same stays in every scenario, so the only difference between scenarios is the policy. Separate random streams for arrivals and for patient characteristics keep this true even when a scenario adds demand:

def rng(replication: int, stream: int) -> np.random.Generator:
    """Independent, reproducible random streams for each replication and purpose."""
    return np.random.default_rng([SEED, replication, stream])

Then compare scenarios replication by replication. For moving a third of discharges before noon, the 95% interval for the change in days with patients waiting was ±0.41 days when paired this way, against ±5.04 if the scenarios had been run independently. For 12-hour waits, it was ±2.0 against ±40.6. A small but real effect that would drown in noise between independent runs becomes clear.

Validating against known behaviour

A simulation can be wrong in ways that look plausible. I check four things before believing any scenario.

Little’s law. The mean number of patients in a bed should equal admissions per hour times mean hours in a bed. In the bed model, it was 109.03 against 109.06. A mismatch means patients are being lost or counted twice.

A textbook special case. Strip the model down to something with a known answer and compare. With constant-rate Poisson arrivals and exponential stays averaging five days, 120 beds is an M/M/c queue, and the Erlang C formula gives its mean wait exactly. The same engine, with a one-line leave, ran 20 replications of ten years at each utilisation, discarding the first year of each:

Data table

Utilisation	Erlang C formula	Simulation	Simulation (low)	Simulation (high)
85%	0.4	0.4	0.4	0.5
90%	1.8	1.7	1.5	1.9
93%	4.7	4.4	4.0	4.8
95%	9.4	9.4	8.1	10.6
97%	21.7	24.2	19.8	28.5

Synthetic data. Source: projects/bed-occupancy-simulation. Whiskers: 95% interval across replications.

The formula falls inside the simulation’s 95% interval at all five utilisations, from 0.36 hours at 85% to 21.66 hours at 97%. The engine is doing what a queue should.

The inputs come back out. Planted inputs should reappear in the outputs: 15.3% of discharges before noon in current practice and 33.3% in the before-noon scenario, weekend discharges at 53% of a weekday rate, and 0.99 of the weekday rate once the weekend delay is switched off.

Real history. With real inputs, the model has to reproduce last year’s occupancy by hour of week, and the frequency of ED waits for a bed, before any scenario gets shown to anyone. Synthetic data can’t do this step, which is why a synthetic model shows method, not answers.

Communicating results to managers

The output of a simulation is a table of distributions. The person who asked wants a decision. A few habits bridge that gap:

• Lead with the answer to the question asked. “Matching weekend discharges to weekdays would have cut 12-hour waits by about 60% in this model. Moving a third of discharges before noon, on its own, by about 8%.”
• Give a range a year can fall in, not a confidence interval. A manager lives through one year at a time, so the 10th to 90th percentile is the honest spread. Keep confidence intervals for the analyst’s own checks and for comparisons between scenarios.
• Round. “About 8%” is more honest than “8.44%” from a model with invented inputs.
• Name the assumption that drives each result. The weekend result assumes weekday discharge rates are achievable at weekends, which needs weekend staffing across several services. The before-noon result assumes stays don’t get shorter, which is probably pessimistic.
• Say what the model leaves out. No escalation beds, no winter, one pool of beds. List these before anyone asks.
• Let people explore. A chart with scenario buttons, like the typical-week chart in the bed occupancy simulation, answers the next three questions without another meeting.

The purpose is a better conversation about the decision, not a number to put in a board paper. If the manager ends up asking sharper questions about weekend pharmacy cover, the simulation did its job.

Reproduce

All the charts and numbers here come from the bed occupancy project:

uv run python projects/bed-occupancy-simulation/run.py