ch13 · Monte Carlo

Builds on ch12.

The question

The sizing model has produced a host count. How sure are we?

ch12 took a stated workload, multiplied along three chains, took the largest of the three answers, and produced a number. Every step was arithmetic you could check by hand. The number is correct. Whether it is right is a different question. Every input to that chain was itself uncertain, and the chain has no way to say so.

This chapter builds the machinery for asking the second question. It assumes you can read code and do arithmetic. It assumes nothing about statistics.

The material

A single number is a bet you did not know you placed

Start with the model as ch12 left it.

Source — web_service-reference · every input on a slider

The first column is what the chain produced: one value per output, from one value per input. The second column is the same model, with each input allowed to be as uncertain as the person who wrote it down is.

Look at the host count. The point estimate is a real number, correctly computed. It sits inside a range that spans an order of magnitude, nearer the low end than the middle. Nothing went wrong. The calculation had no way to mention that its inputs were guesses, so it did not mention it.

Everything below is how that second column was produced.

Instead of one value, a bag of values

The idea is simple. If you do not know what the growth rate will be, do not give the model one growth rate. Give it a bag of plausible growth rates. Run the model once for every value in the bag. You get a bag of answers out, and the bag is the answer.

That is Monte Carlo. Everything else is bookkeeping: how to fill the bag, and how to read it.

Two words, and you will not need many more. The bag of plausible values for an input is its distribution. One value drawn from the bag is a sample. The bag of answers that comes out the other end is the output’s distribution. Reading one is a later section of this chapter.

Where the bag comes from: pick a percentile at random

Filling the bag takes one line of code, and the same line works for every distribution.

Every distribution can be described by a function that answers one question: what value sits at this percentile? Give it 0.5 and it hands back the middle value. Give it 0.9 and it hands back the value that nine tenths of the distribution is below. Call that function the distribution’s percentile function.

Now: pick a percentile uniformly at random, between zero and one, and ask the function what value sits there. Do it a hundred thousand times. You have sampled the distribution.

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)

generator.random(n) draws the percentiles. The percentile function turns them into values. That is the entire sampler, and it is why adding a distribution to this book is three lines rather than a new dependency.

The technique is called inverse transform sampling. Any distribution whose percentile function you can write down, you can sample. Problem 13.1 asks you to write one.

Here is the simplest:

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)

Percentile zero gives the minimum, percentile one gives the maximum, and everything in between is a straight line. You could have guessed that one. The next one needs some thought:

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 the mode, the area under the triangle grows as the square of the distance from the minimum, so inverting it gives a square root. That is the whole derivation. Do it on paper once. Every other distribution in this chapter is the same exercise with different algebra.

Which shape for which input

Choosing a distribution is a judgement, not a technical question. It is a claim about the world, and the first claim a reviewer should argue with.

Lognormal, for prices and growth and anything that compounds. Two properties make it the right shape for those. It cannot go negative, and neither can a price. And a product of several lognormals is another lognormal. A sizing chain is a product, so the uncertainty arriving at the end of one has roughly this shape whether or not anybody chose it.

The parameters here are two percentiles rather than the mean and standard deviation of a logarithm. Nobody has an intuition for the standard deviation of a logarithm. Everybody has one for “I would be surprised if it were under eleven or over nineteen”. That is a sentence a person can say about a price, and it is what the model file records.

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))

Triangular, for an expert’s guess. The least it could be, the most it could be, and the one they would bet on. Most sizing inputs arrive in this shape, because it is the shape of the answer to “what is it, roughly?”. Name its flaw every time you use it: it asserts that nothing outside the bounds can happen, and the bounds came out of somebody’s memory.

Uniform, when the bounds really are all you know. A price capped by a contract. A retention window somebody will pick from a range. Honest there, and dishonest as a default. It says the extremes are as likely as the middle, and almost nothing real is like that.

Normal, for measurement error. In this book it means one thing: the standard error beside a measured constant. A number was measured, the measurement wobbles, and it is as likely to wobble high as low. It is the wrong default for a price: a normal will happily go negative and a price will not.

Choosing badly is not a rounding error. It is a claim about what can happen, made in a model file that will outlive the meeting it came from. So the toolkit refuses an input that is sampled without naming its shape and saying why. Every input in this book records who claimed it, on what basis, and which of the four shapes it is:

Source — web_service-reference · every input on a slider

The tally at the bottom is the honest summary of any model. This many of the numbers are traceable. This many are somebody’s sales material. This many were decided in a room.

Running the bag through the model

Nothing changes.

The model is a graph of quantities, each computed from the ones before it. Evaluating it at a point walks the graph in order, doing arithmetic on numbers. Sampling it walks the same graph in the same order, doing the same arithmetic on arrays. a * b means the same thing whether a and b are two numbers or two hundred thousand. The evaluator in this book has one expression walker, and hands it a different set of functions for each pass.

