ch14 · Correlation and convergence

Builds on ch13.

The question

ch13 produced an interval, and it rested on two things nobody checked: that every input moves on its own, and that a hundred thousand samples was enough to settle the answer.

Both are testable, and this chapter tests both. More samples do not make an interval narrower, which is the opposite of what most people expect.

The material

Inputs that move together

The web service model prices a host and the network port it plugs into. Ask anybody who has bought either: a year when hosts are scarce is usually a year when optics are, because they come through the same supply chain and are quoted in the same quarter.

Drawing them independently is not a neutral choice. It is the claim that one can save you from the other: that a bad host quarter will, on average, be offset by a good network quarter. If that is false, the model is reporting a narrower interval than the evidence supports.

Narrower is the direction that gets a plan approved.

Two inputs that tend to move together are correlated, and how strongly they do is their correlation. A model file can declare it, along with why:

Source — web_service-reference · every input on a slider

The last column is required. A correlation coefficient with no reason attached is a number somebody will copy into the next model without knowing what it was for.

Correlating ranks, not values

The obvious way is to correlate the values: nudge each host price up a little when the network price is up. Do that and you have changed the host price distribution: the thing you carefully chose in ch13, with its own percentiles and its own shape. You set out to encode one belief and quietly overwrote another.

The method this book uses only ever reorders. Every value that was going to be in a column is still in it, in the same quantity. All that changes is which draws line up with which. So each input keeps exactly the distribution the modeller chose. The correlation is expressed in the pairing, not in the numbers.

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

One subtlety in there is easy to skip, and then to be quietly wrong about. The method works by correlating normal scores, and the rank correlation that comes out is weaker than the one that went in, by a known amount. Apply no correction, and every declared correlation lands slightly weaker than it was written. That error is small and consistent, and it is the kind that survives review forever, because nobody expects the number they typed to come back as a different number. So the relation is inverted before use:

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)

Problem 14.2 checks both halves of the claim: that the correlation comes out where it was asked for, and that each input’s own distribution did not move.

What the correlations bought

Both reference models, sampled twice: once with their declared correlations, once with the correlations deleted.

Source — correlation-effect · sizing.mc — Iman-Conover rank correlation

Every row is positive. Assuming independence made every interval narrower, and in the observability model’s ingest chain it made it substantially narrower. That chain has several inputs feeding off the same growth, and pretending they are strangers lets them cancel each other out.

The web service’s rows are modest, and they are in the table deliberately: its two declared pairs are weak ones, and one of them does not reach the five-year total at all. Seeing a correlation that barely matters beside one that does is the fastest way to stop treating the subject as magic.

The general rule: correlation between inputs that push the same way widens the interval. It is not a correction. It is not a refinement. It does not make the model more precise. It removes an assumption that was making the model look better than it was.

How many samples is enough

Now the second thing ch13 assumed. The obvious experiment is to run the model at rising sample counts and watch the interval narrow.

That experiment does not work.

The interval does not narrow. A 90% interval is a property of the distribution the model describes, which is to say of how uncertain the model’s inputs are. More samples do not make that smaller. They converge on it: the answer settles towards the interval the inputs imply. Run the web service model with ten thousand draws and with a million, and the interval is the same width. It was never a function of how hard you looked.

What more samples buy is knowing where that interval is. Two runs of the same model with different random seeds give slightly different answers, and the gap between them shrinks as you draw more. So the experiment has to be run many times over at each sample count, and what gets recorded is the spread between those runs:

Source — convergence-tco · sizing.mc

Two columns, two behaviours. The first settles. The second falls, by about the square root of ten per decade across the range. That is the law measured rather than asserted, and measured with noise, which is the next paragraph.

Look at the individual ratios before you believe the summary, because they wander. Each spread in that column is itself estimated, from a limited number of independent runs. An estimate of a spread is noisy in the way everything else in this chapter is noisy. Measuring how uncertain something is turns out to be an uncertain measurement, and a figure demonstrating that law had better not be the one place in the book that forgets it. The overall rate across the range is far steadier than any single step. That is why it is the number on the last row.

The smallest row is excluded from the law and kept in the table. At that count a 95th percentile is one of the largest handful of draws there were, bounded by the sample itself, and nowhere near the regime the square-root law describes. It stays because it shows the other column at its clearest: that is the one row where the interval has visibly not settled.

To halve the wobble, take four times as many samples. To get it down by a factor of ten, take a hundred times as many. That is brutal, and it is why nobody buys precision this way past a point.

The law also gives a definition of “enough” that is a calculation rather than a habit:

Enough samples is when the answer stops moving between runs at the precision you are going to report it to.

If you are going to write the five-year total to the nearest hundred thousand, you need the run-to-run spread below that, and the table says which sample count gets you there. If you are going to write it to the nearest million, which is the more defensible choice given everything else in this book, you needed far fewer samples than you took.

sizing.mc has the arithmetic both ways round:

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


#: 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

