12. Simulation
Many real systems are too complex for the closed-form models of the previous chapters: arrival rates change during the day, service times are not exponential, and resources interact in ways that no formula captures. Simulation imitates such a system on a computer with random inputs, runs it many times and estimates its performance statistically, which makes it the most flexible tool of OR and, in industry, one of the most widely used. This chapter first shows how to generate random inputs and organize a simulation model (12.1), then how to turn its noisy outputs into trustworthy answers through input modeling, confidence intervals, warm-up handling, variance reduction, validation and optimization (12.2).
<a id="121-simulation-fundamentals"></a>
12.1 Simulation Fundamentals
This subchapter builds a simulator from the ground up. Monte Carlo sampling estimates expectations by averaging random scenarios. Those scenarios need uniform random numbers, and methods that turn uniforms into any distribution (inverse transform, acceptance-rejection). Discrete-event simulation then arranges the samples in time to imitate a dynamic system. The running example is the coffee bar of Section 9.3: 20 customers per hour and a mean service time of 2.5 minutes, so the exact M/M/1 answer $W_q = 12.5$ minutes is available to check the simulation.
Monte Carlo Simulation
Why it matters. Many OR quantities are expectations or probabilities that cannot be computed exactly: the probability that a network stays connected after random failures, the expected cost of an inventory policy, the probability that a project finishes on time (Section 10.2). Monte Carlo replaces the exact computation by an experiment: draw many random scenarios, evaluate each one, and average. Its accuracy does not depend on how many random inputs the model has. That is why it works where enumeration and numerical integration fail, and why it is the basis of every other method in this chapter.
Theory. The goal is $\theta = E[g(X)]$, where $X$ is a random vector (possibly with thousands of components) and $g$ is anything that can be evaluated: a cost, a simulated waiting time, or the 0/1 indicator of an event, in which case $\theta$ is a probability. Draw $n$ independent samples $X_1, \dots, X_n$ and average:
$$ \hat\theta_n = \frac{1}{n} \sum_{i=1}^{n} g(X_i). $$
Two theorems make this work:
- Law of large numbers. $\hat\theta_n \to \theta$ as $n \to \infty$, and $E[\hat\theta_n] = \theta$ for every $n$ (the estimator is unbiased).
- Central limit theorem. For large $n$, $\hat\theta_n$ is approximately normal with standard deviation $\sigma/\sqrt{n}$, where $\sigma^2 = \mathrm{Var}(g(X))$. This gives the confidence interval
$$ \hat\theta_n \pm z_{1-\alpha/2} , \frac{S_n}{\sqrt{n}}, $$
with $S_n$ the sample standard deviation and $z_{0.975} = 1.96$ for 95% confidence. For a probability $p$, $\sigma^2 = p(1-p)$.
The standard error $\sigma/\sqrt{n}$ explains both the strength and the weakness of the method:
- Dimension-free. The error is $O(1/\sqrt{n})$ whatever the dimension of $X$. A grid with $m$ points per dimension needs $m^d$ evaluations, which is hopeless beyond a handful of dimensions.
- Slow. One more correct digit costs 100 times more samples. For a target half-width $\varepsilon$, the required sample size is $n \approx (z\sigma/\varepsilon)^2$.
- Weak on rare events. The relative error of an estimated probability is $\sqrt{(1-p)/(np)}$. Estimating $p = 10^{-6}$ with 10% relative error needs about $10^8$ samples. Importance sampling (Section 12.2) fixes this.
An analogy is measuring the area of a lake on a map by throwing stones at random over a rectangle of known area. The fraction of stones that land in the lake, times the area of the rectangle, estimates the lake's area. The shape of the shoreline does not matter, only the number of throws.
Monte Carlo is the static member of a family of simulation methods:
| Type | Time | Typical OR use |
|---|---|---|
| Monte Carlo (static) | no time, independent scenarios | risk analysis, reliability, project completion, valuation |
| Discrete-event simulation | state changes at event instants | queues, factories, hospitals, logistics |
| Continuous simulation (system dynamics) | differential equations, small time steps | supply chain dynamics, epidemics, policy models |
| Agent-based simulation | many autonomous interacting agents | markets, evacuation, adoption of new products |
Example: reliability of a road network. A distributor must know the probability that its warehouse can still be reached from the depot after a storm. Each of the 12 roads is blocked independently with probability 0.2. Two-terminal reliability is #P-hard in general, but this network has only $2^{12} = 4{,}096$ road states, so the exact value can be computed by enumeration and used to check the estimates.
1 ------------ 4
/ | / \
/ | / \
depot 0 --- 2 6 -------- 7 warehouse
\ \ / /
\ \/ /
3 -- 5 --------
import itertools
import numpy as np
# Road network after a storm: node 0 = depot, node 7 = warehouse.
# Each of the 12 roads is blocked independently with probability q.
edges = [(0, 1), (0, 2), (0, 3), (1, 2), (1, 4), (2, 5),
(3, 5), (4, 6), (5, 6), (4, 7), (6, 7), (5, 7)]
q = 0.2
def connected(up):
"""True if the warehouse (7) can be reached from the depot (0) on open roads."""
reach, stack = {0}, [0]
while stack:
u = stack.pop()
for k, (a, b) in enumerate(edges):
if up[k] and u in (a, b):
v = b if u == a else a
if v not in reach:
reach.add(v); stack.append(v)
return 7 in reach
# Exact answer by enumerating all 2^12 = 4,096 road states (impossible for 100+ roads)
exact = 0.0
for up in itertools.product([0, 1], repeat=len(edges)):
k = sum(up) # number of open roads
exact += connected(up) * (1 - q) ** k * q ** (len(edges) - k)
print(f"exact reliability = {exact:.5f}")
# Monte Carlo: sample road states, average the 0/1 indicator "connected"
rng = np.random.default_rng(2)
X = np.array([connected(rng.random(len(edges)) > q) for _ in range(100_000)])
print(" n estimate 95% half-width true error")
for n in (100, 1_000, 10_000, 100_000):
p_hat = X[:n].mean()
half = 1.96 * np.sqrt(p_hat * (1 - p_hat) / n) # CLT: standard error of a proportion
print(f"{n:8d} {p_hat:.5f} {half:.5f} {p_hat - exact:+.5f}")
# Replications needed for a half-width of 0.001
n_req = 1.96 ** 2 * exact * (1 - exact) / 0.001 ** 2
print(f"n for +/-0.001: {n_req:,.0f}")
# exact reliability = 0.94756
# n estimate 95% half-width true error
# 100 0.95000 0.04272 +0.00244
# 1000 0.95300 0.01312 +0.00544
# 10000 0.94920 0.00430 +0.00164
# 100000 0.94816 0.00137 +0.00060
# n for +/-0.001: 190,900
The warehouse stays reachable with probability 0.948. Every interval contains the exact value, and the half-width shrinks by about $\sqrt{10} \approx 3.2$ for every tenfold increase in samples. A precision of $\pm 0.001$ needs about 191,000 samples. Enumeration needed 4,096 states here, but a regional network with 100 roads has $2^{100} \approx 10^{30}$ states. Monte Carlo needs the same 191,000 samples for the same precision; only the cost of each connectivity check grows.
Key takeaways
- Monte Carlo estimates $E[g(X)]$ by the sample average of independent draws. It is unbiased and converges by the law of large numbers.
- The standard error is $\sigma/\sqrt{n}$: always report a confidence interval, and size $n$ from $n \approx (z\sigma/\varepsilon)^2$.
- The error rate does not depend on the dimension of $X$, but each extra digit costs 100 times more samples.
- Rare-event probabilities have huge relative errors under plain Monte Carlo and need variance reduction.
Practice
- (Beginner) Estimate $\pi$ by drawing $n$ random points in the unit square and counting those with $x^2 + y^2 \le 1$. Use $\sigma^2 = 16,p(1-p)$ with $p = \pi/4$ to compute how many points a half-width of 0.001 requires, then check your answer empirically.
- (Intermediate) The distributor can harden one road so that its blocking probability falls from 0.2 to 0.05. Using the same 100,000 sampled states for all 12 options (reuse each uniform and compare it with the new threshold), estimate the reliability gain of each option and pick the best road. Why does reusing the samples make the ranking more reliable?
- (Interview) "Our risk model ran 10,000 Monte Carlo scenarios and none produced a loss above 10 million, so that probability is zero." How do you respond? Hint: zero events in $n$ trials gives a 95% upper bound of about $3/n$ (the rule of three). Think about the relative error of rare-event estimates and about importance sampling.
Random Number Generation
Why it matters. Every simulation consumes uniform random numbers on $(0, 1)$, often billions of them, and transforms them into arrivals, service times and demands. Computers are deterministic, so they produce pseudo-random numbers: a deterministic sequence that passes statistical tests for randomness. Reproducibility is a feature: the same seed gives the same run, which makes debugging and fair comparisons possible. Bad generators have produced wrong published results, so the generator is not a detail to ignore.
Theory. A good generator is uniform and shows no detectable dependence (including between consecutive pairs and triples), has a very long period, is fast and reproducible, and supports many independent streams.
The classic design is the linear congruential generator (LCG):
$$ x_{k+1} = (a,x_k + c) \bmod m, \qquad u_k = x_k / m. $$
Its period is at most $m$. The Hull-Dobell theorem gives the conditions for the full period $m$ when $c \ne 0$: (i) $c$ and $m$ are coprime, (ii) $a - 1$ is divisible by every prime factor of $m$, (iii) $a - 1$ is divisible by 4 if $m$ is. Multiplicative generators ($c = 0$) with a prime modulus reach period $m - 1$ when $a$ is a primitive root, such as MINSTD with $a = 16807$ and $m = 2^{31} - 1$.
LCGs have a structural flaw (Marsaglia, 1968): consecutive $d$-tuples $(u_k, \dots, u_{k+d-1})$ lie on a small number of parallel hyperplanes. The infamous RANDU ($a = 65539 = 2^{16} + 3$, $c = 0$, $m = 2^{31}$) satisfies $a^2 \equiv 6a - 9 \pmod{2^{31}}$, so
$$ x_{k+2} = 6,x_{k+1} - 9,x_k \pmod{2^{31}}, $$
and every triple of consecutive numbers lies on one of only 15 planes. A three-dimensional simulation driven by RANDU explores a tiny part of the cube.
Modern generators fix these weaknesses:
| Generator | Period | Where it is used | Notes |
|---|---|---|---|
| LCG (MINSTD, RANDU) | $\le 2^{31}$ | legacy code | lattice structure, too short today |
| Mersenne Twister MT19937 | $2^{19937} - 1$ | Python random, NumPy RandomState, many libraries | fast, fails a few modern tests, large state |
| PCG64 | $2^{128}$ | NumPy default_rng | small state, excellent statistical quality |
| MRG32k3a (L'Ecuyer) | about $2^{191}$ | Arena, Simio, simulation libraries | built-in streams and substreams |
| Philox, Threefry (counter-based) | $2^{128}$ and more | GPUs, parallel computing | any position computable directly |
None of them is cryptographically secure. Use the secrets module for keys and tokens.
Generators are tested in two ways:
- Empirical tests compare samples with the uniform distribution: chi-square and Kolmogorov-Smirnov for uniformity, serial tests on pairs and triples, runs and gap tests. Large batteries such as TestU01 (BigCrush) and PractRand run hundreds of them.
- Theoretical tests analyze the generator's structure. The spectral test measures the distance between the hyperplanes of an LCG.
Passing tests never proves randomness, it only fails to detect patterns.
Seeds and streams. Use a separate stream for each source of randomness (arrivals, services, failures). Then a change in one part of the model does not shift the random numbers used elsewhere. This synchronization is what makes common random numbers effective (Section 12.2). For parallel replications, derive independent streams from one master seed (NumPy SeedSequence.spawn, jump-ahead, or counter-based generators) instead of using consecutive seeds with a weak generator. Always record the seeds.
Example: periods, RANDU's planes and independent streams.
import numpy as np
def lcg(a, c, m, seed, n):
"""Linear congruential generator: x_{k+1} = (a x_k + c) mod m, returns x_1..x_n."""
x, out = seed, np.empty(n, dtype=np.int64)
for k in range(n):
x = (a * x + c) % m
out[k] = x
return out
# Period of a toy generator with m = 16
print("a=5, c=3:", lcg(5, 3, 16, 7, 16)) # Hull-Dobell holds: full period 16
print("a=5, c=0:", lcg(5, 0, 16, 7, 16)) # c = 0 breaks it: period 4
# RANDU (IBM, 1960s): a = 65539, c = 0, m = 2^31
n = 300_000
u_randu = lcg(65539, 0, 2 ** 31, 1, n) / 2 ** 31
u_pcg = np.random.default_rng(1).random(n) # NumPy default generator (PCG64)
for name, u in (("RANDU", u_randu), ("PCG64", u_pcg)):
counts = np.histogram(u, bins=20, range=(0, 1))[0]
chi2 = ((counts - n / 20) ** 2 / (n / 20)).sum() # 1-D uniformity, 19 df, 5% critical value 30.14
r1 = np.corrcoef(u[:-1], u[1:])[0, 1] # lag-1 serial correlation
z = 9 * u[:-2] - 6 * u[1:-1] + u[2:] # RANDU: x_{k+2} = 6 x_{k+1} - 9 x_k (mod 2^31)
on = np.abs(z - np.round(z)) < 1e-9 # triple lies on a plane 9u - 6v + w = integer
print(f"{name}: chi2 = {chi2:4.1f} lag-1 corr = {r1:+.4f} "
f"triples on a plane: {on.mean():6.1%} planes: {len(np.unique(np.round(z[on])))}")
# Independent streams for parallel replications: spawn children of one master seed
streams = [np.random.default_rng(s) for s in np.random.SeedSequence(2026).spawn(3)]
print("first draw of each stream:", [round(g.random(), 4) for g in streams])
# a=5, c=3: [ 6 1 8 11 10 5 12 15 14 9 0 3 2 13 4 7]
# a=5, c=0: [ 3 15 11 7 3 15 11 7 3 15 11 7 3 15 11 7]
# RANDU: chi2 = 12.4 lag-1 corr = -0.0004 triples on a plane: 100.0% planes: 15
# PCG64: chi2 = 15.4 lag-1 corr = -0.0006 triples on a plane: 0.0% planes: 0
# first draw of each stream: [0.3827, 0.6711, 0.5508]
The toy generator with $c = 3$ visits all 16 values before repeating, as Hull-Dobell predicts. With $c = 0$ it cycles through only 4 values. RANDU passes the one-dimensional chi-square test (12.4, well below the critical value 30.14) and shows no serial correlation, so simple tests find nothing wrong. Yet 100% of its triples lie on exactly 15 planes, while PCG64 shows no such structure. A generator must be judged in the dimension in which the simulation uses it.
Key takeaways
- Simulations run on pseudo-random numbers: deterministic, reproducible sequences that pass statistical tests.
- LCGs are simple, but their tuples lie on hyperplanes. RANDU's triples lie on 15 planes. Use modern generators such as PCG64, MRG32k3a or Philox.
- One-dimensional tests are not enough. Check the generator in the dimension in which the simulation uses it.
- Use separate, recorded streams per source of randomness and per parallel replication. This gives reproducibility, independence and synchronization for common random numbers.
Practice
- (Beginner) Generate by hand the sequence $x_{k+1} = (5x_k + 1) \bmod 8$ from $x_0 = 0$ until it repeats. Is the period full? Check the three Hull-Dobell conditions.
- (Intermediate) Prove that $65539^2 \equiv 6 \cdot 65539 - 9 \pmod{2^{31}}$. Then run a two-dimensional serial test: split the unit square into a $10 \times 10$ grid, count the pairs $(u_{2k}, u_{2k+1})$ in each cell for RANDU and PCG64, and compute the chi-square statistic (99 degrees of freedom). Does the 2-D test detect RANDU's flaw? Why or why not?
- (Interview) "We run 100 replications on 100 cores and seed each with its core number, 0 to 99. Is that fine?" Hint: it depends on the generator and on how seeds are turned into states. Think about overlapping or correlated sequences, about spawning or jump-ahead, and about reproducing a single replication later.
Inverse Transform Method
Why it matters. A simulation needs samples from exponential, Weibull, empirical and discrete distributions, not just uniforms. The inverse transform method is the default way to get them. It is exact, needs exactly one uniform per sample, and maps uniforms to outputs monotonically. That last property keeps streams synchronized and is what makes common random numbers and antithetic variates work (Section 12.2).
Theory. Let $F$ be a continuous, strictly increasing cdf and $U \sim \text{Uniform}(0, 1)$. Then
$$ X = F^{-1}(U) \quad \text{has cdf } F, \qquad \text{because} \quad P(F^{-1}(U) \le x) = P(U \le F(x)) = F(x). $$
Graphically, the uniform is drawn on the vertical axis and read back through the cdf:
F(x)
1 | _________
| _____/
U | - - - - - - - - - - - -*
| __/ |
| ____/ |
| _____/ |
0 |________/_______________|______________ x
X = F^-1(U)
Steep parts of $F$ (high density) catch many uniforms in a narrow range of $x$, and flat parts spread them out. For a general cdf, including discrete ones, use the generalized inverse $F^{-1}(u) = \min{x : F(x) \ge u}$.
| Distribution | cdf $F(x)$ | Generator $X = F^{-1}(U)$ |
|---|---|---|
| Uniform$(a, b)$ | $(x - a)/(b - a)$ | $a + (b - a)U$ |
| Exponential, rate $\lambda$ | $1 - e^{-\lambda x}$ | $-\ln(1 - U)/\lambda$ |
| Weibull, shape $k$, scale $\eta$ | $1 - e^{-(x/\eta)^k}$ | $\eta,(-\ln(1 - U))^{1/k}$ |
| Pareto, shape $\alpha$, scale $x_m$ | $1 - (x_m/x)^\alpha$ | $x_m (1 - U)^{-1/\alpha}$ |
| Triangular$(a, m, b)$ | piecewise quadratic | $a + \sqrt{U(b-a)(m-a)}$ if $U < \frac{m-a}{b-a}$, else $b - \sqrt{(1-U)(b-a)(b-m)}$ |
| Discrete $x_1 < x_2 < \dots$ | step function | smallest $x_k$ with $F(x_k) \ge U$ |
| Empirical (data) | piecewise linear through the sorted data | linear interpolation between order statistics |
Practical points:
- Discrete distributions need a search through the cumulative probabilities: linear ($O(K)$), binary ($O(\log K)$), or Walker's alias method ($O(1)$ per draw after an $O(K)$ setup) for many draws from a large table.
- Empirical distributions avoid assuming a family, but cannot produce values outside the observed range. Add a fitted tail if extremes matter.
- No closed form (normal, gamma, beta): numerical inversion (
scipy.stats.norm.ppf) keeps all the advantages at a higher cost per sample. Otherwise use the methods of the next topic. - Truncation is easy: to sample $X$ restricted to $[a, b]$, draw $U$ uniform on $[F(a), F(b)]$.
Example: four input distributions from the same uniforms.
import numpy as np
from math import gamma
rng = np.random.default_rng(3)
U = rng.random(200_000) # one uniform per variate
# 1. Exponential service, mean 2.5 min: F(x) = 1 - exp(-x / 2.5)
x_exp = -2.5 * np.log(1 - U)
# 2. Weibull pump lifetime, shape k = 2 (wears out), scale 1000 h: F(x) = 1 - exp(-(x / 1000)^2)
x_wei = 1000 * (-np.log(1 - U)) ** (1 / 2)
# 3. Discrete daily demand for a spare part
values = np.array([0, 1, 2, 3, 4])
cum = np.cumsum([0.1, 0.2, 0.4, 0.2, 0.1]) # F = [0.1, 0.3, 0.7, 0.9, 1.0]
x_dis = values[np.searchsorted(cum, U, side="right")] # smallest k with F(k) > U
# 4. Empirical distribution of 12 observed repair times (piecewise-linear cdf)
obs = np.sort([12, 15, 19, 22, 22, 25, 31, 38, 41, 55, 63, 90])
x_emp = np.interp(U, np.linspace(0, 1, len(obs)), obs) # F^{-1}(u) by linear interpolation
print(f"exponential mean {x_exp.mean():7.2f} (theory {2.5:7.2f}) P(X > 10) {np.mean(x_exp > 10):.4f} (theory {np.exp(-4):.4f})")
print(f"Weibull mean {x_wei.mean():7.2f} (theory {1000 * gamma(1.5):7.2f}) P(X > 2000) {np.mean(x_wei > 2000):.4f} (theory {np.exp(-4):.4f})")
print("discrete frequencies", np.round(np.bincount(x_dis) / len(U), 3), " (theory [0.1 0.2 0.4 0.2 0.1])")
print(f"empirical mean {x_emp.mean():7.2f} min {x_emp.min():.1f} max {x_emp.max():.1f}")
# exponential mean 2.50 (theory 2.50) P(X > 10) 0.0184 (theory 0.0183)
# Weibull mean 886.37 (theory 886.23) P(X > 2000) 0.0184 (theory 0.0183)
# discrete frequencies [0.1 0.198 0.401 0.2 0.1 ] (theory [0.1 0.2 0.4 0.2 0.1])
# empirical mean 34.74 min 12.0 max 90.0
All four generators reproduce their target distributions. The two tail estimates are identical (0.0184) because both events mean the same thing for the shared uniform: $X_{\exp} > 10$ and $X_{\text{Wei}} > 2000$ both happen exactly when $U > 1 - e^{-4}$. Inverse transform couples outputs to inputs monotonically. Variance reduction exploits this in Section 12.2. The empirical generator never leaves the range $[12, 90]$ of the observed repair times, however many samples are drawn.
Key takeaways
- $X = F^{-1}(U)$ has cdf $F$. For discrete and empirical distributions, use the generalized inverse (search or interpolation).
- Closed forms exist for the uniform, exponential, Weibull, Pareto and triangular distributions. Other distributions need numerical inversion or other methods.
- One uniform per variate and a monotone mapping make inverse transform the best choice for common random numbers and antithetic variates.
- Empirical distributions are assumption-free but cannot generate values beyond the observed extremes.
Practice
- (Beginner) A random variable has density $f(x) = 2x$ on $[0, 1]$. Derive its inverse-transform generator and compute the samples for $U = 0.25$ and $U = 0.64$.
- (Intermediate) Implement the triangular generator from the table for an activity duration with $a = 10$, $m = 15$, $b = 30$. Check the sample mean against $(a + m + b)/3$ and the sample histogram against
rng.triangular(10, 15, 30, n).- (Interview) "Why do simulation packages generate exponential service times by inverse transform, even though faster methods exist?" Hint: think about what happens to two simulated systems that use the same seed but different service-time means, and about truncated distributions.
Acceptance-Rejection Method
Why it matters. Many important distributions have no closed-form inverse cdf: the normal, gamma and beta distributions, and custom densities fitted to data. Acceptance-rejection samples them using only the density. The same idea, called thinning, generates arrivals whose rate changes during the day, such as a morning rush followed by a lunch peak. Constant-rate arrivals miss these peaks, so thinning is essential for realistic service models.
Theory. Let $f$ be the target density and $g$ a proposal density that is easy to sample, with a constant $c \ge 1$ such that $f(x) \le c,g(x)$ for all $x$ (the envelope). The algorithm is:
repeat
draw Y ~ g and U ~ Uniform(0, 1), independently
until U <= f(Y) / (c * g(Y))
return X = Y
The points $(Y, U c g(Y))$ are uniform under the curve $c,g$. The accepted ones are uniform under $f$, so their $x$-coordinate has density $f$.
density
| . . . . . .
| . . <- envelope c g(x)
| . x x .
| . _______ x . x rejected point (above f)
|. / o o \ x . o accepted point (below f)
| / o o o \___ x . .
|/_o___o___o___o____\________. . _ x
f(x)
Each trial is accepted with probability $1/c$, so the number of trials is geometric with mean $c$. A good proposal has the same shape as the target, which gives $c$ close to 1. Only the ratio $f/g$ is needed, so $f$ may be known up to a normalizing constant. The same property makes acceptance-rejection a building block of Bayesian sampling.
Other general methods complete the toolbox:
| Method | Idea | Example |
|---|---|---|
| Composition | mixture $f = \sum_i p_i f_i$: pick component $i$ with probability $p_i$, then sample $f_i$ | hyperexponential service times |
| Convolution | sum of independent variables | Erlang = sum of $k$ exponentials, binomial = sum of Bernoullis |
| Box-Muller | two uniforms give two independent normals | $\sqrt{-2\ln U_1}\cos(2\pi U_2)$ |
| Ziggurat | refined acceptance-rejection with layered rectangles | NumPy's normal and exponential generators |
| Thinning | acceptance-rejection in time | non-homogeneous Poisson arrivals |
Thinning (Lewis and Shedler, 1979) generates a non-homogeneous Poisson process with rate $\lambda(t) \le \lambda^$. Generate candidate events from a homogeneous Poisson process with rate $\lambda^$, and keep a candidate at time $t$ with probability $\lambda(t)/\lambda^$. The fraction of candidates kept is the average of $\lambda(t)$ divided by $\lambda^$, so $\lambda^*$ should be a tight bound.
Example: normal variates and a coffee bar day. First, a standard normal is built from a half-normal target with an Exp(1) proposal, for which $c = \sqrt{2e/\pi}$ and the acceptance test simplifies to $U \le e^{-(Y-1)^2/2}$. Second, thinning generates one day of coffee bar arrivals from 7:00 to 15:00, with a base rate of 12 customers per hour, a morning rush peaking near 40 per hour at 8:00 and a lunch peak at 12:30. The rate function and its bound are reused by the discrete-event simulation of the next topic.
import numpy as np
from math import erf, sqrt, pi, e
rng = np.random.default_rng(5)
# 1. Standard normal by acceptance-rejection
# target: half-normal f(x) = sqrt(2/pi) exp(-x^2/2), x >= 0
# proposal: exponential g(x) = exp(-x), with f <= c g for c = sqrt(2e/pi)
def normal_ar(n):
out, tries = [], 0
while len(out) < n:
y = -np.log(1 - rng.random()) # proposal Y ~ Exp(1) by inverse transform
tries += 1
if rng.random() <= np.exp(-(y - 1) ** 2 / 2): # f(y) / (c g(y)) simplifies to this
out.append(y if rng.random() < 0.5 else -y) # random sign: half-normal -> normal
return np.array(out), tries
z, tries = normal_ar(100_000)
print(f"c = {sqrt(2 * e / pi):.4f}, acceptance rate {len(z) / tries:.4f} (theory {1 / sqrt(2 * e / pi):.4f})")
print(f"mean {z.mean():+.4f} var {z.var():.4f} P(Z > 2) {np.mean(z > 2):.4f} (theory {0.5 * (1 - erf(2 / sqrt(2))):.4f})")
# 2. Non-homogeneous Poisson arrivals by thinning (a coffee bar day, 7:00-15:00)
def rate(t):
"""Arrivals per minute, t = minutes after 7:00: base 12/h, rush at 8:00, lunch peak at 12:30."""
return (12 + 28 * np.exp(-((t - 60) / 30) ** 2) + 18 * np.exp(-((t - 330) / 40) ** 2)) / 60
LAM_MAX, T_END = 41 / 60, 480.0 # LAM_MAX >= rate(t) for every t
def arrivals_thinning(rng):
t, out = 0.0, []
while True:
t += rng.exponential(1 / LAM_MAX) # candidate from a Poisson process of rate LAM_MAX
if t > T_END:
return np.array(out)
if rng.random() <= rate(t) / LAM_MAX: # keep it with probability rate(t) / LAM_MAX
out.append(t)
days = [arrivals_thinning(rng) for _ in range(2_000)]
grid = np.linspace(0, T_END, 48_001)
expected = rate(grid).mean() * T_END # integral of rate(t) over the day
print(f"customers/day: simulated {np.mean([len(d) for d in days]):.1f}, expected {expected:.1f}")
hours = np.bincount(np.concatenate(days).astype(int) // 60, minlength=8) / len(days)
print("mean arrivals per hour (7-8, 8-9, ..., 14-15):", np.round(hours, 1))
# c = 1.3155, acceptance rate 0.7605 (theory 0.7602)
# mean +0.0019 var 1.0013 P(Z > 2) 0.0230 (theory 0.0228)
# customers/day: simulated 141.9, expected 142.0
# mean arrivals per hour (7-8, 8-9, ..., 14-15): [24.4 24.2 12.1 11.9 15. 27. 15.2 12.1]
The normal generator accepts 76% of its proposals, exactly $1/c$, and reproduces the mean, variance and tail of $N(0, 1)$. Thinning produces 142 customers per day as the integral of $\lambda(t)$ predicts, with the rush visible at 7:00 to 9:00 and the lunch peak at 12:00 to 13:00. Only 43% of the candidates are kept ($142$ out of $\lambda^* \cdot 480 = 328$), which costs little because candidates are cheap. A model with a constant 17.75 customers per hour would match the daily total but miss both peaks, and the peaks are what create queues.
Key takeaways
- Acceptance-rejection samples a density $f$ through a proposal $g$ with $f \le c,g$. Each trial succeeds with probability $1/c$, so choose $g$ to resemble $f$.
- Only $f/g$ is needed, so $f$ can be unnormalized.
- Composition, convolution, Box-Muller and the ziggurat complete the toolbox for standard distributions.
- Thinning generates time-varying Poisson arrivals: generate at the peak rate $\lambda^$ and keep each event with probability $\lambda(t)/\lambda^$.
Practice
- (Beginner) Sample from the density $f(x) = 6x(1-x)$ on $[0, 1]$ (a Beta(2, 2) distribution) with a uniform proposal. Find the smallest $c$, the acceptance probability and the expected number of trials per sample.
- (Intermediate) Sample a gamma distribution with shape 2.5 and scale 1 using an exponential proposal with mean 2.5. Compute $c = \max_x f(x)/g(x)$ numerically on a grid, implement the sampler, and compare its acceptance rate, mean and variance with theory. Why would an exponential proposal with mean 1 fail?
- (Interview) "Our call-center model generates Poisson arrivals at the daily average rate, and it says we are overstaffed. Agents say they are overwhelmed every morning." What is wrong, and how would you fix the model? Hint: averaging the rate removes the peaks. Consider a non-homogeneous Poisson process generated by thinning, and day-to-day variation in the rate itself.
Discrete-Event Simulation
Why it matters. Monte Carlo handles static questions. Queues, factories, hospitals and warehouses evolve over time, with entities competing for resources and decisions depending on the current state. Discrete-event simulation (DES) is the standard way to model them, and it is the engine inside commercial tools such as Arena, Simio, AnyLogic and FlexSim and the Python library SimPy. In these systems the state changes only at discrete instants, the events, so the simulator jumps from one event to the next instead of ticking through every second.
Theory. A DES model has these components:
| Component | Meaning | Coffee bar |
|---|---|---|
| Entities | objects that flow through the system | customers |
| Resources | limited capacity that entities compete for | baristas |
| State | variables that describe the system at time $t$ | busy baristas, waiting customers |
| Events | instantaneous changes of state | arrival, departure |
| Simulation clock | current simulated time | minutes after 7:00 |
| Future event list (FEL) | scheduled events, ordered by time | heap of (time, event type) |
| Statistical counters | accumulated outputs | waits, area under the queue length |
The next-event time advance loop drives everything:
initialize state, counters, clock = 0; schedule the first arrival
while the FEL is not empty (and no stop condition holds):
remove the earliest event (t, type) from the FEL
update time-persistent counters for the interval [clock, t]
clock = t
run the event routine: change the state, schedule future events, record observations
report statistics
For a FCFS queue with $c$ servers, the two event routines are:
- Arrival: schedule the next arrival. If a server is free, make it busy and schedule this customer's departure (wait 0). Otherwise put the customer at the end of the queue.
- Departure: if the queue is not empty, take the first waiting customer, record its wait and schedule its departure. Otherwise the server becomes idle.
Outputs come in two kinds:
- Observation-based (tally) statistics average over entities, such as the mean wait over customers: $\frac{1}{N}\sum_i W_{q,i}$.
- Time-persistent statistics average over time, such as the mean queue length $\frac{1}{T}\int_0^T Q(t),dt$. Because $Q(t)$ is a step function that only changes at events, the integral is a sum of rectangles accumulated just before each event:
Q(t)
2 | +-------+
1 | +--------+ +----+
0 |------+ +---------
0 3.1 5.0 7.2 8.0 t
area = 1(5.0-3.1) + 2(7.2-5.0) + 1(8.0-7.2) = 7.1
Three world views organize the same logic in different ways:
| World view | The modeler writes | Used by |
|---|---|---|
| Event scheduling | one routine per event type that updates the state | hand-coded simulators, this example |
| Process interaction | the life of an entity as a sequence of steps (arrive, request barista, hold, release) | SimPy, Arena, Simio, most commercial tools |
| Activity scanning (three-phase) | conditions that start activities | some manufacturing tools |
Fixed-increment time advance (ticking the clock by $\Delta t$) is the alternative. It wastes effort when events are sparse and introduces discretization errors and artificial ties. It is the right choice for continuous dynamics (system dynamics) and for many agent-based models.
Example: one day at the coffee bar. The bar opens at 7:00 and closes its doors at 15:00, and customers inside are still served. Arrivals follow the rush and lunch pattern of the previous topic (thinning), and service times are exponential with mean 2.5 minutes. The simulator below uses event scheduling with a binary heap as the FEL. Later topics reuse coffee_bar and erlang_c_wq.
import heapq
from collections import deque
from math import factorial
import numpy as np
# rate, LAM_MAX, T_END: defined in the Acceptance-Rejection Method example
def coffee_bar(c, seed, rate=rate, lam_max=LAM_MAX, t_end=T_END,
service=lambda g: g.exponential(2.5), trace=0):
"""Event-scheduling DES of a FCFS coffee bar with c baristas. Doors close at t_end;
customers already inside are still served. Times in minutes after 7:00."""
rng = np.random.default_rng(seed)
fel, seq = [], 0 # future event list: heap of (time, seq, kind)
def schedule(t, kind):
nonlocal seq
heapq.heappush(fel, (t, seq, kind)); seq += 1 # seq breaks ties in FIFO order
def next_arrival(t): # thinning, as in the previous topic
while True:
t += rng.exponential(1 / lam_max)
if t > t_end:
return None
if rng.random() <= rate(t) / lam_max:
return t
clock, busy, queue = 0.0, 0, deque() # state: clock, busy baristas, waiting customers
area_q, waits, arr, max_q, n_ev = 0.0, [], [], 0, 0 # statistical counters
schedule(next_arrival(0.0), "arrival")
while fel:
t, _, kind = heapq.heappop(fel) # next event: jump the clock straight to it
area_q += len(queue) * (t - clock) # integral of Q(t): time-persistent statistic
clock = t
if kind == "arrival":
nxt = next_arrival(clock)
if nxt is not None:
schedule(nxt, "arrival") # bootstrapping: each arrival schedules the next
if busy < c:
busy += 1; waits.append(0.0); arr.append(clock) # free barista: no wait
schedule(clock + service(rng), "departure")
else:
queue.append(clock) # remember the arrival time to measure the wait
max_q = max(max_q, len(queue))
else: # departure frees a barista
if queue:
a = queue.popleft()
waits.append(clock - a); arr.append(a)
schedule(clock + service(rng), "departure")
else:
busy -= 1
n_ev += 1
if n_ev <= trace:
print(f" {7 + int(clock // 60):02d}:{clock % 60:04.1f} {kind:9s} busy={busy} queue={len(queue)}")
w, a = np.array(waits), np.array(arr)
rush = (a >= 30) & (a < 90) # customers arriving 7:30-8:30
return dict(served=len(w), wq=w.mean(), wq_rush=w[rush].mean(), p5=np.mean(w > 5),
lq=area_q / clock, max_q=max_q, area_q=area_q, sum_w=w.sum())
print("Trace of the first events, 1 barista:")
coffee_bar(1, seed=3, trace=8)
def erlang_c_wq(lam, mu, c): # steady-state M/M/c wait (Section 9.3)
a, rho = lam / mu, lam / (c * mu)
tail = a ** c / factorial(c) / (1 - rho)
return tail / (sum(a ** k / factorial(k) for k in range(c)) + tail) / (c * mu - lam)
print("\nbaristas Wq Wq 7:30-8:30 P(wait>5) Lq max queue | M/M/c formula, mean rate")
for c in (1, 2, 3):
runs = [coffee_bar(c, seed) for seed in range(200)] # 200 independent days
avg = {k: np.mean([r[k] for r in runs]) for k in runs[0]}
print(f"{c:8d} {avg['wq']:5.2f} {avg['wq_rush']:10.2f} {avg['p5']:9.3f} {avg['lq']:5.2f}"
f" {avg['max_q']:9.1f} | Wq = {erlang_c_wq(142 / 480, 1 / 2.5, c):.2f}")
# Trace of the first events, 1 barista:
# 07:00.2 arrival busy=1 queue=0
# 07:02.7 departure busy=0 queue=0
# 07:03.3 arrival busy=1 queue=0
# 07:07.1 departure busy=0 queue=0
# 07:13.2 arrival busy=1 queue=0
# 07:13.2 departure busy=0 queue=0
# 07:13.6 arrival busy=1 queue=0
# 07:15.3 arrival busy=1 queue=1
#
# baristas Wq Wq 7:30-8:30 P(wait>5) Lq max queue | M/M/c formula, mean rate
# 1 12.68 18.16 0.563 3.84 16.5 | Wq = 7.10
# 2 0.92 2.09 0.058 0.28 5.9 | Wq = 0.40
# 3 0.12 0.29 0.003 0.04 2.8 | Wq = 0.05
The trace shows the clock jumping from event to event: nothing is computed between 07:07.1 and 07:13.2. The customer who arrived at 07:13.2 had a very short service and left within the same tenth of a minute.
The steady-state M/M/c formula, fed with the daily mean rate of 17.75 customers per hour, understates the average wait by a factor of about 2 at every staffing level. It says nothing about the rush, where the wait with two baristas is 2.1 minutes, five times the formula's value. Delay grows convexly with load, so the peaks contribute far more waiting than the quiet hours save. With one barista the morning rush exceeds capacity (40 against 24 customers per hour), and the longest queue of the day averages 16.5 customers. This is the kind of question that only simulation answers directly: two baristas are enough on average, and the rush-hour column suggests that a third barista during the morning rush would remove most of the remaining waiting.
Key takeaways
- DES models a system by its state, events, a simulation clock and a future event list. The clock jumps from event to event.
- Event routines change the state and schedule future events. Arrivals schedule the next arrival (bootstrapping).
- Observation-based statistics average over entities. Time-persistent statistics integrate the state over time, accumulated as rectangles between events.
- Event scheduling, process interaction and activity scanning are three ways to write the same logic. Most commercial tools use process interaction.
- Time-varying, non-exponential systems are where DES beats steady-state formulas.
Practice
- (Beginner) Simulate by hand a single-server FCFS queue with arrivals at times 0, 1, 2, 6 and 7 and service times 3, 2, 4, 1 and 2. Build the event table (clock, event, queue length, server status, FEL) and compute the mean wait and the time-average queue length until the last departure.
- (Intermediate) Extend
coffee_barwith a third barista who works only from 7:30 to 9:00 (add "shift start" and "shift end" events), and with balking: an arriving customer who sees 6 or more people waiting leaves. Report the mean wait, the rush-hour wait and the number of lost customers per day for both extensions.- (Interview) "Why not just advance the simulation clock in steps of one second and update everything at each step?" Hint: compare the work done when events are sparse, think about events that fall between ticks and about artificial ties, and name the kinds of models for which time stepping is the natural choice.
<a id="122-simulation-analysis"></a>
12.2 Simulation Analysis
A simulation is a statistical experiment: every output is a random variable, and a single run is a single sample. This subchapter covers the analysis around the model: building input distributions from data, turning outputs into confidence intervals, handling the warm-up period of steady-state studies, getting more precision per run, checking that the model is right, and using simulation to optimize decisions. The examples keep using the coffee bar, both the M/M/1 version with its exact $W_q = 12.5$ minutes and the coffee_bar simulator of Section 12.1.
Input Data Modeling
Why it matters. A simulation can only be as good as its inputs. A model with perfect logic but the wrong service-time distribution gives precise answers to the wrong question. Section 9.3 showed that queue lengths depend on the variability of service times through $c^2$, so the shape of the distribution matters as much as its mean. Input modeling turns raw data into the distributions that drive the simulation.
Theory. Input modeling follows six steps:
- Collect the right data. Measure the quantity the model needs: the service time, not the time between departures. Aim for at least 100 observations. Watch for censored data: recorded sales understate demand when there are stockouts.
- Check the i.i.d. assumption. Plot the data over time and compute autocorrelations. Trends, shift effects and day-of-week patterns call for separate distributions per period or a time-varying model (thinning, Section 12.1).
- Hypothesize families. Use the mechanism behind the data (Section 9.1), summary statistics ($c^2$, skewness) and the histogram. Positive, right-skewed times suggest the gamma, lognormal or Weibull distributions. A $c^2$ near 1 suggests the exponential.
- Estimate parameters, usually by maximum likelihood (MLE). The exponential MLE is $1/\bar x$. The gamma, lognormal and Weibull MLEs need a numerical solve.
- Assess the fit.
- Graphically: Q-Q and P-P plots, which show where a fit fails, especially in the tails that drive queues and stockouts.
- With tests: chi-square, Kolmogorov-Smirnov (KS) and Anderson-Darling (more sensitive in the tails).
- Across candidate families: log-likelihood or $\text{AIC} = 2k - 2\ln L$, where $k$ is the number of parameters.
- Without data, use expert estimates: a triangular or PERT-beta distribution from minimum, most likely and maximum values, or a uniform distribution if only a range is known. Follow up with sensitivity analysis.
Tests need care:
- With parameters estimated from the same data, the standard KS p-values are too optimistic. Use the Lilliefors correction or a parametric bootstrap.
- With thousands of observations every family is rejected for tiny deviations, and with 20 observations almost nothing is. Use tests as a screen, not as a verdict.
| Option | Advantages | Disadvantages |
|---|---|---|
| Fitted theoretical distribution | smooth, compact, generates values beyond the data, easy to vary in what-if studies | may be the wrong family |
| Empirical distribution (resampling) | no assumption about the family | bounded by observed extremes, noisy tails |
| Trace-driven (replay the recorded data) | reproduces history exactly, best for validation | cannot explore new conditions, limited length |
Input uncertainty. Fitted parameters are estimates too. With 200 observations, the fitted standard deviation can easily be 10% off, and the simulated output inherits that error. That error does not shrink however many replications are run. The bootstrap quantifies it: resample the data, refit, rerun, and look at the spread.
Example: fitting service times. The coffee bar recorded 200 service times. The hidden truth is lognormal with mean 2.5 and standard deviation 1.5 minutes. The code fits four families and checks how each choice changes the M/G/1 waiting time given by the Pollaczek-Khinchine formula $W_q = \lambda E[S^2] / (2(1-\rho))$ (Section 9.3).
import numpy as np
from scipy import stats
# 200 observed service times (minutes). In practice they come from timestamps; here the
# hidden truth is lognormal with mean 2.5 and standard deviation 1.5.
rng = np.random.default_rng(11)
sig = np.sqrt(np.log(1 + (1.5 / 2.5) ** 2))
data = rng.lognormal(np.log(2.5) - sig ** 2 / 2, sig, 200)
print(f"n = {len(data)}, mean = {data.mean():.2f}, sd = {data.std(ddof=1):.2f}, "
f"cv^2 = {data.var(ddof=1) / data.mean() ** 2:.2f}, lag-1 autocorr = {np.corrcoef(data[:-1], data[1:])[0, 1]:+.3f}")
candidates = {"exponential": stats.expon, "gamma": stats.gamma,
"lognormal": stats.lognorm, "Weibull": stats.weibull_min}
lam, true_ES2 = 1 / 3, 2.5 ** 2 + 1.5 ** 2 # 20 customers/hour; E[S^2] = mean^2 + var
print("\ndistribution params log-lik AIC KS stat KS p-value P-K Wq (min)")
for name, dist in candidates.items():
params = dist.fit(data, floc=0) # maximum likelihood, location fixed at 0
k = len(params) - 1 # free parameters (loc is not estimated)
ll = dist.logpdf(data, *params).sum()
ks = stats.kstest(data, dist.cdf, args=params) # p-value is optimistic: params were fitted
ES, ES2 = dist.mean(*params), dist.moment(2, *params)
wq = lam * ES2 / (2 * (1 - lam * ES)) # Pollaczek-Khinchine, Section 9.3
print(f"{name:12s} {k:4d} {ll:8.1f} {2 * k - 2 * ll:6.1f} {ks.statistic:7.3f} {ks.pvalue:10.3f} {wq:9.2f}")
print(f"{'truth':12s} {'':4s} {'':8s} {'':6s} {'':7s} {'':10s} {lam * true_ES2 / (2 * (1 - lam * 2.5)):9.2f}")
# Q-Q check in the upper tail, where queues are decided
q = [0.5, 0.9, 0.95, 0.99]
print("\nquantile 0.50 0.90 0.95 0.99")
print("data ", np.round(np.quantile(data, q), 2))
for name in ("exponential", "lognormal"):
dist = candidates[name]
print(f"{name:13s}", np.round(dist.ppf(q, *dist.fit(data, floc=0)), 2))
# n = 200, mean = 2.47, sd = 1.34, cv^2 = 0.29, lag-1 autocorr = +0.052
#
# distribution params log-lik AIC KS stat KS p-value P-K Wq (min)
# exponential 1 -380.5 762.9 0.274 0.000 11.37
# gamma 2 -308.2 620.4 0.049 0.704 7.11
# lognormal 2 -304.0 612.1 0.042 0.855 7.39
# Weibull 2 -319.2 642.3 0.081 0.136 7.51
# truth 8.50
#
# quantile 0.50 0.90 0.95 0.99
# data [2.17 4.24 4.97 6.9 ]
# exponential [ 1.71 5.68 7.39 11.35]
# lognormal [2.16 4.17 5.02 7.11]
The data show no autocorrelation, and their $c^2 = 0.29$ is far below the exponential's 1. The exponential fit is rejected by every criterion: its KS statistic is 0.274, its AIC is 150 points worse, and its 99% quantile is 11.4 minutes against 6.9 in the data. The gamma and lognormal fits are both acceptable, and the lognormal has the best AIC. The choice has real consequences. The exponential assumption predicts a wait of 11.4 minutes, 54% more than the lognormal fit, which could lead to hiring a barista who is not needed. Even the best fit gives 7.4 minutes against the true 8.5, because the sample standard deviation (1.34) underestimates the true value (1.5). That gap is input uncertainty from having only 200 observations, and more replications would not remove it.
Key takeaways
- Collect the right quantity, check independence and stationarity, fit candidate families by maximum likelihood, and compare them graphically and with AIC.
- Goodness-of-fit tests are screens. Their p-values are optimistic with fitted parameters, and their power depends heavily on the sample size.
- The tails and $c^2$ of service and demand distributions drive congestion and stockouts. A wrong family can bias outputs by 50% or more.
- Fitted parameters carry input uncertainty, which more replications cannot remove. Quantify it with the bootstrap.
Practice
- (Beginner) A sample of 150 interarrival times has mean 4.0 and standard deviation 4.1 minutes. Which family is a natural first candidate, and what is its MLE? Which plot and which test would you use to check it?
- (Intermediate) Quantify input uncertainty in the example with the bootstrap: resample the 200 observations 500 times, refit the lognormal each time, compute the Pollaczek-Khinchine $W_q$, and report a 90% interval. Does it contain the true 8.50?
- (Interview) "The analyst fitted a normal distribution to three years of daily sales, the KS test passed, and the inventory simulation now uses it. Any concerns?" Hint: think about seasonality and trends, negative values, stockouts that censor demand, promotions, autocorrelation, and p-values with fitted parameters.
Output Analysis and Confidence Intervals
Why it matters. A single simulation run is a single random sample, so its output may be far from the true mean. The most common mistake is to treat the thousands of observations inside one run as independent. Waits of successive customers are strongly correlated, so such confidence intervals are far too narrow, and decisions are then taken with unjustified confidence. Output analysis gives intervals that mean what they claim.
Theory. The method depends on the type of simulation:
| Terminating simulation | Non-terminating (steady-state) simulation | |
|---|---|---|
| Ends | at a natural event (closing time, project completion) | never; interest is in long-run behavior |
| Initial conditions | part of the model (the bar opens empty) | a nuisance that biases early output |
| Measure | e.g. mean wait over a day | e.g. long-run mean wait $W_q$ |
| Method | independent replications | replication/deletion, batch means |
Independent replications (terminating). Run $R$ replications with independent streams. Each gives one output $Y_r$, such as the mean wait of one day. The $Y_r$ are i.i.d., so the $t$ interval applies:
$$ \bar Y \pm t_{R-1,,1-\alpha/2},\frac{S}{\sqrt{R}}. $$
Why within-run data fail. For a stationary output sequence with variance $\sigma^2$ and lag-$k$ autocorrelation $\rho_k$,
$$ \mathrm{Var}(\bar Y_n) = \frac{\sigma^2}{n}\Big[1 + 2\sum_{k=1}^{n-1}\Big(1 - \frac{k}{n}\Big)\rho_k\Big]. $$
The bracket can be in the hundreds for a busy queue, so the naive variance $\sigma^2/n$ is far too small. The effective number of independent observations is $n$ divided by that bracket.
Steady-state methods.
- Replication/deletion. Run $R$ replications, delete the warm-up period of each (next topic), and treat the averages of the remaining data as i.i.d. outputs. It is simple and robust, but pays for a warm-up in every replication.
- Batch means. Make one long run, delete the warm-up, and split the rest into $b$ batches of size $m$. If $m$ is much larger than the correlation length, the batch means are nearly independent and the $t$ interval with $b - 1$ degrees of freedom applies. Typical choices are 10 to 30 batches.
- Others: the regenerative method (cycles between returns to an empty system are i.i.d.), standardized time series and spectral methods.
Precision and run length. The half-width $h$ shrinks like $1/\sqrt{R}$. From a pilot with $R_0$ replications and half-width $h_0$, the target $\varepsilon$ needs about $R \approx R_0 (h_0/\varepsilon)^2$. Pilot variance estimates are noisy, especially for skewed outputs, so a sequential procedure is safer: add replications until $h \le \varepsilon$, or until the relative half-width $h/|\bar Y| \le \gamma$.
Comparing systems. Build a confidence interval for the difference. Use a paired-$t$ interval on paired outputs with common random numbers (next topics), or Welch's interval for independent runs. Several systems need a Bonferroni correction or ranking-and-selection procedures (Simulation-Based Optimization).
What to report. Report the interval, not just the point estimate. Also report service-level measures, such as $P(\text{wait} > 5)$ (estimated as proportions), and quantiles.
Example: interval coverage and a sequential procedure. Part 1 checks how often a nominal 95% interval for the M/M/1 coffee bar contains the true $W_q = 12.5$ minutes. It repeats the whole experiment 200 times. Part 2 estimates the mean daily wait of the two-barista bar of Section 12.1, a terminating simulation, to within $\pm 0.1$ minute.
import numpy as np
from scipy import stats
# coffee_bar: defined in the Discrete-Event Simulation example
def mm1_waits(n, rng, lam=1 / 3, mu=1 / 2.5, w0=0.0):
"""Queue waits of n successive M/M/1 customers (Lindley recursion), first one waits w0."""
a, s = rng.exponential(1 / lam, n), rng.exponential(1 / mu, n)
w = np.empty(n); w[0] = w0
for k in range(1, n):
w[k] = max(0.0, w[k - 1] + s[k - 1] - a[k]) # wait = previous wait + its service - gap
return w
def t_interval(x, level=0.95):
"""Mean and half-width of a t confidence interval, treating x as i.i.d."""
return x.mean(), stats.t.ppf((1 + level) / 2, len(x) - 1) * x.std(ddof=1) / np.sqrt(len(x))
# Part 1, steady state: does the 95% interval contain the true Wq = 12.5 min?
rng = np.random.default_rng(0)
hits, widths = {"naive": 0, "batch means": 0}, {"naive": [], "batch means": []}
for _ in range(200): # 200 independent experiments
w = mm1_waits(55_000, rng)[5_000:] # delete a warm-up of 5,000 customers
for name, x in (("naive", w), ("batch means", w.reshape(20, -1).mean(axis=1))):
m, h = t_interval(x) # 50,000 waits, or 20 batches of 2,500
hits[name] += abs(m - 12.5) <= h; widths[name].append(h)
print(f"lag-1 autocorrelation of successive waits: {np.corrcoef(w[:-1], w[1:])[0, 1]:.3f}")
for name in hits:
print(f"{name:12s} coverage {hits[name] / 200:5.1%} mean half-width {np.mean(widths[name]):.2f} min")
# Part 2, terminating: mean wait over a whole day, 2 baristas. Sequential procedure:
# add 10 days at a time until the 95% half-width is at most 0.1 min.
days = []
while True:
days += [coffee_bar(2, seed)["wq"] for seed in range(len(days), len(days) + 10)]
m, h = t_interval(np.array(days))
if len(days) in (10, 50) or h <= 0.1:
print(f"{len(days):4d} days: Wq = {m:.2f} +/- {h:.2f} min")
if h <= 0.1:
break
# lag-1 autocorrelation of successive waits: 0.971
# naive coverage 17.0% mean half-width 0.13 min
# batch means coverage 91.5% mean half-width 1.48 min
# 10 days: Wq = 0.62 +/- 0.28 min
# 50 days: Wq = 0.83 +/- 0.13 min
# 180 days: Wq = 0.90 +/- 0.10 min
Successive waits have a correlation of 0.97. The naive interval, which treats 50,000 waits as independent, is more than 10 times too narrow and contains the truth in only 17% of the experiments. Batch means give honest intervals of about $\pm 1.5$ minutes with 91.5% coverage, slightly below 95% because 20 batches still retain a little correlation. In Part 2, the pilot of 10 days suggested about $10 \times (0.28/0.1)^2 \approx 79$ days. Daily mean waits are skewed, though: a few bad days dominate. The pilot underestimated both the mean and the variance, and the sequential procedure needed 180 days.
Key takeaways
- Terminating simulations use independent replications. Each replication gives one i.i.d. observation for a $t$ interval.
- Observations within a run are autocorrelated. Treating them as i.i.d. gives intervals that are far too narrow.
- Steady-state studies use replication/deletion or batch means. Batches must be long compared with the correlation length.
- Pilot-based sample sizes are rough. A sequential procedure guarantees the target precision.
- Compare systems with intervals on differences, preferably paired with common random numbers.
Practice
- (Beginner) Five replications of a day give mean waits of 4.2, 5.1, 3.8, 6.0 and 4.9 minutes. Compute the 95% confidence interval ($t_{4,,0.975} = 2.776$) and estimate how many replications a half-width of 0.3 minute requires.
- (Intermediate) Repeat the coverage experiment with 5, 20 and 100 batches over the same total run length. Explain the trade-off between correlated batches (too many small ones) and few degrees of freedom (too few large ones).
- (Interview) "Our simulation report says the mean wait is 4.3 minutes. We simulated a million customers, so it is very accurate." What questions do you ask? Hint: one run or several? How was the interval computed, and was autocorrelation accounted for? Terminating or steady state? Warm-up? And statistical precision says nothing about whether the model is valid.
Warm-Up Period and Steady State
Why it matters. A steady-state simulation must start in some state, usually empty and idle, which is not typical of the long run. Early observations are therefore biased: a queue that starts empty shows short waits at first. Averaging them in biases the estimate. Discarding too much wastes data. Choosing the warm-up period is the main practical decision in a steady-state study.
Theory. Let $Y_k$ be the $k$-th output, for example the wait of customer $k$. Initialization bias means that $E[Y_k]$ differs from the steady-state mean $\nu$ for small $k$, and converges to $\nu$ as $k$ grows. The speed depends on the system. Heavily loaded queues converge slowly: for the M/M/1 queue the relaxation time grows roughly like $1/(1-\rho)^2$, so a system at 95% load needs about 9 times longer to warm up than one at 85%.
The standard remedy is deletion: discard the first $d$ observations and estimate $\nu$ from the rest,
$$ \bar Y(n, d) = \frac{1}{n - d} \sum_{k=d+1}^{n} Y_k. $$
Increasing $d$ reduces bias but increases variance, so the aim is a small mean squared error. Four approaches are common:
- Welch's graphical method. Make $R$ replications, average the $k$-th observation across replications to get $\bar Y_k$, smooth this curve with a moving average of window $w$, and choose $d$ where the curve flattens. It is simple and visual, but needs judgment.
- MSER-5 (White, 1997). Group the data into batches of 5, then choose $d$ (at most half of the data) to minimize $$ \text{MSER}(d) = \frac{1}{(n-d)^2} \sum_{k=d+1}^{n} \big(Y_k - \bar Y(n, d)\big)^2, $$ the squared standard error of the truncated mean. It is automatic and works well on averaged data, but is noisy on single runs.
- Long runs. The bias of $\bar Y_n$ decreases like $1/n$, the standard error only like $1/\sqrt{n}$. In a run far longer than the warm-up, the bias becomes negligible relative to the noise.
- Better initial conditions. Starting near a typical state (for example with the mean number in the system) shortens the warm-up. A very unusual start lengthens it.
Terminating simulations need no deletion when their initial conditions are real, such as a bar that really opens empty. A model of a hospital ward for next month must start from today's real occupancy, or from a warm-up long enough to reach a realistic state.
Example: warm-up of the M/M/1 coffee bar. 2,000 replications of 3,000 customers start either empty or crowded (the first customer waits 60 minutes, as after an equipment failure). The code applies Welch's method and MSER-5 and measures the bias of runs of 1,000 customers.
import numpy as np
# mm1_waits: defined in the Output Analysis example (true steady-state Wq = 12.5 min)
rng = np.random.default_rng(1)
R, n = 2_000, 3_000
empty = np.array([mm1_waits(n, rng) for _ in range(R)]) # start empty and idle
crowded = np.array([mm1_waits(n, rng, w0=60.0) for _ in range(R)]) # first customer waits 60 min
# Welch's method: average customer k across replications, then smooth with a moving average
def moving_avg(y, win=50):
w = lambda k: min(k, win) # Welch: shrink the window near the start
return np.array([y[k - w(k):k + w(k) + 1].mean() for k in range(len(y) - win)])
welch_e, welch_c = moving_avg(empty.mean(axis=0)), moving_avg(crowded.mean(axis=0))
print("customer k 1 25 50 100 200 400 800 1600 2950")
for name, y in (("start empty ", welch_e), ("start crowded", welch_c)):
print(name, " ".join(f"{y[k]:6.2f}" for k in (0, 24, 49, 99, 199, 399, 799, 1599, 2949)))
def mser5(y):
"""MSER-5 truncation point: batch means of 5, minimize variance of the rest / (n - d)^2."""
b = y[: len(y) // 5 * 5].reshape(-1, 5).mean(axis=1)
stat = [b[d:].var() / (len(b) - d) for d in range(len(b) // 2)] # = sum sq. dev / (n-d)^2
return 5 * int(np.argmin(stat))
print(f"\nMSER-5 cut on the Welch averages: empty {mser5(empty.mean(axis=0))}, crowded {mser5(crowded.mean(axis=0))}")
cut = [mser5(w) for w in crowded] # one truncation point per run
print(f"MSER-5 on single crowded runs: median cut {np.median(cut):.0f}, "
f"10th-90th percentile {np.percentile(cut, 10):.0f}-{np.percentile(cut, 90):.0f}")
# Bias of a short run of 1,000 customers with and without deletion
for cut_d in (0, 100, 300):
est_e, est_c = empty[:, cut_d:1_000].mean(axis=1), crowded[:, cut_d:1_000].mean(axis=1)
print(f"delete {cut_d:3d}: mean estimate {est_e.mean():5.2f} (empty start) {est_c.mean():6.2f} (crowded start) true 12.50")
# customer k 1 25 50 100 200 400 800 1600 2950
# start empty 0.00 6.85 8.74 11.06 12.23 12.45 12.30 12.64 12.09
# start crowded 60.00 48.47 39.38 25.95 15.92 13.04 12.18 12.45 12.65
#
# MSER-5 cut on the Welch averages: empty 200, crowded 540
# MSER-5 on single crowded runs: median cut 110, 10th-90th percentile 20-976
# delete 0: mean estimate 12.03 (empty start) 16.09 (crowded start) true 12.50
# delete 100: mean estimate 12.39 (empty start) 13.52 (crowded start) true 12.50
# delete 300: mean estimate 12.50 (empty start) 12.60 (crowded start) true 12.50
Starting empty, the 25th customer waits 6.9 minutes on average, about half the steady-state value. The curve reaches 12.5 only after about 200 customers. Starting crowded, the bias is positive and much larger, and it takes about 400 customers to fade. MSER-5 on the Welch averages suggests deleting 200 and 540 customers. On single runs it is erratic (from 20 to 976), which is why averaging across replications helps. Without deletion, a 1,000-customer run is biased by $-4%$ from an empty start and by $+29%$ from a crowded one. Deleting 300 customers removes both biases. A run of 50,000 customers, as in the previous topic, would make even the crowded start almost harmless.
Key takeaways
- Steady-state runs start in an atypical state. Early outputs are biased, and the bias lasts longer at high load.
- Delete a warm-up period $d$. Welch's method (averaged, smoothed curves) and MSER-5 are the standard ways to choose it.
- Single-run warm-up estimates are noisy. Use several replications, and be generous with $d$ when data are cheap.
- Long runs and realistic initial states reduce the problem. Terminating simulations with real initial conditions need no deletion.
Practice
- (Beginner) For each case, decide whether a warm-up deletion is needed: (a) a store that opens empty at 7:00, simulated until closing; (b) the long-run mean wait at a 24/7 call center; (c) a hospital ward over the next month, starting the simulation empty.
- (Intermediate) Repeat the Welch experiment for arrival rates giving $\rho = 0.9$ and $\rho = 0.95$. How does the point where the curve flattens grow? Compare its growth with $1/(1-\rho)^2$.
- (Interview) "Our factory model always runs for one simulated week, starting empty, and shows great throughput." What would you check? Hint: an empty factory has no work-in-process to block anything. Plot the output over time, compare the warm-up with the run length, and consider a longer run with deletion.
Variance Reduction Techniques
Why it matters. The precision of a simulation estimate is proportional to $\sigma/\sqrt{n}$. Halving the variance is worth as much as doubling the computing budget. Variance reduction techniques (VRTs) get that precision almost for free by using the structure of the model. They are essential when comparing alternatives, where common random numbers often give large gains, and for rare events, where importance sampling can turn an impossible experiment into an easy one.
Theory. The main techniques are:
| Technique | Idea | Requirement | Typical gain |
|---|---|---|---|
| Common random numbers (CRN) | use the same random numbers for every alternative compared | synchronized streams, monotone response | large, for differences |
| Antithetic variates | pair each run on $U$ with a run on $1 - U$ | output monotone in the inputs (inverse transform) | 1.5-3x |
| Control variates | correct the output with a correlated quantity of known mean | a control with known expectation | $1/(1 - \rho_{YC}^2)$ |
| Importance sampling | sample from a distribution that makes the rare event common, then reweight | a good proposal (a bad one increases variance) | orders of magnitude for rare events |
| Stratified sampling, Latin hypercube | force the samples to cover strata evenly | known stratum probabilities | moderate |
| Conditioning | replace an output by its conditional expectation given part of the inputs | an analytic $E[Y \mid Z]$ | moderate to large |
Common random numbers. For the difference between two systems,
$$ \mathrm{Var}(Y_1 - Y_2) = \mathrm{Var}(Y_1) + \mathrm{Var}(Y_2) - 2,\mathrm{Cov}(Y_1, Y_2). $$
Driving both systems with the same random numbers makes the covariance positive, so the noise largely cancels and the difference reflects the change in design, not luck. Synchronization matters: each source of randomness needs its own stream, so the $k$-th customer gets the same service draw in both systems (Section 12.1). Analyze CRN outputs with a paired interval on the differences.
Antithetic variates. Run the model once with uniforms $U$ and once with $1 - U$, and average the pair. With $\rho = \mathrm{Corr}(Y, Y')$,
$$ \mathrm{Var}\Big(\frac{Y + Y'}{2}\Big) = \frac{\sigma^2}{2}(1 + \rho), $$
which beats two independent runs when $\rho < 0$. A long service time in one run becomes a short one in its partner.
Control variates. Let $C$ be a quantity measured in the same run with known mean $\mu_C$, such as the average simulated service time. Then
$$ Y_{cv} = Y - \beta (C - \mu_C), \qquad \beta^* = \frac{\mathrm{Cov}(Y, C)}{\mathrm{Var}(C)}, \qquad \mathrm{Var}(Y_{cv}) = (1 - \rho_{YC}^2),\mathrm{Var}(Y). $$
If a run happened to draw unusually long services, $C$ is above its mean, and the correction lowers $Y$ accordingly. Estimating $\beta$ from the same runs introduces a small bias, negligible for large $n$.
Importance sampling. For a rare event, sample $X$ from a proposal density $h$ instead of $f$ and reweight:
$$ \theta = E_f[g(X)] = E_h\Big[g(X),\frac{f(X)}{h(X)}\Big]. $$
A good $h$ makes the event frequent while keeping the weights $f/h$ stable. For sums of exponentials, exponential tilting (a larger mean) is the classic choice.
Example: four VRTs at the coffee bar. The output is the mean wait of the first 100 customers after opening (a terminating measure), generated by inverse transform so that antithetic variates are possible. The control variate is the average service time minus the average interarrival time, whose mean $2.5 - 3 = -0.5$ is known. CRN compares the current espresso machine with a faster one (mean service 2.2 minutes). Importance sampling estimates the probability that 100 orders take more than 350 minutes of barista time.
import numpy as np
from scipy import stats
def avg_wait(ua, us, mean_s=2.5, mean_a=3.0):
"""Mean queue wait of the first 100 customers (start empty), one row per replication.
Inputs are uniforms, turned into exponentials by inverse transform."""
a, s = -mean_a * np.log(1 - ua), -mean_s * np.log(1 - us)
w = np.zeros_like(a)
for k in range(1, a.shape[1]): # Lindley, vectorized over replications
w[:, k] = np.maximum(0, w[:, k - 1] + s[:, k - 1] - a[:, k])
return w.mean(axis=1), s.mean(axis=1) - a.mean(axis=1) # output and control variate
rng = np.random.default_rng(4)
N, K = 2_000, 100 # budget: 2,000 runs of 100 customers
# Crude Monte Carlo
y, _ = avg_wait(rng.random((N, K)), rng.random((N, K)))
crude_var = y.var(ddof=1) / N
# Antithetic variates: 1,000 pairs (U, 1 - U), each pair averaged
ua, us = rng.random((N // 2, K)), rng.random((N // 2, K))
pair = (avg_wait(ua, us)[0] + avg_wait(1 - ua, 1 - us)[0]) / 2
anti_var = pair.var(ddof=1) / (N // 2)
# Control variate: C = mean service - mean interarrival, E[C] = 2.5 - 3 = -0.5 is known
y2, c = avg_wait(rng.random((N, K)), rng.random((N, K)))
beta = np.cov(y2, c)[0, 1] / c.var(ddof=1) # optimal coefficient, estimated
ycv = y2 - beta * (c - (-0.5))
cv_var = ycv.var(ddof=1) / N
print("method estimate std. error variance reduction")
for name, est, v in (("crude", y.mean(), crude_var), ("antithetic", pair.mean(), anti_var),
("control var.", ycv.mean(), cv_var)):
print(f"{name:12s} {est:8.3f} {np.sqrt(v):10.3f} {crude_var / v:8.1f}x")
print(f"corr(output, control) = {np.corrcoef(y2, c)[0, 1]:.3f}, beta = {beta:.2f}")
# Common random numbers: is a faster machine (mean service 2.2 min) worth it?
ua, us = rng.random((N, K)), rng.random((N, K))
d_crn = avg_wait(ua, us)[0] - avg_wait(ua, us, mean_s=2.2)[0] # same uniforms for both
d_ind = avg_wait(ua, us)[0] - avg_wait(rng.random((N, K)), rng.random((N, K)), mean_s=2.2)[0]
for name, diff in (("independent", d_ind), ("CRN", d_crn)):
print(f"{name:12s} saving {diff.mean():.3f} +/- {1.96 * diff.std(ddof=1) / np.sqrt(N):.3f} min")
# Importance sampling: P(100 orders take more than 350 min), service ~ Exp(mean 2.5)
exact = stats.gamma.sf(350, 100, scale=2.5)
n = 10_000
x = rng.exponential(2.5, (n, K)).sum(axis=1)
print(f"\nP(total > 350): exact {exact:.3e} crude {np.mean(x > 350):.3e}")
z = rng.exponential(3.5, (n, K)) # sample from a tilted law, mean 3.5
lr = np.exp(K * np.log(3.5 / 2.5) - z.sum(axis=1) * (1 / 2.5 - 1 / 3.5)) # likelihood ratio f/g
est = (z.sum(axis=1) > 350) * lr
print(f"importance sampling {est.mean():.3e} +/- {1.96 * est.std(ddof=1) / np.sqrt(n):.1e}")
# method estimate std. error variance reduction
# crude 8.908 0.150 1.0x
# antithetic 8.729 0.119 1.6x
# control var. 8.732 0.089 2.8x
# corr(output, control) = 0.754, beta = 11.72
# independent saving 3.573 +/- 0.321 min
# CRN saving 3.631 +/- 0.132 min
#
# P(total > 350): exact 1.611e-04 crude 1.000e-04
# importance sampling 1.646e-04 +/- 6.9e-06
With the same budget of 2,000 runs, antithetic variates cut the variance by a factor of 1.6 and the control variate by 2.8. Theory predicts $1/(1 - 0.754^2) = 2.3$ for the control variate, and the difference is noise in the variance estimates. The mean wait of the first 100 customers is about 8.7 minutes, below the steady-state 12.5 because the bar opens empty (previous topic). CRN shrinks the half-width of the saving from 0.32 to 0.13 minute, a sixfold variance reduction, the same as running six times as many independent replications. The faster machine saves about 3.6 minutes per customer in the first 100. For the rare event, crude sampling saw a single hit in 10,000 trials, an estimate with about 80% relative error. Importance sampling with the same 10,000 samples estimates $1.65 \times 10^{-4}$ within $\pm 4%$, consistent with the exact $1.61 \times 10^{-4}$.
Key takeaways
- VRTs buy precision with structure instead of computing time: CRN for comparisons, antithetics and control variates for single estimates, importance sampling for rare events.
- CRN needs synchronized streams and a paired analysis. Antithetics need outputs that are monotone in the inputs, which inverse transform provides.
- A control variate with correlation $\rho$ to the output cuts variance by $1 - \rho^2$. Input averages with known means are natural controls.
- Importance sampling reweights by the likelihood ratio. A poorly chosen proposal can make variance worse, so always check the estimated error.
Practice
- (Beginner) An output has variance 4. Its antithetic partner has correlation $-0.5$ with it. Compute the variance of the antithetic pair average and of the average of two independent runs. How many independent runs is one antithetic pair worth?
- (Intermediate)
coffee_bardraws arrivals and service times from a single generator, so changing the number of baristas desynchronizes the streams. Give arrivals and services separate generators (spawned from one seed), then compare 2 against 3 baristas over 100 days with and without CRN. Report the half-width of the paired difference in each case.- (Interview) "With CRN, the new layout is 0.2 minute better with a tight confidence interval. With independent seeds, the interval includes zero. Which result should we trust?" Hint: both estimate the same difference. Ask which interval was built on paired differences, and whether the streams were really synchronized.
Model Verification and Validation
Why it matters. A simulation always produces precise-looking numbers, whether it is right or not. Verification asks "did we build the model right?": does the program implement the intended logic and assumptions? Validation asks "did we build the right model?": does the model represent the real system accurately enough for the decision at hand? Together with credibility, the decision makers' trust in the model, they decide whether a simulation study is worth anything.
Theory. Sargent's framework relates the real system, the conceptual model and the computer model:
+---------------------+
+--------->| Real system |<----------+
| +---------------------+ |
conceptual model operational validation
validation (compare outputs with
(are the assumptions real-system data)
and logic right?) |
| |
+---------------------+ verification +---------------------+
| Conceptual model |----------------->| Computer model |
| assumptions, logic | (is the code | simulation program |
+---------------------+ correct?) +---------------------+
Verification techniques:
- Build and test the model in modules, and use code review and version control.
- Use traces of the event sequence (Section 12.1) and animation to follow individual entities.
- Run degenerate and extreme cases: no arrivals, infinitely many servers, deterministic times, where the answer is known.
- Compare against analytic models in the special cases where they apply (M/M/c, Jackson networks), with a statistical tolerance.
- Check conservation laws and identities: arrivals = departures + entities still inside. For a system that starts and ends empty, $\int_0^T Q(t),dt$ equals the sum of all queue waits exactly. This is the sample-path form of Little's law.
- Use assertions in the code, and fixed seeds so that a failure can be reproduced.
Validation techniques:
- Face validity: experts who know the system examine the logic, animation and outputs.
- Validation of assumptions: input modeling (earlier topic) and structural assumptions such as the queue discipline and staffing rules.
- Operational (results) validation: run the model under the conditions of a historical period and compare its outputs with the real data, using confidence intervals or tests on the difference. A Turing test asks experts whether they can tell real reports from simulated ones.
- Sensitivity analysis: check that conclusions are stable over plausible ranges of the uncertain inputs.
- Predictive validation: predict a future period, then compare.
Validation is always relative to a purpose and to a range of conditions. A model validated at today's demand is not automatically valid at double the demand. Failing to reject the model with little data is weak evidence. A statistically significant difference may still be too small to matter for the decision.
Example: verifying and validating the coffee bar simulator. Three verification checks test coffee_bar against known answers: an exact identity, a degenerate case and analytic queueing formulas. For validation, 20 days of the "real" bar are observed, where service is really lognormal (the truth of the Input Data Modeling example). Two models are compared with them: the exponential model of Section 12.1 and a model using the fitted lognormal.
import numpy as np
from scipy import stats
# coffee_bar, erlang_c_wq: defined in the Discrete-Event Simulation example
# ---- Verification: is the program a correct implementation of the model? ----
r = coffee_bar(2, seed=0)
assert abs(r["area_q"] - r["sum_w"]) < 1e-6 * r["sum_w"] # integral of Q(t) = sum of waits
print("1. area under Q(t) = sum of waits:", round(r["area_q"], 3), round(r["sum_w"], 3))
r = coffee_bar(200, seed=0) # more baristas than customers ever present
print("2. 200 baristas -> mean wait", r["wq"])
const = dict(rate=lambda t: 1 / 3, lam_max=1 / 3, t_end=100_000.0) # constant 20/h, no rush
for c, mean_s in ((1, 2.5), (2, 5.0)): # same load rho = 0.83, other formula
wq = np.array([coffee_bar(c, seed, service=lambda g, m=mean_s: g.exponential(m), **const)["wq"]
for seed in range(10)]) # 10 independent long runs
half = stats.t.ppf(0.975, 9) * wq.std(ddof=1) / np.sqrt(10)
exact = erlang_c_wq(1 / 3, 1 / mean_s, c)
print(f"3. M/M/{c}: simulated Wq {wq.mean():5.2f} +/- {half:.2f}, formula {exact:5.2f} -> "
f"{'pass' if abs(wq.mean() - exact) <= half else 'INVESTIGATE'}")
# ---- Validation: is the model a good enough description of the real coffee bar? ----
# The real bar's service is lognormal (mean 2.5, sd 1.5). 20 days of its daily mean wait are observed.
sig = np.sqrt(np.log(1 + 0.6 ** 2))
real = lambda g: g.lognormal(np.log(2.5) - sig ** 2 / 2, sig)
observed = np.array([coffee_bar(2, seed, service=real)["wq"] for seed in range(1000, 1020)])
fit_sig = np.sqrt(np.log(1 + (1.34 / 2.47) ** 2)) # lognormal fitted in Input Data Modeling
models = {"exponential(2.5)": lambda g: g.exponential(2.5),
"fitted lognormal": lambda g: g.lognormal(np.log(2.47) - fit_sig ** 2 / 2, fit_sig)}
print(f"\nobserved: mean daily wait {observed.mean():.2f} min (20 days)")
for name, service in models.items():
sim = np.array([coffee_bar(2, seed, service=service)["wq"] for seed in range(200)])
t = stats.ttest_ind(sim, observed, equal_var=False) # Welch two-sample t-test
print(f"model {name}: {sim.mean():.2f} min, difference {sim.mean() - observed.mean():+.2f}, p-value {t.pvalue:.3f}")
# 1. area under Q(t) = sum of waits: 58.89 58.89
# 2. 200 baristas -> mean wait 0.0
# 3. M/M/1: simulated Wq 13.01 +/- 0.70, formula 12.50 -> pass
# 3. M/M/2: simulated Wq 11.43 +/- 0.68, formula 11.36 -> pass
#
# observed: mean daily wait 0.62 min (20 days)
# model exponential(2.5): 0.92 min, difference +0.30, p-value 0.002
# model fitted lognormal: 0.60 min, difference -0.02, p-value 0.801
All three verification checks pass. The area under $Q(t)$ equals the sum of the waits to machine precision, so the time-persistent and observation-based counters agree. With 200 baristas nobody waits. With constant arrivals, the simulator reproduces the M/M/1 and M/M/2 formulas within statistical tolerance. That exercises the thinning, the event list and the multi-server logic all at once. The exponential model passes verification but fails validation: it predicts 0.92 minutes against 0.62 observed, a 48% overestimate that 20 days of data detect clearly ($p = 0.002$). With the fitted lognormal, the difference shrinks to 0.02 minute and is not significant. The program was correct all along. The error was in an input assumption, which only validation could reveal.
Key takeaways
- Verification checks the code against the conceptual model. Validation checks the model against reality, for a specific purpose.
- Verify with traces, degenerate cases, analytic special cases (with statistical tolerance), conservation laws and assertions.
- Validate assumptions and results: compare outputs with real data under the same conditions, using intervals or tests on the difference.
- A model can be correctly programmed and still invalid. Validity holds only for the conditions tested.
Practice
- (Beginner) Classify each activity as verification or validation: (a) checking that customers entered = customers served + customers still inside; (b) showing the animation to the store manager; (c) comparing the simulated mean wait with last month's measurements; (d) running with deterministic service times and checking the waits by hand.
- (Intermediate) Plant a bug in
coffee_bar: in the departure routine, start the next service without appending the customer's wait towaits. Which of the three verification checks catches it? Then plant a second bug, wherebusyis not decremented when the queue is empty, and design a conservation check (customers arrived = customers served) that catches it.- (Interview) "Our model matched last year's average wait within 2%, so it is validated. Now we will use it to evaluate doubling the demand." Is this sound? Hint: the model was validated under particular conditions, so doubling demand is an extrapolation. A match on averages can hide errors in distributions or compensating errors. Consider validating several measures and other periods.
Simulation-Based Optimization
Why it matters. A simulation evaluates a design that someone proposes. Managers want the best design: how many baristas per shift, which $(s, S)$ levels, how large each buffer should be. Simulation-based optimization searches over designs when the objective is an expectation that can only be estimated with noise. It is how simulation turns from a what-if tool into a decision tool.
Theory. The problem is
$$ \min_{x \in X} ; f(x) = E[F(x, \xi)], $$
where $F(x, \xi)$ is the output of one simulation run with decision $x$ and random inputs $\xi$. Four difficulties set it apart from deterministic optimization:
- Noise. Each evaluation is an estimate, so a design can look best just by luck.
- Cost. One evaluation may take minutes or hours.
- No gradients. Derivatives must be estimated, or avoided.
- Optimizer's curse. The minimum of noisy estimates is biased low. The selected design looks better than it is, and its estimated cost from the search should not be reported.
| Method | Decision space | Idea |
|---|---|---|
| Ranking and selection (Rinott, KN, OCBA) | a few to hundreds of alternatives | allocate replications so that $P(\text{correct selection}) \ge 1 - \alpha$ whenever the best is at least $\delta$ better |
| Sample average approximation (SAA) | any; uses a deterministic solver | fix $n$ scenarios (CRN) and optimize the sample average exactly (Section 9.4) |
| Stochastic approximation (Robbins-Monro, Kiefer-Wolfowitz, SPSA) | continuous | gradient steps with noisy gradient estimates and decreasing step sizes |
| Gradient estimation (IPA, likelihood ratio) | continuous | derivative of the output estimated from the same run |
| Metamodels (response surfaces, kriging, Bayesian optimization) | continuous, expensive runs | fit a surrogate to simulated points and optimize the surrogate |
| Metaheuristics with simulation (OptQuest, racing of Section 11.3) | large discrete or mixed | search with simulated fitness, more replications for promising designs |
SPSA (simultaneous perturbation stochastic approximation) estimates the whole gradient with only two simulations per iteration, whatever the dimension $p$. It perturbs every component at once with a random $\pm 1$ vector $\Delta_k$:
$$ \hat g_{k,i} = \frac{F(x_k + c_k \Delta_k) - F(x_k - c_k \Delta_k)}{2,c_k,\Delta_{k,i}}, \qquad x_{k+1} = x_k - a_k,\hat g_k. $$
Good practice applies to all methods:
- Use CRN across candidates.
- Screen many designs cheaply, then select among the few survivors with more replications.
- Re-estimate the chosen design on fresh, independent replications, and report that estimate with a confidence interval.
Example: tuning an $(s, S)$ policy by simulation. The inventory item of Section 10.1 has Poisson weekly demand with mean 10, order cost $K = 64$, holding cost 1 and backorder cost 9. Its exact optimum is $(s, S) = (6, 40)$ with cost 35.02 per week, which allows grading the simulation optimizers. Each of 1,180 candidate policies is simulated for 300 weeks. Three methods are compared over 50 independent experiments: independent demand streams per policy, a common demand stream (CRN and SAA), and a two-stage procedure that re-tests the 10 best CRN candidates on 5,000 fresh weeks.
from functools import lru_cache
import numpy as np
# avg_cost(s, S): exact long-run cost, defined in the (s, S) Policies example (Section 10.1)
lam, K, h, b = 10.0, 64.0, 1.0, 9.0 # same item; avg_cost reads these globals
pol = np.array([(s, S) for s in range(0, 20) for S in range(s + 1, s + 60)]) # 1,180 policies
s_arr, S_arr = pol[:, 0], pol[:, 1]
def sim_cost(demand):
"""Average weekly cost of every policy at once. demand has shape (weeks,) for common
random numbers, or (weeks, policies) for independent streams."""
x, total = S_arr.astype(float), np.zeros(len(pol))
for t in range(demand.shape[0]):
order = x <= s_arr # review: order up to S if x <= s
y = np.where(order, S_arr, x)
x = y - demand[t] # demand is backordered if x < 0
total += K * order + h * np.maximum(x, 0) + b * np.maximum(-x, 0)
return total / demand.shape[0]
exact = lru_cache(None)(lambda s, S: avg_cost(int(s), int(S))[0])
OPT = exact(6, 40) # 35.02, the known optimum
rng = np.random.default_rng(0)
weeks = 300 # simulation budget per policy
res = {"independent streams": [], "common random numbers": [], "two-stage (CRN + re-test)": []}
for rep in range(50): # 50 independent optimization experiments
for name in res:
D = rng.poisson(10, (weeks, len(pol)) if name.startswith("indep") else weeks)
est = sim_cost(D)
if name.startswith("two"): # re-evaluate the 10 best on 5,000 fresh weeks
top = np.argsort(est)[:10]
fresh = sim_cost(rng.poisson(10, 5_000))[top]
k, est_k = top[np.argmin(fresh)], fresh.min()
else:
k, est_k = int(np.argmin(est)), est.min() # sample average approximation: take the best
res[name].append((exact(*pol[k]), est_k))
print("method mean true cost worst mean gap estimate - true")
for name, r in res.items():
true, est = np.array(r).T
print(f"{name:26s} {true.mean():10.2f} {true.max():8.2f} {true.mean() / OPT - 1:+6.1%} {np.mean(est - true):+10.2f}")
near = sum(exact(s, S) <= 1.02 * OPT for s in range(0, 14) for S in range(25, 56))
print(f"policies within 2% of the optimum {OPT:.2f}: {near}")
# method mean true cost worst mean gap estimate - true
# independent streams 35.26 36.01 +0.7% -1.32
# common random numbers 35.21 35.73 +0.5% -0.91
# two-stage (CRN + re-test) 35.07 35.20 +0.1% -0.07
# policies within 2% of the optimum 35.02: 62
All three methods find good policies, because the cost surface is flat near the optimum: 62 policies are within 2% of it. The differences show up in reliability and honesty. With independent streams, the chosen policy is on average 0.7% worse than optimal, 2.8% in the worst experiment, and its estimated cost is 1.32 too low. That gap is the optimizer's curse: the winner among 1,180 noisy estimates is mostly the luckiest one. CRN makes the comparison fairer, because all policies face the same demand, but the reported cost is still optimistic. The two-stage procedure is the most reliable (never worse than 0.5%), and its fresh estimate is almost unbiased ($-0.07$). The same pattern holds for staffing, buffer sizing and any other design chosen by simulation: screen with CRN, then confirm on new random numbers.
Key takeaways
- Simulation-based optimization minimizes an expectation that is observed only with noise. It needs methods that are robust to noise and economical with runs.
- Ranking and selection suits a few alternatives, SAA turns the problem into a deterministic one, stochastic approximation (SPSA) handles continuous variables, and metamodels or metaheuristics handle large spaces.
- The optimizer's curse biases the winner's estimate low. Re-evaluate the chosen design on fresh replications before reporting it.
- CRN, screening followed by selection, and flat optima make good designs easy to find, so focus on reliability and honest reporting.
Practice
- (Beginner) Ten staffing plans are each simulated with 5 replications, and the best looks 3% cheaper than the second best. Explain why the "best" may not be best, and why its estimated cost is probably too low.
- (Intermediate) Implement SPSA for the $(s, S)$ problem, treating $s$ and $S - s$ as continuous and rounding inside the simulator. Start at the heuristic $(13, 49)$ of Section 10.1, use CRN within each pair of evaluations, step sizes $a_k = a/(k+10)^{0.602}$ and perturbations $c_k = c/k^{0.101}$. Compare the policy after 200 iterations with the exact optimum.
- (Interview) "We must choose among 5 warehouse layouts. One replication takes an hour, and we have 100 hours of computing. How do you allocate the runs?" Hint: think about ranking and selection with an indifference zone, about sending more replications to the close contenders (OCBA), about CRN, and about what difference in performance actually matters to the business.
<a id="13-computational-tools-and-practice"></a>