So uncertainty propagates for free. An input with a distribution becomes an array. Everything downstream of it becomes an array. Everything else stays a single number and broadcasts.

Every node in that graph has a distribution once the model has been sampled, not just the ones at the end. The interactive version of this figure will show you any of them. That is the fastest way to find out where an interval got wide: you walk the chain until the histograms stop being narrow.

Reading the answer

Two words, and then we are done with vocabulary.

A percentile is the value a given fraction of the bag is below. The 5th percentile is the value only one sample in twenty came in under.

An interval is the gap between two of them. This book reports the gap between the 5th and the 95th percentile, calls it the 90% interval, and deliberately does not call it a confidence interval. That phrase means something precise to a statistician and something vaguer to everybody else. What is meant here is the plain reading: the model put nine tenths of its belief in this range.

Why two percentiles rather than the smallest and the largest answer in the bag: the ends are properties of how many answers you collected, not of the model. Collect ten times as many and the largest gets larger, every time, because the unlucky combinations had more chances to turn up. The middle settles, and ch14 measures how quickly.

Note where the red line sits relative to the middle of the distribution. For a chain of multiplications with skewed inputs, the answer you get from the average inputs is not the average answer, and it is not the middle one either. There is a theorem behind that. You do not need it. You need to have seen it happen once.

How often each ceiling is breached

The ceilings are what the machinery was built for, and a spreadsheet has no equivalent of this table.

Source — web_service-reference · every input on a slider

At the point estimate, every ceiling but one is fine. That is not surprising. ch12 sized the fleet from the point estimates, so of course it satisfies the ceilings it was sized against. The one it was not sized against is the one that is not fine. The last two columns are the same model asked a different question: across everything this model thinks could happen, how often is this limit breached?

That is a sizing answer. Not “you need this many hosts”, but “at this many hosts, this is how often the thing you were trying to avoid happens anyway”. Somebody can take responsibility for the second. Nobody can take responsibility for the first, because it does not say anything.

Once the question is in that form, it has a price. Here is the same model with a bigger fleet bought:

Source — web_service-sized_for_growth · every input on a slider

What the extra capital buys is the difference between two percentages. Whether it is worth it is not a modelling question, and ch21 is about how to put it to the person whose decision it is.

The seed

Every sampled result in this book records the seed its random numbers came from, and the sample count, and a hash of the sampler’s source. Run it again with those three and you get the same numbers to the last digit.

That is not fastidiousness. An unseeded simulation is a measurement nobody can repeat, and a figure nobody can repeat is a figure nobody can check. This book refuses that everywhere else, and has no reason to start allowing it here.

What this cannot tell you

Whether the model has the right shape. Everything above takes the structure as given and asks what the inputs are worth. If a cost line is missing, if a ceiling was never declared, or if two quantities were multiplied that should have been added, sampling will propagate the error beautifully and report a confident interval around the wrong answer. That is structural error. It is invisible to every technique in this chapter, and it is the subject of ch20 · The missing node.

Whether the shapes were chosen honestly. A triangular with generous bounds and a lognormal with tight ones will give different intervals for the same input, and nothing here can tell you which was right. The distribution is an assumption like any other, and this book makes you write it in a file with your name on it for that reason.

How the inputs were drawn together. The sampler above draws each input on its own. The intervals above were not produced that way. This model declares two pairs that move together: a host’s price and the network’s, quoted by the same supply chain; and the busy hour and what a request costs, because a busier service is a slower one per request. The evaluator applies the pairing after the draw. So every figure on this page already carries it, and this chapter has not said so until now. Drawing those pairs independently would make every interval here narrower, which is the direction that gets a plan approved. ch14 names the pairs, and measures what they were worth.

Whether a hundred thousand samples was enough. This chapter assumed it and did not establish it. ch14 has the argument, and the way to work it out for a model of your own.

How likely any of this is. The interval is a statement about the model’s declared inputs. It is not a forecast, it carries no track record, and its 95th percentile is not a promise. It is the best available summary of what you have written down. That is worth a great deal more than a single number, and a great deal less than knowledge.

Problems

Four. The first three are in tests/monte_carlo/ and are graded against definitions the tests compute for themselves. The fourth has no test and no known answer.

13.1 — Add a distribution. Implement the percentile function for a shape this book does not have: a quantity known to within a factor, whose density is proportional to one over the value. Derive it as the chapter derived the triangular’s. The test grades it against the density itself, integrated at test time, so there is nothing to look up.

