Primer: probability and statistics

Read first

• Nothing. If you can compute an average and read a short Python loop, you have everything this primer needs.

The question this primer answers: when a routing method uses random numbers, how do you describe its results honestly, and how do you tell whether one method really beats another?

Intuition

From Post 15 on, the methods in this series roll dice: random starting tours, random kicks, random decisions about whether to keep a worse tour. Run the same method twice and you can get two different answers. This primer is about living with that.

We use the same eight stops as everywhere else. Distances are in units of 100 m, rounded, and the optimal tour A-B-C-D-H-E-F-G-A costs 315.

A random tour

The simplest random experiment in routing: start at A and visit the other seven stops in a random order, every order equally likely. Python’s random.Random(1) is a pseudo-random generator started from the seed 1. It produces numbers that look random but are completely determined by the seed, so the three tours below come out the same on every machine. Step through them. As in other posts, the diagram measures unrounded kilometres, so its lengths can differ by a few tenths from the matrix totals in the captions.

The whole distribution

With eight stops we don’t have to guess what a random tour looks like: there are only 2,520 distinct tours, so we can price them all. The table counts how many fall in each band of 50. That table is the distribution of the length of a random tour.

The average of all 2,520 lengths is exactly 488, 55% above the optimum. The longest tour, A-D-F-B-E-C-G-H-A, costs 620. Only 6 tours are within 10% of the optimum. A random tour is a poor tour, which is why every method in this series does better than rolling dice.

Samples

On berlin52 there are about 7.76 × 10⁶⁵ tours, far too many to price. So we sample: draw some tours at random and describe those. Here are samples of growing size, all from seed 1:

The sample mean wanders around the true 488, and wanders less as the sample grows. The standard error, defined below, measures by how much. And 1,000 samples found the optimum zero times. With probability 1/2,520 per draw, the expected number of hits in 1,000 draws is 1,000/2,520 ≈ 0.4, so finding none is the most likely outcome.

The formal idea

Probability and uniform random numbers. A probability is a number between 0 and 1 saying how often an event happens in the long run. The building block of every randomized method here is a number u drawn uniformly from [0, 1): every interval of length p inside it has probability p of containing u. So “do X with probability p” is written as one line of code: draw u, and do X if u < p.

flowchart TD
  P["Event with probability p"] --> U["Draw u uniformly from [0, 1)"]
  U --> Q{"u < p?"}
  Q -- yes --> Y["Event happens"]
  Q -- no --> N["It doesn't"]

Post 18 uses exactly this to keep a worse tour with probability e^(−Δ/T). For a tour 35 longer at T = 20 that is e^(−35/20) ≈ 0.174.

Expectation. If an event has probability p and you try N times, independently, the expected number of times it happens is p·N. With p = 0.174 and N = 1,000 that is 174. The script below draws 1,000 uniform numbers and counts 167: close to 174, not equal. Expectation describes the long-run average, not any single run.

More generally, the expected value of a random number X with possible values x₁, x₂, … and probabilities p₁, p₂, … is

E[X] = Σ xᵢ · pᵢ

For a uniformly random tour of the eight stops, each of the 2,520 tours has probability 1/2,520, so E[X] is simply the average of all 2,520 lengths: 488.

Variance and standard deviation. The variance measures spread around the mean μ = E[X], and the standard deviation σ is its square root, in the same units as X:

Var(X) = E[(X − μ)²]        σ = √Var(X)

For random tours, σ = 49.0. A sample of N values has a sample mean x̄ and a sample standard deviation s (computed with N − 1 in the denominator, which corrects a small bias). The standard error s/√N estimates how far x̄ typically lies from the true mean. Quadrupling the sample only halves it: the √N is why more runs help, but slowly.

Pseudo-random generators and seeds. Python’s random module uses the Mersenne Twister of Matsumoto and Nishimura (1998), which has a period of 2¹⁹⁹³⁷ − 1. Reusing a seed reproduces the same sequence, which is what makes every table in this series reproducible. Knuth’s warning applies when choosing a generator: “random numbers should not be generated with a method chosen at random.” For experiments, use one well-tested generator and record the seeds.

