Distribution Comparison

The previous lesson chose a single number and wrote its margin. Most of the questions in a report ask not for a single number but for a comparison: is the consumption distribution of two regions different from each other. By habit this comparison is made over five numbers — the minimum, the first quartile, the median, the third quartile, the maximum — and the five numbers are put side by side as a box summary.

In this lesson the box summary and the density display are not drawn; both are built as data structures and printed as numbers. Drawing, axis and design decisions are the job of M26/K06. The question here comes earlier: do five numbers determine a distribution. The measure is how many peaks two datasets whose five numbers are identical to each other give in their bin counts.

• AS6. The box summary is five numbers and two fences. The interquartile range is the difference between the third quartile and the first quartile; the fences stand outside by 1.5 times this range. In M26/K03 these five numbers were an elimination decision; here there is no elimination, the numbers are used to describe the shape of the distribution.
• AS7. Five numbers do not determine a distribution. There are several datasets that give the same five numbers, and the difference among them is not reflected in the box summary at all.
• AS8. A bin count is the numeric counterpart of density. The value axis is divided into equal-width bins and the observations falling into each bin are counted; the number of peaks is read from these counts.
• AS9. The notch is the margin of the median and is produced with bootstrap: the median is recomputed in 1,000 draws, and the middle 95% of the resulting values is the notch.
• AS10. The width of the notch depends on the number of observations and on the density around the median. The same median gives a much wider notch in a dataset where the neighborhood of the median is sparse.
• AS11. Whether two notches intersect is a display cue, not a test; what the intersection of intervals means is the subject of the next lesson.

The Box Summary Is a Data Structure

# dataset.py — MODEL: the measurement network of K03-K04 is regenerated with the same seed.
# A long-format reading table of 3,199 rows, 1,329 subscribers, three periods, five regions.
import math

SEED, RAW, M32 = 20260218, 1400, 0xFFFFFFFF
PERIODS = ["2026-01", "2026-02", "2026-03"]
REGION = [("north", 0.28, 21), ("south", 0.22, 17), ("east", 0.18, 26),
          ("west", 0.14, 14), ("central", 0.18, 23)]


def rng(t):
    x = ((t ^ (t >> 16)) * 2246822507) & M32
    x = ((x ^ (x >> 13)) * 3266489909) & M32
    s = [(x ^ (x >> 16)) & M32]

    def next_value():
        s[0] = (s[0] * 1664525 + 1013904223) & M32
        return s[0] / 4294967296
    return next_value


def pick(u, w):
    t = 0.0
    for i, x in enumerate(w):
        t += x
        if u < t:
            return i
    return len(w) - 1


_kept, _dropped = [], []
for i in range(RAW):
    r = rng(SEED + i)
    b = REGION[pick(r(), [x[1] for x in REGION])]
    pick(r(), [0.06, 0.24, 0.30, 0.24, 0.11, 0.05])
    pick(r(), [0.08, 0.14, 0.27, 0.34, 0.17])
    (_dropped if r() < 0.046 else _kept).append(
        {"subscriber_id": 10001 + i, "region": b[0], "base": b[2]})
VALID = {a["subscriber_id"] for a in _kept}

READINGS = []
for a in _kept + _dropped:
    r = rng(SEED + 7000 + a["subscriber_id"])
    for d in range(pick(r(), [0.05, 0.12, 0.21, 0.62])):
        t = 0.0 if r() < 0.038 else math.floor(
            a["base"] * math.exp((r() + r() + r() - 1.5) * 0.62)
            * (1 - d * 0.05) * 100 + 0.5) / 100
        if a["subscriber_id"] not in VALID:
            r()
            continue
        READINGS.append({"subscriber_id": a["subscriber_id"], "region": a["region"],
                         "period": PERIODS[d], "consumption_m3": t})
        r()


def percentile2(v, p):
    s = sorted(v)
    i = p * (len(s) - 1)
    a = int(i)
    b = min(a + 1, len(s) - 1)
    return s[a] + (s[b] - s[a]) * (i - a)


def five_numbers(v):
    return [min(v), percentile2(v, 0.25), percentile2(v, 0.50),
            percentile2(v, 0.75), max(v)]


def box(v):
    b = five_numbers(v)
    iqr = b[3] - b[1]
    return b, iqr, b[1] - 1.5 * iqr, b[3] + 1.5 * iqr


north = sorted(o["consumption_m3"] for o in READINGS if o["region"] == "north")
b, iqr, fence_low, fence_high = box(north)
print(f"north: {len(north)} readings")
print("five-number summary: " + " | ".join(f"{x:.2f}" for x in b))
print(f"interquartile range {iqr:.2f} m3, fences {fence_low:.3f} and {fence_high:.3f} m3")
print(f"{sum(1 for x in north if x < fence_low or x > fence_high)} readings outside the fences")

north: 921 readings
five-number summary: 0.00 | 15.74 | 19.66 | 24.63 | 45.47
interquartile range 8.89 m3, fences 2.405 and 37.965 m3
50 readings outside the fences