def log_uniform_ppf(u: np.ndarray, minimum: float, maximum: float) -> np.ndarray:
    """Problem 13.1 — the percentile function of a shape this book does not have.

    A quantity you know to within a factor and no better: as likely to sit anywhere in its range
    on a multiplicative scale, so that doubling is as plausible as halving wherever you start. Its
    density is proportional to one over the value, between ``minimum`` and ``maximum`` and zero
    outside them::

        density(x) = 1 / (x * ln(maximum / minimum))     for minimum <= x <= maximum

    That is the whole specification. It is the shape for a ballpark: how many distinct label
    values a service will turn out to carry, the size of a table nobody has counted.

    Your job is the percentile function: given percentiles in (0, 1), return the values that sit
    at them. Derive it the way the chapter derived the triangular's: the area under the density
    from the minimum up to a value, set equal to the percentile, then inverted. It comes out as
    one line.

    The test grades it against the density rather than against a formula. It integrates the
    density numerically up to each value you return and checks that the area is the percentile
    you were given. So there is nothing to look up, and a function that is monotonic and bounded
    but the wrong shape fails.
    """
    raise NotImplementedError("problem 13.1")

13.2 — Sample a model by hand. The test hands you two of the web service model’s inputs, each with the distribution the model file declares for it. Sample both yourself, without the toolkit’s sampler, work the cost they feed through by hand, and reproduce the interval this book publishes for it to within sampling error. The point is to discover how small the machinery is.

def sample_two_inputs(declared: dict[str, dict], seed: int, samples: int) -> dict[str, np.ndarray]:
    """Problem 13.2 — sample a model by hand.

    ``declared`` holds two of the web service model's inputs, ``staff_fte`` and
    ``fully_loaded_salary``, each keyed to the distribution the model declares for it, as the
    file writes it. One is a triangular::

        {"triangular": {"minimum": ..., "likely": ..., "maximum": ...}}

    and one is a lognormal::

        {"lognormal": {"p10": ..., "p90": ...}}

    The test reads both out of ``models/web_service/model.yaml`` and hands them to you, so a
    change to the model changes what you are given instead of going stale in a number typed
    here. Sample what you are handed.

    Return a dictionary with three keys. ``staff_fte`` and ``fully_loaded_salary`` each hold
    ``samples`` draws from their distribution, seeded with ``seed``; ``annual_staff_cost`` holds
    what those draws imply, worked draw by draw through the model's own formula for it, which
    you read out of the same file.

    Do it the way the chapter did: pick the percentiles uniformly at random, then ask each
    distribution what value sits at each one. Use ``numpy`` and ``sizing.normal`` if you like.
    Do **not** use ``sizing.mc`` or ``sizing.evaluate``: the point of this one is to find out how
    little machinery there is between a declared distribution and a bag of numbers, and between
    the bag and an interval the book publishes.

    The test compares your inputs' percentiles against what the declared distributions say, and
    your cost's interval against the one the book publishes for the reference scenario, each to
    within the sampling error the count allows — which means it also checks that you understood
    what "to within sampling error" has to mean here. These two inputs are declared independent
    of everything else in the model, which is what makes the interval reachable by hand; ch14 is
    about the ones that are not.
    """
    raise NotImplementedError("problem 13.2")

13.3 — Where the point estimate sits. For each output of the web service model, find the share of the sampled answers that fall below the point estimate. For a cost it is near a half. For the recommended host count it is not, and the chain says why. Name the step that moved it.

def where_the_point_sits(
    samples: dict[str, np.ndarray], point: dict[str, float]
) -> dict[str, float]:
    """Problem 13.3 — where the point estimate sits among the sampled answers.

    ``samples`` holds every output's draws and ``point`` the same outputs' point estimates: the
    value the model computes from every input's middle. Return, for every output in ``point``,
    the fraction of that output's draws that fall below its point estimate.

    For a quantity that is one input passed straight through, the answer is a half. For a product
    of the model's shapes it is close to a half. For the web service's recommended host count it
    is not, by more than sampling error, and the chain says why: at least two of its steps are
    not multiplications, a maximum over three chains and a rounding up inside each, and a point
    estimate walks through them as if they were.

    Then finish the last line of this docstring with one sentence naming the step that moves the
    point furthest from the middle, and in which direction. The test grades the fractions, not
    the sentence; whoever reads your model will grade the sentence.

    What moved it:
    """
    raise NotImplementedError("problem 13.3")

13.4 — Defend a distribution. No test: the invoices are yours, and so is the shape you would defend.

Take a price you pay, find two years of invoices for it, and decide which of the four shapes in this chapter you would use and why. Then check what the last two years would have looked like under your choice. If the answer embarrasses you, that is the exercise working.

A good answer names the shape, the reason it and not another, and the two ends of the band, as a sentence somebody could disagree with. What would show it wrong is the invoices: more than a couple of the twenty-four outside the band you declared means the band was too narrow, and none anywhere near its ends means it was too wide to be a claim at all.

Where to go next

Metropolis and Ulam’s original paper [Metropolis & Ulam (1949)] is seven pages, is readable without any statistics, and is a useful corrective to the idea that this is a modern technique.

numpy.random’s documentation on generators and seeding is worth twenty minutes, particularly the part about why default_rng exists and what it replaced.

ch14 picks up the two things this chapter used without establishing: the correlations the intervals above already carry, and the sample count.