Appendix B · The Monte Carlo module, read end to end

ch13 and ch14 quote pieces of this module where they need them. This page is all of it, in the order it is written, for a reader who wants to see that there is nothing else in it.

There is no simulation framework underneath this and no statistics package beside it. numpy for arrays, one rational approximation for the inverse normal, and nothing else. The reason is not minimalism: a reader who cannot see the sampler cannot check the interval, and an interval nobody can check is decoration.

The one idea

## One idea

Every distribution here is sampled the same way, and it is the only sampling idea in the book:

    draw a percentile uniformly at random, and ask the distribution what value sits at it.

That is *inverse transform sampling*. It is why each distribution below needs exactly one
function — its percentile function, ``ppf`` — and why adding a distribution to this book is
three lines rather than a new dependency. Problem 13.1 asks the reader to add one.

Each distribution needs exactly one function: the value at a given percentile. Sampling it is drawing percentiles at random and looking the values up. Everything below is arrangement around that.

The seed

def rng(seed: int) -> np.random.Generator:
    """The generator, seeded.

    One per run, created here so that every result in the book can record the seed that produced
    it and be reproduced exactly. An unseeded Monte Carlo is a measurement nobody can repeat,
    which is the same thing this repository refuses everywhere else.

    PCG64 rather than the legacy Mersenne Twister: numpy's modern generator is what
    ``default_rng`` gives you, and pinning it by name means a numpy upgrade cannot silently
    change which stream a stamped seed refers to.
    """
    return np.random.Generator(np.random.PCG64(seed))

Every stamped result computed from a model records the seed that produced it, and bench/stamp.py refuses one that does not. An unseeded run is a measurement nobody can repeat, and this repository refuses those everywhere else. Where an experiment uses many seeds, as the convergence table on this page does with one per replicate, the result records the seed they are all derived from, and the rule that derives them.

The four percentile functions

Uniform

def uniform_ppf(u: np.ndarray, minimum: float, maximum: float) -> np.ndarray:
    """Every value between two bounds, equally likely.

    Honest when the bounds are genuinely all you know — a contract that caps a price, a retention
    window somebody will choose from a range. Dishonest as a default, because it says the bounds
    are as likely as the middle, and almost nothing real is like that.
    """
    return minimum + u * (maximum - minimum)

Triangular

def triangular_ppf(u: np.ndarray, minimum: float, likely: float, maximum: float) -> np.ndarray:
    """An expert's guess: the least it could be, the most, and the one they would bet on.

    The shape most sizing inputs arrive in, because it is the shape of the answer to "what is it,
    roughly?". Its flaw is worth stating every time it is used: it asserts that nothing outside
    the bounds can happen, and the bounds came from somebody's memory.

    The two branches meet at the mode. Below it the area grows as the square of the distance from
    the minimum, which inverts to a square root — the whole derivation, and problem 13.1 asks for
    it again for a distribution that is not here.
    """
    if not minimum <= likely <= maximum:
        raise ValueError(
            f"triangular needs min <= likely <= max, got {minimum}, {likely}, {maximum}"
        )
    if maximum == minimum:
        return np.full_like(u, minimum)
    width = maximum - minimum
    at_mode = (likely - minimum) / width
    below = minimum + np.sqrt(u * width * (likely - minimum))
    above = maximum - np.sqrt((1.0 - u) * width * (maximum - likely))
    return np.where(u < at_mode, below, above)

Two branches meeting at the mode. Below it the area under the triangle grows as the square of the distance from the minimum, so inverting it is a square root; above it, the same thing from the other end. That derivation is all there is to adding a distribution to this book, and problem 13.1 asks for it again for a shape that is not here.

Lognormal