The box summary here is not a drawing but a record of seven numbers: five numbers, the interquartile range and two fences. The 921 readings of the north region came down to these seven. The median is 19.66 m³, the interquartile range is 8.89 m³, the fences are between 2.405 and 37.965 m³, and 50 readings stay outside these limits.

That the smallest value of the dataset is zero suggests that most of the readings below the lower fence are zero readings. M26/K03 had dealt with these readings as a cleaning decision; here none of them is eliminated, only how many of them stay outside the fences is counted. The fence is not an elimination threshold but a measure of shape.

Seven numbers take little space in a report and speed up the comparison. Their price is measured in the next section: these seven numbers say much less than one would think about how the 921 readings are arranged.

The Same Five Numbers, Another Dataset

# twin.py — SETUP: the sorted values of north are pushed, inside each quarter slice, toward the end
# of the slice. Since the 0th, 230th, 460th, 690th and 920th values in the sorted array stay in
# place, the five-number summary is preserved exactly, the inner distribution changes.
def push(part, end, p=4):
    a, z = part[0], part[-1]
    u = [(x - a) / (z - a) for x in part]
    return [z - (z - a) * (1 - t) ** p if end == "upper" else a + (z - a) * t ** p
            for t in u]


C = [0, 230, 460, 690, len(north) - 1]
twin = (push(north[C[0]:C[1] + 1], "upper")[:-1]
        + push(north[C[1]:C[2] + 1], "lower")[:-1]
        + push(north[C[2]:C[3] + 1], "upper")[:-1]
        + push(north[C[3]:], "lower"))
b2, iqr2, low2, high2 = box(twin)
print(f"twin dataset: {len(twin)} values (fictional, derived from the sorted values of north)")
print("five-number summary: " + " | ".join(f"{x:.2f}" for x in b2))
print(f"interquartile range {iqr2:.2f} m3, fences {low2:.3f} and {high2:.3f} m3")
print(f"{sum(1 for x in twin if x < low2 or x > high2)} values outside the fences")
print(f"are the two box summaries the same: {[round(x, 2) for x in b] == [round(x, 2) for x in b2]}")

twin dataset: 921 values (fictional, derived from the sorted values of north)
five-number summary: 0.00 | 15.74 | 19.66 | 24.63 | 45.47
interquartile range 8.89 m3, fences 2.405 and 37.965 m3
35 values outside the fences
are the two box summaries the same: True

The twin dataset is fictional and how it was built is in the code: the sorted array is divided into four quarter slices, and the values inside each slice are pushed toward one end of the slice. The push does not break the ordering and the values that coincide with the slice ends stay in place, so the five numbers are preserved exactly. Because the push direction of the lower and the upper half was chosen opposite, the values pile up on the two sides of the median.

The two box summaries are the same: the same median, the same quartiles, the same interquartile range, the same fences. A report that puts the box summaries side by side cannot report any difference between the two datasets. The only thing it can report is in the last line, and that is not among the five numbers either: the number of observations outside the fences has dropped from 50 to 35, that is, the same fences leave a different number of values outside in the two datasets.

This is not a defect of the box summary but its scope. The five numbers report five points of the ordering; they say nothing about how the values are distributed between those points, and they were not designed to. The information the box summary carries is independent of the position in the ordering of all but four of the 921 readings: every arrangement inside a slice gives the same seven numbers.

It must be written once more that the twin dataset is a fictional setup with no counterpart in the real network. Its purpose is not to imitate a region but to show the size of the difference the box summary cannot distinguish. The same effect appears in real datasets in a more moderate form, and it usually arises from two different subscriber groups being collected into a single row.

Bin Counts

# bins.py — the density display as numbers: the value axis is divided into 4 m3 bins,
# the observations in each bin are counted. The text bar is one mark per ten observations.
WIDTH, N = 4.0, 12


def counts(v):
    c = [0] * N
    for x in v:
        c[min(int(x / WIDTH), N - 1)] += 1
    return c


def peaks(c):
    return [i for i in range(1, N - 1)
            if c[i] > c[i - 1] and c[i] > c[i + 1] and c[i] >= 0.3 * max(c)]


cn, ct = counts(north), counts(twin)
for i in range(N):
    print(f"{i * WIDTH:5.0f}-{(i + 1) * WIDTH:<3.0f} {cn[i]:4d} {'#' * round(cn[i] / 10):<34}"
          f" {ct[i]:4d} {'#' * round(ct[i] / 10)}")
print(f"peak count: north {len(peaks(cn))}, twin {len(peaks(ct))}")
print(f"the peaks of the twin are in the {[int(i * WIDTH) for i in peaks(ct)]} m3 bins, "
      f"while the median {b2[2]:.2f} m3 is in the trough")

    0-4     32 ###                                  32 ###
    4-8      0                                       0
    8-12    41 ####                                  0
   12-16   175 ##################                  318 ################################
   16-20   234 #######################             114 ###########
   20-24   182 ##################                  110 ###########
   24-28   130 #############                       329 #################################
   28-32    66 #######                              11 #
   32-36    37 ####                                  4
   36-40    16 ##                                    0
   40-44     5                                       2
   44-48     3                                       1
