Univariate Analysis

The previous lesson looked at the difference between two region means and wrote its margin of variation next to it: the 4.463 m³ difference appeared in none of 2,000 permutations. The mean itself was never questioned. The 24.674 m³ written for the east region is 557 readings compressed into a single number, and it says nothing about how those 557 numbers are arranged.

This lesson covers how a column is examined on its own, and its unit of measure is this: the shape that the summary hides. Two distributions that give the same mean and the same standard deviation are built, both are reported with the same two numbers, and it is shown that the difference between them never enters the report.

• EA6. Univariate analysis asks a column three questions: where the center is, how scattered the values are, and what the shape of the distribution is. Unless all three are written together, the column has not been described.
• EA7. The five-number summary is the minimum, the first quartile, the median, the third quartile and the maximum. The interquartile range is the distance between Q3 and Q1 and was defined in the Data Collection and Preparation course.
• EA8. Skewness is measured in this course by the range between the mean and the median. If the range is positive, the long tail is on the right; the mean follows the tail, the median does not.
• EA9. Bin counts give the shape of the distribution as numbers: the limits of each bin and the number of values falling into it are printed. The bin width is a decision and it changes the shape.
• EA10. Outliers are not removed in this course. Removal was a cleaning decision of the Data Collection and Preparation course and its unit of measure was the row count; here the same values stay in the set, and only their effect on the summary is counted.

Center, Dispersion and Shape

# distribution.py — MODEL: the fictional measurement network of K03 and K04 is rebuilt
# from the same seed. The dataset is fictional; no real subscriber, meter or region appears.
import math

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


def rng(t: int):
    """Linear congruential generator: seed written in the code, same output every run."""
    x = ((t ^ (t >> 16)) * 2246822507) & M32
    x = ((x ^ (x >> 13)) * 3266489909) & M32
    s = [(x ^ (x >> 16)) & M32]

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


def pick(u: float, w: list[float]) -> int:
    t = 0.0
    for i, x in enumerate(w):
        t += x
        if u < t:
            return i
    return len(w) - 1


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

READINGS = []
for k in _subscribers + _dropped:
    r = rng(SEED + 7000 + k["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(k["base"] * math.exp(
            (r() + r() + r() - 1.5) * 0.62) * (1 - d * 0.05) * 100 + 0.5) / 100
        if k["subscriber_id"] not in _valid:
            r()
            continue
        READINGS.append({"subscriber_id": k["subscriber_id"], "region": k["region"],
                         "period": PERIODS[d], "consumption_m3": t,
                         "household_size": k["household_size"],
                         "satisfaction": k["satisfaction"],
                         "tier": 1 if t < TARIFFS[0] else 2 if t < TARIFFS[1]
                         else 3 if t < TARIFFS[2] else 4})
        r()


def mean(v: list[float]) -> float:
    return sum(v) / len(v)


def std_dev(v: list[float]) -> float:
    """Standard deviation; denominator n-1, because what is at hand is a sample."""
    m = mean(v)
    return math.sqrt(sum((x - m) ** 2 for x in v) / (len(v) - 1))