def lognormal_ppf(u: np.ndarray, p10: float, p90: float) -> np.ndarray:
    """Multiplicative uncertainty: the shape of prices, growth rates and anything compounding.

    Parameterised by two percentiles rather than by the mean and standard deviation of the
    logarithm, because nobody has an intuition for the second and everybody has one for the first.
    "I would be surprised if it were under 11 or over 19" is a sentence a person can say about a
    price, and it is exactly ``p10=11, p90=19``.

    Two properties earn it its place. It cannot go negative, and neither can a price. And a
    product of several of them is another one, which is what a chain of multiplications in a
    sizing model *is* — so the uncertainty that arrives at the end of the chain has this shape
    whether or not anybody chose it.
    """
    if not 0 < p10 < p90:
        raise ValueError(f"lognormal needs 0 < p10 < p90, got p10={p10}, p90={p90}")
    # Solve for the two parameters of the underlying normal from the two stated percentiles.
    log_median = (np.log(p10) + np.log(p90)) / 2.0
    log_spread = (np.log(p90) - np.log(p10)) / (2.0 * Z90)
    return np.exp(log_median + log_spread * normal_ppf(u))

Parameterised by two percentiles, not by the mean and standard deviation of the logarithm. Nobody has an intuition for the mean of a logarithm. Everybody has one for percentiles: I would be surprised if it were under this, or over that is a sentence a person can say about a price.

Normal

def normal_ppf_scaled(u: np.ndarray, mean: float, sd: float) -> np.ndarray:
    """Symmetric error around a central value.

    In this book it means one thing: the measurement uncertainty of a ``measured`` node (ch03).
    A constant was measured, the measurement has a standard error, and that error is as likely to
    be high as low. It is the wrong default for a price — see :func:`lognormal_ppf` — and the
    wrong shape for anything that cannot go negative, which it happily will.
    """
    if sd < 0:
        raise ValueError(f"normal needs sd >= 0, got {sd}")
    return mean + sd * normal_ppf(u)


#: The distributions a model file may declare, by the key it declares them under. Adding one is a
#: percentile function and a line here; there is no registration machinery and no plugin system,
#: because four shapes cover every input in both reference models and a fifth should have to
#: argue for itself.

One job in this book: the measurement error of a measured node. It is the wrong default for a price, because it will happily go negative.

The registry

SHAPES: dict[str, Callable[..., np.ndarray]] = {
    "uniform": uniform_ppf,
    "triangular": triangular_ppf,
    "lognormal": lognormal_ppf,
    "normal": normal_ppf_scaled,
}

No registration machinery and no plugin system. Four shapes cover every input in both reference models, and a fifth should have to argue for itself (Appendix C).

Drawing

def sample(spec: dict, n: int, generator: np.random.Generator) -> np.ndarray:
    """Draw ``n`` values from a declared distribution.

    Two lines, and they are the two lines of the whole chapter: draw percentiles, look up values.
    """
    shape, parameters = one_shape(spec)
    return SHAPES[shape](generator.random(n), **parameters)

Two lines, and they are the two lines of the entire subject.

def one_shape(spec: dict) -> tuple[str, dict]:
    """Pull the single distribution out of a declaration, refusing ambiguity.

    A node that declares two shapes is not a node with a shape to be guessed at; it is a model
    file somebody edited without deleting the old line, and picking either one for them hides it.
    """
    declared = [key for key in spec if key in SHAPES]
    if len(declared) != 1:
        known = ", ".join(sorted(SHAPES))
        raise ValueError(
            f"a distribution declares exactly one shape, got {sorted(spec)!r}. Known shapes: {known}"
        )
    shape = declared[0]
    return shape, dict(spec[shape])

Inputs that move together

def correlation_matrix(names: list[str], pairs: list[dict]) -> np.ndarray:
    """A full correlation matrix from the pairs a model bothered to declare.

    Everything not named is left at zero, which is an assumption and not a fact — ch14 is mostly
    about how much that assumption costs. Stating it here rather than hiding it in a default is
    the point of building the matrix explicitly.
    """
    index = {name: i for i, name in enumerate(names)}
    matrix = np.eye(len(names))
    for pair in pairs:
        a, b, rho = pair["a"], pair["b"], float(pair["rho"])
        if a not in index or b not in index:
            missing = a if a not in index else b
            raise ValueError(f"correlation names {missing!r}, which is not a sampled input")
        if not -1.0 <= rho <= 1.0:
            raise ValueError(f"correlation between {a} and {b} is {rho}, outside [-1, 1]")
        matrix[index[a], index[b]] = matrix[index[b], index[a]] = rho
    return matrix