peak count: north 1, twin 2
the peaks of the twin are in the [12, 24] m3 bins, while the median 19.66 m3 is in the trough

The two datasets that the box summary could not separate were separated with twelve pairs of numbers. In north the counts rise to a single peak in the 16–20 m³ bin and fall off to both sides. In the twin there are two peaks: 318 values in the 12–16 m³ bin, 329 values in the 24–28 m³ bin, and between them stands a trough dropping to 114 and 110.

The place of the trough matters: the median is 19.66 m³ and falls exactly into this trough. In the twin dataset the values closest to the median are the sparsest ones. The sentence “the typical subscriber consumes 19.66 m³” is defensible for north, but for the twin dataset it is true and still misleading — there is almost nobody around that value.

The bin width is a decision, and this decision itself moves the number: if 10 m³ had been chosen instead of 4 m³, the two peaks would fall into a single bin and the difference would again be invisible. The visual effect and the choice of bin width are the subject of M26/K06; the only thing recorded here is that the chosen width is written under the output.

Notch Interval

# notch.py — the margin of the median: drawn 1,000 times with replacement, the median is
# recomputed at every draw, and the middle 95% of the resulting values is the notch interval.
DRAWS = 1000


def percentile(s, p):
    i = p * (len(s) - 1)
    a = int(i)
    z = min(a + 1, len(s) - 1)
    return s[a] + (s[z] - s[a]) * (i - a)


def notch(v, seed, draws=DRAWS):
    rand = rng(seed)
    n = len(v)
    d = sorted(percentile2([v[int(rand() * n)] for _ in range(n)], 0.50)
               for _ in range(draws))
    return percentile(d, 0.025), percentile(d, 0.975)


central = [o["consumption_m3"] for o in READINGS if o["region"] == "central"]
for name, v in (("north", north), ("twin", twin), ("central", central)):
    a, z = notch(v, SEED + 211)
    print(f"{name:7s} n={len(v):4d} median {percentile2(v, 0.50):6.3f} "
          f"notch {a:.3f} - {z:.3f} width {z - a:.3f}")
na, nz = notch(north, SEED + 211)
ca, cz = notch(central, SEED + 211)
print(f"do the notches of north and central intersect: {na < cz and ca < nz}")

north   n= 921 median 19.660 notch 19.230 - 20.150 width 0.920
twin    n= 921 median 19.660 notch 18.203 - 21.349 width 3.146
central n= 600 median 21.535 notch 20.795 - 22.330 width 1.535
do the notches of north and central intersect: False

The notch is the margin of variation of the median and is produced with the very same bootstrap as in the previous lesson. The medians of north and the twin dataset are the same number and so are their observation counts, yet their notches are 0.920 and 3.146 m³ wide. The three and a half times in between come not from the number of observations but from the density around the median: a median standing in a trough goes much farther from draw to draw.

This is the numeric expression of what the box summary lacks. While two boxes that put the five numbers side by side cannot be told apart, the margin of the same median differs by more than three times. A report that writes the notch next to the median reports, even without drawing anything, that the two distributions are not the same. The notch here is the sixth and seventh number added on top of the five, and its cost is a single run of 1,000 draws.

The last line is a real comparison: the notches of the north and central regions do not intersect. Whether this means “the medians of the two regions are different” is not decided here. The interpretation built on intersection and non-intersection is the first question of the next lesson.

Summary

• The box summary is a record, not a drawing: five numbers, the interquartile range and two fences; the 921 readings of north come down to the values 0.00 | 15.74 | 19.66 | 24.63 | 45.47 and an interquartile range of 8.89 m³.
• A second dataset that gives the same five numbers can be built; the two box summaries are exactly the same and the difference shows only in the number of observations outside the fences, 50 against 35.
• Bin counts separate the two datasets: in north a single peak at 16–20 m³, in the twin two peaks in the 12–16 and 24–28 m³ bins, and the median is inside the trough with counts of 114 and 110.
• The notch of the median is produced with bootstrap and comes out 0.920 m³ and 3.146 m³ wide for the same median; the width depends on the number of observations and on the density around the median.
• The notches of north and central do not intersect; this is a display cue and by itself not a test result.

Next Step

In this lesson intervals were produced twice and widths were written twice, but what an interval says was never defined. The confidence interval was used as a tool in M15/K02; the next lesson takes the interval itself as its subject. How the width depends on the number of observations is measured — what happens to the width when the number of observations is quadrupled is counted — and what the interval does not say is written out plainly. In the end the intersection of two intervals is dealt with: intersection does not mean no difference, and the proof of that is computing the difference’s own interval.