What settling is not

A model whose answer has stopped moving between runs has settled its arithmetic. That is all it has settled.

You can run a model at a million samples, watch the interval stabilise to four significant figures, and present it with complete confidence, while a whole cost line is missing from the model. The sampling converged beautifully on the wrong number. Convergence is a statement about the calculation, never about the thing being calculated. ch20’s problem 20.3 hands you a model in that state, and an invoice it cannot reach.

What this cannot tell you

Which correlations exist. Everything above takes the declared pairs as given. Nothing here discovers a correlation, and nothing here warns you about one you failed to declare. From inside the model, an undeclared correlation is indistinguishable from a correlation of zero. That is the same failure as a missing node, and it belongs to ch20.

Whether the coefficient is right. A rank correlation declared as moderate rather than strong is a guess with the same standing as any other assumption in the model. The because field beside it in the model file is the only thing standing behind it.

Anything about correlations that are not monotonic. Rank correlation describes two quantities that tend to move in the same direction. Two that move together up to a point and then diverge, which is what a ceiling does to everything downstream of it, are not described by a single coefficient at all. This book does not pretend otherwise.

That more samples are ever the answer to a wide interval. They are not. More samples tell you where the interval is, not how wide it is. A wide interval means the inputs are uncertain. The only things that narrow it are measuring something or deciding something: ch19 · Which input to go and measure.

Problems

Four. The first two are graded, in tests/correlation_and_convergence/. The last two are not, and say why.

14.1 — Show the square-root law. Run one output of the web service model at several sample counts, with several independent seeds at each, and assert that the run-to-run spread falls as one over the square root of the count. The test hands you the run; the seeds, the percentile of each run and the spread between them are yours. The tolerance is itself a sampling question: the test derives it from the number of runs, and it runs as many as the chapter’s own experiment did, because a dozen is not enough to assert a ratio of two spreads.

def spread_at(run: Callable[[int, int], np.ndarray], samples: int, replicates: int) -> float:
    """Problem 14.1 — show the square-root law yourself.

    ``run(samples, seed)`` is one run of the web service model: that many draws of its five-year
    total, from that seed. The test builds it, so you do not have to touch the sampler; the seeds
    are yours to choose.

    Return the *run-to-run spread* of the 95th percentile of the total: call ``run``
    ``replicates`` times at ``samples`` draws, each with a different seed, take the p95 of each
    run, and return how far those p95 values spread from run to run, measured as their standard
    deviation, since that is the quantity the law is about.

    That is the quantity one-over-root-n governs. The interval itself is not — ch14 is mostly
    about the difference, and this problem is where you convince yourself.

    The test calls this at several sample counts and asserts that the value you return falls as
    one over the square root of the count, to a tolerance it computes rather than one that is
    written down.
    """
    raise NotImplementedError("problem 14.1")

14.2 — Correlate two inputs, and change neither. The page shows one entry for a model’s correlations block, with the coefficient and the reason left empty. Fill them in. The test adds it to a copy of the web service model with every other correlation taken out and asserts two things: that the interval on a shared output widens, and that neither input’s own distribution moves. The second is the property that makes the method trustworthy, and it is one line to check.

# Problem 14.2 - one entry for a model's `correlations:` block, in the file's own form.
#
# The test starts from the web service model with every correlation taken out, adds this entry,
# and checks three things: the interval on capital widens, neither price's own distribution
# moves, and the rank correlation that comes out is the one declared here. Two fields are empty.
# Declare a rank correlation of at least a half, since the section says why a weak one barely
# shows, and say why the two move together: the chapter calls the reason the required column.
- a: host_price
  b: network_price_per_host
  rho:
  because:

14.3 — Break the convergence experiment. No test. The experiment in this chapter uses a different seed for every replicate. Change it so that every replicate at a given sample count shares one seed, re-run it, and explain what the figure now shows and why it is worthless. Then say what else in this repository would have to be wrong for that mistake to survive review.

14.4 — Which of your inputs move together. No test: nothing here can see which of your quantities move together.

Nothing in this book discovers a correlation; they are all declared. Go through your own inputs in pairs and find the ones that are not independent: the growth rate and the peak ratio, the price and the quantity, the compression ratio and the kind of data.

For each pair, say which direction and roughly how strongly. Then say what it does to your answer. Correlated inputs moving the same way widen the result, and treating them as independent is the commonest way a model quietly reports less doubt than it has.

A good answer names at least one pair and says whether ignoring it makes your interval too narrow or too wide. If you find no pairs at all in a chain of six quantities about one system, look again. Independence is a strong claim, and it is rarely true.

Where to go next

Iman and Conover’s paper [Iman & Conover (1982)] is the method in this chapter, and is unusually readable for a statistics paper of its era. The section on what the method does not guarantee is the part to read twice.

ch19 is the question this chapter keeps deferring: given that the interval is wide, which single input should you go and measure?

ch20 is the failure that neither this chapter nor ch13 can see.