def median(v: list[float]) -> float:
    s = sorted(v)
    n = len(s)
    return s[n // 2] if n % 2 else (s[n // 2 - 1] + s[n // 2]) / 2


def quartile(v: list[float], p: float) -> float:
    """Linear-interpolation quartile; same as the definition in K03."""
    s = sorted(v)
    i = p * (len(s) - 1)
    a = int(i)
    return s[a] if a + 1 >= len(s) else s[a] + (s[a + 1] - s[a]) * (i - a)


def region_values(name: str) -> list[float]:
    return [o["consumption_m3"] for o in READINGS if o["region"] == name]


east, north = region_values("east"), region_values("north")
print(f"fictional dataset (seed {SEED}): {len(READINGS)} rows, {len(_valid)} subscribers")
print(f"east: {len(east)} rows, mean {mean(east):.3f} m3, "
      f"median {median(east):.2f} m3, standard deviation {std_dev(east):.3f}")
print(f"mean - median = {mean(east) - median(east):.3f} m3, "
      f"{100 * (mean(east) - median(east)) / std_dev(east):.1f} percent of the standard deviation")

fictional dataset (seed 20260218): 3199 rows, 1329 subscribers
east: 557 rows, mean 24.674 m3, median 24.12 m3, standard deviation 9.198
mean - median = 0.554 m3, 6.0 percent of the standard deviation

A column cannot be described with a single number. Central tendency gives two separate numbers: the mean is 24.674 m³ and the median 24.12 m³. Dispersion is the third number: the standard deviation is 9.198, that is, the readings wander around the center by more than nine cubic meters on average. These three numbers do not yet tell the shape, but they give the first sign of it.

The sign is the parting of the mean and the median. The median looks at the ranking and does not care how far away a value is; the mean carries the distance of every value. The 0.554 m³ range between the two says the long tail is on the right: a few large readings have pulled the mean up and left the median in place. The range is 6.0% of the standard deviation and can be reported directly as the measure of skewness.

Same Mean, Different Distribution

# equivalent.py — a second distribution with the same mean and the same standard deviation.
# It is fictional: the sum of three uniform draws is symmetric, then it is fitted to the mean
# and standard deviation of the east. The aim is comparison, not simulation.
def symmetric(n: int, target_mean: float, target_std: float, seed: int) -> list[float]:
    r = rng(seed)
    raw = [r() + r() + r() - 1.5 for _ in range(n)]
    m, s = mean(raw), std_dev(raw)
    v = [round(target_mean + (x - m) * target_std / s, 2) for x in raw]
    shift = target_mean - mean(v)  # rounding shifts the mean a little
    return [round(x + shift, 2) for x in v]


MODEL = symmetric(len(east), mean(east), std_dev(east), SEED + 601)


def five_numbers(v: list[float]) -> list[float]:
    return [min(v), quartile(v, 0.25), median(v), quartile(v, 0.75), max(v)]


HEADERS = ["minimum", "Q1", "median", "Q3", "maximum", "IQR", "mean", "skewness"]
WIDTHS = [17] + [10] * 8
y = lambda h: "".join(str(v).ljust(WIDTHS[i]) if not i else str(v).rjust(WIDTHS[i])
                      for i, v in enumerate(h))
print(y(["distribution"] + HEADERS))
for name, v in [("east (observed)", east), ("model (symmetric)", MODEL)]:
    b = five_numbers(v)
    print(y([name] + [f"{x:.2f}" for x in b] + [f"{b[3] - b[1]:.2f}",
             f"{mean(v):.3f}", f"{mean(v) - median(v):.3f}"]))
print(f"\nboth distributions have a mean of {mean(east):.3f} and a standard deviation of {std_dev(east):.3f}")
print(f"mean - median range: observed {mean(east) - median(east):.3f}, "
      f"model {mean(MODEL) - median(MODEL):.3f}")

distribution        minimum        Q1    median        Q3   maximum       IQR      mean  skewness
east (observed)        0.00     19.23     24.12     29.85     51.90     10.62    24.674     0.554
model (symmetric)      0.66     18.37     24.61     31.43     47.79     13.06    24.674     0.064

both distributions have a mean of 24.674 and a standard deviation of 9.198
mean - median range: observed 0.554, model 0.064

Two rows are reported with the same two numbers. The mean is the same to the third decimal, and so is the standard deviation. The sentence “Mean consumption in the east region is 24.674 m³, standard deviation 9.198” is true for both and does not tell them apart.

The five-number summary does. The minimum of the observation is 0.00, that of the model 0.66. The third quartile of the observation is 29.85, that of the model 31.43. The real difference is in the interquartile range: 10.62 in the observation, 13.06 in the model. Of two distributions that give the same standard deviation, one has a markedly narrower middle half. This means that what enlarges the standard deviation of the observation is not the dispersion in the middle but a few values at the ends.

The last column reduces the measure to a single number. Skewness is 0.554 in the observation and 0.064 in the model — more than eight times. In a symmetric distribution the mean and the median are expected to coincide, and the model does this; the remaining 0.064 is the residue of the draw. In the observation the range is no longer a residue, it is the shape of the distribution itself.

Shape as Numbers

# bins.py — the shape of the distribution is printed as numbers: how many values in each bin.
# No chart is drawn; the bin width is a decision and it changes the shape.
def count_bins(v: list[float], width: int, top: int = 55) -> list[int]:
    k = [0] * (top // width)
    for x in v:
        k[min(int(x // width), len(k) - 1)] += 1
    return k


WIDTH = 5
a5, m5 = count_bins(east, WIDTH), count_bins(MODEL, WIDTH)
print(f"{'bin':>12}{'east':>8}{'model':>8}   east bar")
for i, (x, y2) in enumerate(zip(a5, m5)):
    print(f"[{i * WIDTH:2d},{(i + 1) * WIDTH:2d})".rjust(12) + f"{x:8d}{y2:8d}   "
          + "#" * round(x / 4))
a10 = count_bins(east, 2 * WIDTH)
print(f"\nwidth {WIDTH}: first two bins {a5[0]} and {a5[1]}; "
      f"width {2 * WIDTH}: first bin {a10[0]}")
print(f"readings written as zero {sum(1 for x in east if x == 0)}, "
      f"real readings in the [{WIDTH},{2 * WIDTH}) bin {a5[1]}")

         bin    east   model   east bar
     [ 0, 5)      20       7   #####
     [ 5,10)       1      32
     [10,15)      35      51   #########
     [15,20)     111      83   ############################
     [20,25)     135     112   ##################################
     [25,30)     119      98   ##############################
     [30,35)      71      92   ##################
     [35,40)      35      57   #########
     [40,45)      16      23   ####
     [45,50)      11       2   ###
     [50,55)       3       0   #

width 5: first two bins 20 and 1; width 10: first bin 21
readings written as zero 20, real readings in the [5,10) bin 1

Bin counts show what the five numbers cannot see. There are 20 values in the first bin and 1 value in the second. This is not a tail, it is a break: between the 20 records heaped at zero and the region of real readings sits an almost empty bin. The count says this because the 20 records heaped at zero are rows where zero was written since the meter could not be read; they are not a measured consumption but a record format. The model has no such break: 7 and 32 are consecutive and run smoothly.

The bin width is a decision, and changing this decision changes the shape. When width 10 is chosen, the first bin holds 21 values; 20 and 1 melt into a single number, and the break disappears from the table. The same data, the same procedure, two different shapes. This is why the bin width is written inside the report, not in a footnote.

Outliers Are Not Removed, They Are Measured

# outlier.py — values outside the fences are NOT removed here. In the Data Collection and
# Preparation course removal was a cleaning decision; here the values stay, and the effect is measured.
def fences(v: list[float]) -> tuple[float, float]:
    q1, q3 = quartile(v, 0.25), quartile(v, 0.75)
    iqr = q3 - q1
    return q1 - 1.5 * iqr, q3 + 1.5 * iqr


WIDTHS = [9, 6, 12, 12, 10, 10, 10, 11, 8]
row = lambda h: "".join(str(v).ljust(WIDTHS[i]) if not i else str(v).rjust(WIDTHS[i])
                        for i, v in enumerate(h))
print(row(["dataset", "n", "lower fence", "upper fence", "below", "above", "mean",
           "kept mean", "diff"]))
kept = {}
for name, v in [("east", east), ("north", north), ("model", MODEL)]:
    lower, upper = fences(v)
    kept[name] = [x for x in v if lower <= x <= upper]
    print(row([name, len(v), f"{lower:.2f}", f"{upper:.2f}",
               sum(1 for x in v if x < lower), sum(1 for x in v if x > upper),
               f"{mean(v):.3f}", f"{mean(kept[name]):.3f}",
               f"{mean(kept[name]) - mean(v):+.3f}"]))
print(f"\neast / north difference: with all rows {mean(east) - mean(north):.3f} m3, "
      f"if the outside-fence values were removed {mean(kept['east']) - mean(kept['north']):.3f} m3")
print(f"east standard deviation {std_dev(east):.3f} -> {std_dev(kept['east']):.3f}, "
      f"median {median(east):.2f} -> {median(kept['east']):.2f}")

dataset       n lower fence upper fence     below     above      mean  kept mean    diff
east        557        3.30       45.78        20        10    24.674     25.154  +0.480
north       921        2.41       37.96        32        18    20.210     20.532  +0.321
model       557       -1.22       51.02         0         0    24.674     24.674  +0.000

east / north difference: with all rows 4.463 m3, if the outside-fence values were removed 4.622 m3
east standard deviation 9.198 -> 7.416, median 24.12 -> 24.22

The values outside the fences are not deleted here; the table only counts what would happen if they were. In the east there are 30 rows outside the fences, and 20 of them are outside the lower fence. Removal does not lower the mean, it raises it from 24.674 to 25.154, because most of the removed rows are readings written as zero. This is where counting an outlier as a “large value” misleads.

The same decision moves the mean in the north by 0.321 m³, and since the two movements do not go the same way, the regional difference changes too: 4.622 m³ instead of 4.463 m³. The difference measured in the previous lesson has shifted by 0.159 m³ with no new data, through a single decision about the fence. This is a variation separate from the margin that comes from the sample, and added to it.

The model row makes clear what the measure is. In a symmetric distribution there is no value outside the fences, and removal leaves the mean exactly where it is. The number of values outside the fences is a consequence not of the data being “faulty” but of the shape of the distribution; the number 30 measures the tails of the east, not its quality.

The last row shows why the median is called resistant. When thirty rows are taken out, the standard deviation falls from 9.198 to 7.416, that is, by 19.4%; the median moves from 24.12 to 24.22, that is, by 0.10. The measure of dispersion depends on the ends, the measure of ranking does not.

Summary

• A column is described with three measures: for the east, a mean of 24.674 m³, a median of 24.12 m³ and a standard deviation of 9.198; the 0.554 m³ range between them is the first sign of shape.
• A symmetric model distribution that gives the same mean and the same standard deviation diverges in the five-number summary: an interquartile range of 10.62 against 13.06, skewness of 0.554 against 0.064.
• Bin counts give what the five numbers cannot see: the observation has 20 values in the first bin and 1 in the second, and this break comes from the readings written as zero.
• The bin width is a decision: when width 10 is chosen, the first bin becomes 21 and the break is not visible in the table.
• The 30 rows outside the fences are not removed; removal would have moved the east mean to 25.154 and the difference between east and north from 4.463 to 4.622 m³, and lowered the standard deviation by 19.4%.

Next Step

This lesson looked at a single column and wrote its shape as numbers. Most of the questions in the dataset cannot be answered with a single column, though: do consumption and household size rise together, is there a link between tariff tier and region, how does the satisfaction level relate to consumption. Examining two columns together brings a new summary — the correlation coefficient — and that coefficient, just like the mean, compresses everything into a single number. The next lesson covers bivariate analysis: two different relationships that give the same coefficient are built, a cross tabulation shows what the coefficient cannot see, and it is written why rank correlation is the right measure on an ordinal scale.