Skewed distributions, medians and percentiles. Not every distribution is a symmetric hump. Post 19 ran ILS on berlin52 from 100 seeds and recorded the iteration at which each run first found 7,542:

Most runs finish within 100 iterations, but a few take hundreds. That long right tail pulls the mean (62.8) well above the median (41). The median is the middle value when the results are sorted; with an even count, the average of the two middle values, as Python’s statistics.median computes it. The p-th percentile is the value below which p% of the results fall; with 100 runs, the 90th percentile is simply the 90th value in sorted order, 84. For run lengths and running times, report the median and a high percentile, not just the mean. One unlucky seed moved this mean by 7.7 iterations on its own.

Success within a budget. If one try succeeds with probability p, independently of the others, the chance that k tries all fail is (1 − p)ᵏ, so

P(at least one success in k tries) = 1 − (1 − p)ᵏ

In Post 15, 12 of 1,000 random restarts on berlin52 found 7,542, so p ≈ 0.012:

Expecting 1.2 successes in 100 tries still leaves a 30% chance of none.

Comparing two methods. Randomness means two methods’ averages differ even when the methods are equally good. A fair comparison runs both methods on the same instances with the same list of seeds, reports the spread for each, and only calls a difference real when it is clearly larger than that spread. Post 18 did this on a 500-city instance: one rule’s best and worst seeds differed by up to 0.9%, so gaps under 0.6% between rules were called ties. Barr and colleagues (1995) set out guidelines for exactly this kind of reporting: state the experimental conditions in full and make the results reproducible.

flowchart TD
  Q["Question: does method X beat Y?"] --> I["Fix instances and a list of seeds"]
  I --> R["Run X and Y on every instance and seed"]
  R --> S["Summarise each: median, spread, percentiles"]
  S --> C{"Difference larger than the spread?"}
  C -- yes --> W["Report a difference, with the numbers"]
  C -- no --> T["Report a tie, or run more seeds"]

Common misconceptions

“A fixed seed makes the result fair.” A fixed seed makes the result reproducible. One seed is still one sample: Post 19’s seed 31 found 7,542 in 1 iteration and seed 62 in 820, for the same algorithm. Report many seeds, chosen before you look at the results.

“The average run is the typical run.” For skewed results it isn’t. The mean of Post 19’s run lengths, 62.8, is larger than 72 of the 100 runs. The median says what a typical run does; the high percentiles say how bad it can get.

“If p = 1/100, then 100 tries guarantee a success.” They give an expected one success, and a probability of 1 − 0.99¹⁰⁰ ≈ 63% of at least one. A guarantee needs a method that doesn’t depend on luck.

“Method X averaged 0.3% better, so it is better.” Only if 0.3% is large compared with the spread between seeds. In Post 18, differences of that size between acceptance rules were within the seed-to-seed noise.

Try it in Python

The script prices all 2,520 tours of the eight stops to get the exact distribution, draws seeded samples of random tours, checks that an event with probability p happens about p·N times, and tabulates the chance of at least one success in k tries. It uses only the standard library and needs no data files.

import math
import random
import statistics
from itertools import permutations

stops = {"A": (2, 2), "B": (6, 1), "C": (9, 3), "D": (8, 7),
         "E": (5, 9), "F": (1, 8), "G": (4, 5), "H": (10, 9)}
d = {a: {b: round(10 * math.dist(stops[a], stops[b])) for b in stops} for a in stops}
OPTIMAL = 315

def cost(tour):
    return sum(d[tour[i]][tour[(i + 1) % len(tour)]] for i in range(len(tour)))

