Software 16.5.1 (Random numbers and Monte Carlo). The Park–Miller generator, checked for mean and variance; exponential variates by the inverse transform; and a Monte Carlo estimate of \(\int _0^1e^{-v^2}\,dv\) with its \(95\%\) half-width, against the exact value \(\frac {\sqrt \pi }2\operatorname {erf}(1)\):

>>> import numpy as np
>>> def lcg(seed, n, a=16807, m=2**31 - 1):     # Park and Miller's generator
...     x, out = seed, []
...     for _ in range(n):
...         x = (a * x) % m; out.append(x / m)
...     return np.array(out)
...
>>> u = lcg(12345, 100000)
>>> round(u.mean(), 4), round(u.var(), 4)   # uniform: 1/2 and 1/12 = 0.0833
(0.5, 0.0834)
>>> x = -np.log(1 - u) / 2.0         # inverse transform: exponential, rate 2
>>> round(x.mean(), 4), round(x.std(), 4)  # both 1/2
(0.5008, 0.502)
>>> f = lambda v: np.exp(-v**2)      # Monte Carlo: int_0^1 exp(-v^2) dv
>>> fx = f(u); est = fx.mean(); half = 1.96 * fx.std(ddof=1) / np.sqrt(len(u))
>>> round(est, 5), round(half, 5)          # estimate and 95% half-width
(0.74684, 0.00125)
>>> from scipy.special import erf
>>> round(np.sqrt(np.pi) / 2 * erf(1), 5)  # the exact value
0.74682

Reading the output. The \(100\,000\) uniforms have mean \(0.5000\) and variance \(0.0834\) (exact \(\frac 1{12} = 0.0833\)); the exponentials with rate \(2\) have mean and standard deviation \(0.50\), as they should. The integral is estimated as \(0.74684 \pm 0.00125\), and the exact \(0.74682\) lies well inside.

MATLAB/Octave: rand, -log(1-u)/2, mean and std.

Software 16.5.2 (A discrete-event simulation). An event-driven simulation of the bank of Example 14.2.2 with \(3\) tellers, the event list kept as a heap; \(30\) independent replications of \(2000\) hours, and a \(95\%\) confidence interval for the mean wait in minutes, using \(t_{0.975, 29} = 2.045\):

>>> import heapq
>>> def bank(c, lam, mu, T, rng):          # event-driven M/M/c: mean wait
...     events = [(rng.exponential(1 / lam), 'arrive')]
...     free, queue, waits = c, [], []
...     while events:
...         t, kind = heapq.heappop(events)
...         if t > T: break
...         if kind == 'arrive':
...             nxt = t + rng.exponential(1 / lam)
...             heapq.heappush(events, (nxt, 'arrive'))
...             if free:                   # a teller is free: serve at once
...                 free -= 1; waits.append(0.0)
...                 done = t + rng.exponential(1 / mu)
...                 heapq.heappush(events, (done, 'end'))
...             else: queue.append(t)
...         elif queue:                    # a service ends: take the next
...             waits.append(t - queue.pop(0))
...             heapq.heappush(events, (t + rng.exponential(1 / mu), 'end'))
...         else: free += 1
...     return np.mean(waits)
...
>>> rng = np.random.default_rng(7)
>>> reps = np.array([bank(3, 24, 10, 2000, rng) for _ in range(30)]) * 60
>>> m, h = reps.mean(), 2.045 * reps.std(ddof=1) / np.sqrt(len(reps))
>>> round(m, 3), round(h, 3)               # minutes; Erlang C gives 6.47
(6.417, 0.127)

Reading the output. The mean wait is \(6.42 \pm 0.13\) minutes, consistent with the Erlang C value \(6.47\). With the program verified, it could now be changed in ways no formula covers: tellers taking breaks, customers giving up after a while, or a separate counter for business accounts.

MATLAB/Octave: a priority queue of events, or the SimEvents toolbox.

Software 16.5.3 (Variance reduction). Estimating \(\mathbb E[e^U]\) from \(2n\) independent uniforms against \(n\) antithetic pairs (the same number of function evaluations), with the variances of the two estimates; then the difference in expected newsvendor profit between orders of \(195\) and \(210\), estimated with independent and with common demands:

>>> rng = np.random.default_rng(3)
>>> g = lambda v: np.exp(v)                # estimate E exp(U) = e - 1
>>> n = 10000; U = rng.uniform(size=n)
>>> plain = g(rng.uniform(size=2 * n))     # 2n independent draws
>>> anti = (g(U) + g(1 - U)) / 2           # n antithetic pairs, same effort
>>> np.e - 1, round(plain.mean(), 5), round(anti.mean(), 5)
(1.718281828459045, 1.72028, 1.7194)
>>> round(plain.var() / (2 * n), 9), round(anti.var() / n, 9)   # variances
(1.2057e-05, 4.01e-07)
>>> d1 = rng.normal(200, 40, 5000); d2 = rng.normal(200, 40, 5000)
>>> def prof(q, d):                        # newsvendor profit of Section 13
...     return 10 * np.minimum(d, q) + np.maximum(q - d, 0) - 6 * q
...
>>> indep = prof(195, d1) - prof(210, d2)  # compare two orders, independently
>>> crn = prof(195, d1) - prof(210, d1)    # ... and with the same demands
>>> round(indep.mean(), 2), round(indep.std() / np.sqrt(5000), 3)
(6.0, 4.448)
>>> round(crn.mean(), 2), round(crn.std() / np.sqrt(5000), 3)
(11.88, 0.905)

Reading the output. Both estimates of \(e - 1 = 1.71828\) are close, but the antithetic one has variance \(4.0 \times 10^{-7}\) against \(1.2 \times 10^{-5}\), a factor of \(30\). For the policy comparison, independent samples give a difference of \(6.0\) with standard error \(4.4\), which cannot even tell which order is better; common random numbers give \(11.9\) with standard error \(0.9\), clearly favouring the order of \(195\), the one near the critical fractile.

MATLAB/Octave: the same vector operations.