def rank_to_score_correlation(rank_rho: np.ndarray) -> np.ndarray:
    """Convert the correlation a modeller declares into the one the machinery needs.

    A subtlety that is easy to skip and then be quietly wrong about. :func:`correlate` works by
    correlating *normal scores* and then reordering, and the rank correlation that comes out is
    not the number that went in — it is attenuated, by a known amount:

        rank correlation  =  (6 / pi) * arcsin(score correlation / 2)

    Declare 0.8 and, uncorrected, you get 0.785. Small, consistent, and exactly the kind of error
    that survives review forever because nobody expects the number they typed to be a different
    number. So the relation is inverted here: a model file's ``rho`` means the rank correlation
    the modeller wants, and this works out what to ask the scores for.
    """
    return 2.0 * np.sin(np.pi * rank_rho / 6.0)

Inducing a correlation on ranks and then reading it back on values does not return the number you asked for. It comes back attenuated, by an amount that depends only on the coefficient: ask for a strong correlation, measure the result, and it is visibly weaker. This is the correction that is easy to leave out and hard to find afterwards, and the fix is one line of trigonometry before the sort rather than an apology in the documentation after it.

def correlate(
    columns: np.ndarray, target: np.ndarray, generator: np.random.Generator
) -> np.ndarray:
    """Induce a rank correlation between columns without touching their distributions.

    Iman and Conover's method @imanconover1982. It is worth understanding rather than importing,
    because what it does *not* do is the reason it is the right tool: it only ever **reorders**
    each column. Every value that was going to be in a column is still in it, so each input keeps
    exactly the distribution the modeller chose, and the only thing that changes is which draws
    line up with which. Correlating the values instead — the obvious first attempt — quietly
    changes the marginals, and then the model is answering a question nobody asked.

    Problem 14.2 checks both halves of that claim: the interval widens, and the marginals do not
    move.

    ``columns`` is (samples, inputs); ``target`` is the square matrix from
    :func:`correlation_matrix`, holding rank correlations.
    """
    n, k = columns.shape
    if target.shape != (k, k):
        raise ValueError(f"correlation matrix is {target.shape}, expected {(k, k)}")
    if k == 0 or np.allclose(target, np.eye(k)):
        return columns

    # A reference set with the right shape and no correlation: the normal scores, independently
    # shuffled per column, from the run's own generator so the whole thing reproduces from one
    # stamped seed. Working in scores rather than in the data is what makes the method
    # indifferent to what the marginals actually are.
    scores = normal_ppf(np.arange(1, n + 1) / (n + 1))
    reference = np.column_stack([generator.permutation(scores) for _ in range(k)])

    wanted_scores = rank_to_score_correlation(target)
    np.fill_diagonal(wanted_scores, 1.0)
    try:
        wanted = np.linalg.cholesky(wanted_scores)
    except np.linalg.LinAlgError as exc:
        raise ValueError(
            "the declared correlations are not mutually consistent — no set of inputs can have "
            "all of them at once. Check for a triangle of strong correlations that disagree."
        ) from exc
    have = np.linalg.cholesky(np.corrcoef(reference, rowvar=False))
    shaped = reference @ np.linalg.solve(have, wanted).T

    # Reorder each column to follow the shaped scores' ranking. `argsort` twice gives the rank of
    # every element; sorting the column and indexing by rank puts the largest value where the
    # largest score is, and so on down.
    out = np.empty_like(columns)
    for j in range(k):
        ranks = np.argsort(np.argsort(shaped[:, j]))
        out[:, j] = np.sort(columns[:, j])[ranks]
    return out

Iman–Conover [Iman & Conover (1982)], which is short enough to read: build a reference sample with the correlation you want, rank it, and shuffle each input column into the same rank order. Every column keeps its own distribution exactly, because every value that was drawn is still there, and only the pairing between columns changes. That is why it works on all four shapes without knowing anything about them.