# 1. The exact distribution: every distinct tour once (A first, mirror images dropped)
everything = [cost("A" + "".join(p)) for p in permutations("BCDEFGH") if p[0] < p[-1]]
print(f"all {len(everything)} tours: mean {statistics.mean(everything)}, "
      f"standard deviation {statistics.pstdev(everything):.1f}, "
      f"min {min(everything)}, max {max(everything)}, "
      f"within 10% of optimal {sum(c <= 1.1 * OPTIMAL for c in everything)}")
for low in range(300, 650, 50):
    print(f"  {low}-{low + 49}: {sum(low <= c < low + 50 for c in everything)} tours")

# 2. Sampling: random tours from a seeded generator, then the sample's mean and spread
def random_tour(rng):
    rest = list("BCDEFGH")
    rng.shuffle(rest)  # every order equally likely
    return "A" + "".join(rest)

for n in (10, 100, 1000):
    rng = random.Random(1)  # same seed, same tours, on any machine
    sample = [cost(random_tour(rng)) for _ in range(n)]
    sd = statistics.stdev(sample)
    print(f"sample of {n:>4}: mean {statistics.mean(sample):.1f}, standard deviation {sd:.1f}, "
          f"standard error {sd / math.sqrt(n):.1f}, optimal tours {sample.count(OPTIMAL)}")

# 3. An event with probability p happens about p * N times in N tries
p = math.exp(-35 / 20)  # Post 18: chance of keeping a tour 35 longer at T = 20
rng = random.Random(1)
kept = sum(rng.random() < p for _ in range(1000))  # rng.random() is uniform on [0, 1)
print(f"p = {p:.3f}: expected {1000 * p:.0f} of 1,000, observed {kept}")

# 4. At least one success in k independent tries, each with probability p
p = 12 / 1000  # Post 15: 12 of 1,000 random restarts found 7,542
for k in (10, 50, 100, 250, 500):
    print(f"k = {k:>3}: expected hits {p * k:.1f}, P(at least one) = {1 - (1 - p) ** k:.3f}")

It prints the tables above. In summary:

Three lessons sit in this output. Sample means approach the true mean, but slowly: the error shrinks with √N. Counts of random events scatter around p·N rather than hitting it. And rare events, like drawing the optimal tour, routinely fail to show up even in large samples.

Exercise

Change the seed in part 2 of the script from 1 to 2, 3, 4 and 5, keeping the sample size at 100. Write down the five sample means. How far do they spread around 488, and how does that compare with the standard error of about 4.5?

Hint: if the standard error describes the typical distance of a sample mean from the true mean, most of your five means should lie within one or two standard errors of 488.

Where to go next

• Post 15: Local optima runs 1,000 random restarts on berlin52 and asks how often they find the optimum.
• Post 18: Acceptance criteria uses u < e^(−Δ/T) to keep worse tours and compares rules across seeds.
• Post 19: Putting ILS together reports a 100-seed run-length distribution with medians and percentiles.
• Post 20: Simulated annealing builds a whole method on random acceptance with a falling temperature.
• Primer E: Python, and reading Rust is the next primer, for the code in Part 5 and later.

References

• Knuth, D. E. (1997). The Art of Computer Programming, Volume 2: Seminumerical Algorithms, 3rd ed. Addison-Wesley. Chapter 3, “Random Numbers”.
• Matsumoto, M., and Nishimura, T. (1998). “Mersenne Twister: A 623-Dimensionally Equidistributed Uniform Pseudo-Random Number Generator.” ACM Transactions on Modeling and Computer Simulation 8(1): 3–30.
• Python Software Foundation. “random — Generate pseudo-random numbers.” Python 3 documentation.
• Wasserman, L. (2004). All of Statistics: A Concise Course in Statistical Inference. Springer.
• Barr, R. S., Golden, B. L., Kelly, J. P., Resende, M. G. C., and Stewart, W. R. (1995). “Designing and Reporting on Computational Experiments with Heuristic Methods.” Journal of Heuristics 1(1): 9–32.
• Birattari, M. (2009). Tuning Metaheuristics: A Machine Learning Perspective. Studies in Computational Intelligence, Springer.

See the whole series · ← All writing