Reading the answer

PERCENTILES = (5, 25, 50, 75, 95)


def summarise(x: np.ndarray) -> dict:
    """What a bag of values is reported as.

    Interpolated percentiles, unlike the book's measurement chapters, and for a reason worth
    naming: these are 100,000 draws from a continuous distribution, not eleven timings of a real
    program. There is no "sample that actually happened" to prefer here — every one of them is
    something the model made up.
    """
    x = np.asarray(x, dtype=float)
    values = np.percentile(x, PERCENTILES)
    return {
        "min": float(np.min(x)),
        **{f"p{p}": float(v) for p, v in zip(PERCENTILES, values, strict=True)},
        "max": float(np.max(x)),
        "mean": float(np.mean(x)),
        "sd": float(np.std(x, ddof=1)) if x.size > 1 else 0.0,
    }

def interval(x: np.ndarray, lo: float = 5, hi: float = 95) -> tuple[float, float]:
    """The range the model puts ``hi - lo`` per cent of its belief in.

    Not a confidence interval and not a guarantee. It is a statement about this model's inputs,
    and it is exactly as good as they are — which is what ch14's closing section is about.
    """
    low, high = np.percentile(np.asarray(x, dtype=float), [lo, hi])
    return float(low), float(high)


def half_width(x: np.ndarray, lo: float = 5, hi: float = 95) -> float:
    """Half the width of an interval — the number that falls as one over the square root of n."""
    low, high = interval(x, lo, hi)
    return (high - low) / 2.0

def samples_needed(observed_half_width: float, at_n: int, target_half_width: float) -> int:
    """How many draws to get the interval down to a width you would report.

    Straight from the square-root law: the width falls as one over the square root of the sample
    count, so to halve it you need four times as many. ch14 shows the law rather than asserting
    it, and this is the arithmetic that follows once you believe it.

    This is a statement about *sampling noise only* — how much the answer wobbles because it was
    made of a finite number of draws. It says nothing about whether the model is right, and the
    distinction is the most important one in the chapter.
    """
    if target_half_width <= 0:
        raise ValueError("target half-width must be positive")
    ratio = observed_half_width / target_half_width
    return int(np.ceil(at_n * ratio * ratio))

The arithmetic of “enough”, both ways round: the width you have at the draws you took, and the draws you need for the width you want. People usually settle this with a habit. The arithmetic is unforgiving — a factor of ten less wobble costs a hundred times the samples.

Source — convergence-tco · sizing.mc

Two columns, two behaviours, and confusing them is the most common misunderstanding in the subject. The interval settles — it is a property of how uncertain the model’s inputs are, and more samples converge on it rather than shrinking it. The run-to-run spread falls, at close to the square root of ten per decade, because that is a property of how hard you looked.

Shipping a distribution

#: Above this ratio between the 95th percentile and the 5th, equal-width bins stop describing the
#: quantity: most of the draws fall in the first bin and the rest of the picture is empty. A queue
#: near saturation does this. Deliberately measured across the interval rather than across the
#: extremes, because one stray draw from a long tail should not change how everything is binned.
LOG_BINNING_SPAN = 100.0


def histogram(x: np.ndarray, bins: int = 64) -> dict:
    """A node's distribution, small enough to ship to a browser for every node in the graph.

    Counts and edges rather than the draws themselves: forty numbers a node instead of a hundred
    thousand, which is what makes it affordable to let the reader click any node in the DAG and
    watch a narrow input distribution turn into a wide output one further down the chain.

    The bins are equal in width, unless the draws span more than :data:`LOG_BINNING_SPAN`, in
    which case they are equal in *ratio* and ``spacing`` says so. Whatever draws the histogram
    needs to know which it has: the same counts mean different things on the two axes.
    """
    x = np.asarray(x, dtype=float)
    low, high = float(x.min()), float(x.max())
    p5, p95 = (float(v) for v in np.percentile(x, [5, 95]))
    if low > 0 and p5 > 0 and p95 / p5 >= LOG_BINNING_SPAN:
        counts, edges = np.histogram(x, bins=np.geomspace(low, high, bins + 1))
        return {"counts": counts.tolist(), "edges": edges.tolist(), "spacing": "log"}
    counts, edges = np.histogram(x, bins=bins)
    return {"counts": counts.tolist(), "edges": edges.tolist(), "spacing": "linear"}

Counts and edges rather than the draws: a few dozen numbers a node instead of a hundred thousand. That is cheap enough for the interactive page to let a reader click any node and see its distribution, not only the outputs. Watching a narrow input turn into a wide output three steps down the chain is the fastest way to understand a sizing model, and it costs almost nothing to ship.

The bins are equal in width unless the quantity spans orders of magnitude, in which case they are equal in ratio and the payload says so. A queue near saturation does this: half the draws land in the first equal-width bin and the picture becomes a spike beside an empty page. The figure reads that flag and labels its axis accordingly; ch05’s concurrency figure is one that does.

What is deliberately absent

## What is deliberately absent

No variance reduction, no quasi-random sequences, no importance sampling. Every one of them
narrows an interval for the same number of draws, and every one of them also makes the interval
harder to explain to the person who has to sign for the money. This book's bottleneck is never
compute; it is whether the reader believes the answer.

No fitted distributions either. Nothing here reads data and tells you which shape it is. Choosing
a shape is an editorial act with provenance attached (ch03), and a function that guesses it for
you produces a model whose central assumption nobody ever wrote down.

Variance reduction, quasi-random sequences and importance sampling all narrow an interval for the same number of draws, and all of them make the interval harder to explain to the person who has to sign for the money. This book’s bottleneck was never compute.

Fitting is absent for a different reason. A function that reads your data and tells you which shape it is produces a model whose central assumption nobody ever wrote down. Choosing a shape is an editorial act with provenance attached (ch03), and it belongs in the model file where a reviewer can argue with it.

The inverse normal

Two of the four shapes need the inverse normal CDF, and there is no closed form. sizing/normal.py is Acklam’s rational approximation [Acklam (2003)], implemented here and checked against Python’s own statistics.NormalDist across the range:

def normal_ppf(u: np.ndarray | float) -> np.ndarray:
    """The value below which a fraction ``u`` of a standard normal distribution lies.

    ``normal_ppf(0.5)`` is 0, ``normal_ppf(0.9)`` is about 1.28, and ``normal_ppf(0.975)`` is the
    1.96 that turns up in every textbook margin of error.

    Raises on 0 and 1 rather than returning an infinity. Both are real bugs when they happen — a
    uniform draw is never exactly 0 or 1 under numpy's generator, so an endpoint here means a
    percentile was computed rather than drawn, and silently returning an infinity would put it in
    a sum and turn a whole model's output into ``nan`` several steps later.
    """
    u = np.asarray(u, dtype=float)
    if np.any((u <= 0.0) | (u >= 1.0)):
        raise ValueError("normal_ppf needs percentiles strictly between 0 and 1")

    out = np.empty_like(u)

    lower = u < _TAIL
    upper = u > 1.0 - _TAIL
    central = ~(lower | upper)

    # The middle 95%, as a rational function of the distance from the median.
    q = u[central] - 0.5
    r = q * q
    out[central] = _poly(_CENTRAL_NUM, r) * q / (_poly(_CENTRAL_DEN, r) * r + 1.0)

    # Both tails, in the variable that makes them well behaved. The upper tail is the lower one
    # reflected, which is worth doing explicitly: it is the only reason the function is accurate
    # at 0.999 as well as at 0.001, and a reader checking the code should be able to see that.
    for mask, tail_u, sign in ((lower, u[lower], 1.0), (upper, 1.0 - u[upper], -1.0)):
        q = np.sqrt(-2.0 * np.log(tail_u))
        out[mask] = sign * _poly(_TAIL_NUM, q) / (_poly(_TAIL_DEN, q) * q + 1.0)

    return out

Running it

python3 -m pytest tests/test_mc.py           # the sampler, against properties it must have
python3 -m bench.run_uncertainty             # the correlation and convergence experiments