Aggregation and the axis Argument

Context

An aggregation collapses many values into one: a total, a mean, an extreme, a spread, a middle value. On a table it can collapse everything at once or only one direction, and the axis argument names the dimension that disappears, so the same call returns one value per feature or one per listing.

Summaries come first with any dataset: each feature's mean and spread to see its scale, the median price, the cutoff for the cheapest quarter, the position of the most expensive listing. They come last too, because a model's error is a mean over its residuals, and standardizing features, the starting point of Module 3's scaling lesson, is a per-feature mean and standard deviation applied back to the table.

The wrong axis returns a plausible array and no error, of the wrong length, or on a square table of the right one. A single missing value turns NumPy's mean into nan, and NumPy's standard deviation divides by the number of values where the table library of Section 2.4 divides by one less, so one column can report two spreads. The previous lesson stretched small arrays to meet a table, and this one collapses a table into small arrays.

Goal

• Summarize an array with sum, mean, min, max, var, std, np.median and np.percentile, and find an extreme's position with argmax and argmin.
• Collapse a table per feature with axis=0 or per listing with axis=1, and predict each result's shape before running it.
• Subtract a per-row summary from its table with keepdims=True, and standardize every feature of a table in one line.
• Count missing values with np.isnan and skip them with np.nanmean, np.nansum and np.nanmedian.
• Choose between ddof=0 and ddof=1, and state what each one divides by.

Prerequisites

Nothing here is required, and the lesson works without it. Both are here if you want to see an axis collapse a small table before you start.

• YouTube - "Complete Python NumPy Tutorial (Creating Arrays, Indexing, Math, Statistics, Reshaping)" by Keith Galli, the "Statistics" chapter at 42:19. np.min on a small 2 by 3 array, then np.max and np.sum with an axis.
• Official docs - NumPy, "the absolute basics for beginners", "More useful array operations". Aggregating a whole array, then each column with axis=0.

Table of Contents

• Whole-Array Aggregation
• The axis Argument
• keepdims and Standardizing Columns
• Missing Values and ddof

Theory

Whole-Array Aggregation

What. An aggregation reduces an array to one number. sum, mean, min and max are methods of every array, var is the mean of the squared distances from the mean, and std is its square root, back in the data's own units. np.median is the middle of the sorted values, np.percentile and np.quantile read any cut point, and argmin and argmax return a position, not a value. Why. A column of prices says little at a glance, and its total, center and spread say a lot. A position matters as much as a value, because argmax tells you which listing holds the highest price, so the rest of that listing can be read off the table. When. Whenever one number must describe a whole array, a total, a typical value, a spread or a cutoff, and argmax or argmin whenever the question is which row holds the extreme.

import numpy as np

# Sale prices of six listings, in thousands.
price = np.array([305, 248, 410, 200, 265, 342])

print(price.sum(), price.mean())       # 1770 295.0
print(price.min(), price.max())        # 200 410
print(price.var())                     # 4611.333333333333
print(price.std())                     # 67.90679887414318
print(np.median(price))                # 285.0
print(np.percentile(price, [25, 75]))  # [252.25 332.75]
print(np.quantile(price, 0.25))        # 252.25
print(price.argmax(), price.argmin())  # 2 3
print(price[price.argmax()])           # 410

var is in thousands squared, which is why it reads as 4611, and std takes the square root to return to thousands, so prices here spread about 68 thousand around their mean.

np.median sorts the six prices and averages the middle two, 265 and 305. np.percentile(price, 25) reads the point 25 percent of the way along the sorted positions, which lands a quarter of the way from 248 to 265, at 252.25, and the 50th percentile is the median. np.quantile takes the same cut as a fraction, 0.25 where np.percentile takes 25.

argmax returned 2, the position of the most expensive listing, and indexing with it returns the price itself.

DON'T:

• DON'T read argmax as the largest value. price.argmax() is the position 2, and on a 2-D table without an axis it counts positions through the flattened table, so the (4, 3) table of the next subtitle gives 8.
• DON'T call price.median(). Arrays have no such method, and the line raises AttributeError: 'numpy.ndarray' object has no attribute 'median'.
• DON'T hand np.percentile a fraction. np.percentile(price, 0.25) is the 0.25th percentile, 200.6, nearly the minimum, where np.quantile(price, 0.25) is 252.25.

The axis Argument

What. Every aggregation takes an axis argument that names the dimension to collapse. On an (n, k) table, axis=0 collapses the rows and returns one value per column, per feature, and axis=1 collapses the columns and returns one value per row, per listing. The axis you name is the one that disappears from the shape. Why. A table has two summaries, and both are routine: each source's average estimate down its column, and each listing's consensus across its row. Without axis, an aggregation collapses everything into one number that mixes listings and sources. When. Whenever a table's summary should be per feature, axis=0, as for the means and spreads that scaling needs, or per listing, axis=1, as for a listing's consensus or total.

The table holds four listings and three price estimates for each, the estimates Lesson 2.2.4 corrected and compared:

import numpy as np

# Four listings, three estimates each, in thousands:
# agent's appraisal, tax assessment, online estimate.
est = np.array([[310, 280, 325],
                [245, 231, 250],
                [420, 365, 442],
                [199, 186, 209]])

print(est)
print(est.shape, est.mean())
print(est.mean(axis=0), est.mean(axis=0).shape)
print(est.mean(axis=1), est.mean(axis=1).shape)

[[310 280 325]
 [245 231 250]
 [420 365 442]
 [199 186 209]]
(4, 3) 288.5
[293.5 265.5 306.5] (3,)
[305. 242. 409. 198.] (4,)

axis=0 names the first dimension, the four listings, and removes it, leaving one average per source: the online estimate averages 306.5 and the tax assessment 265.5. axis=1 removes the three sources, leaving one consensus per listing. Every aggregation from the first subtitle takes axis the same way: est.max(axis=1) is each listing's highest estimate, [325 250 442 209], and est.argmax(axis=1) is the source that gave it, [2 2 2 2], the online estimate every time.

A wrong axis raises nothing. Asking for each source's average with axis=1 returns four consensus values where three averages were wanted, plausible numbers of the wrong length, and on a square table even the length is right. Add these lines to the same script:

est3 = est[:3]            # three listings, square
print(est3.mean(axis=0))  # [325. 292. 339.]
print(est3.mean(axis=1))  # [305. 242. 409.]

Both results are (3,) and both look like prices. Only the axis says that the first is one average per source and the second one consensus per listing.

DON'T:

• DON'T read axis=0 as "along the columns, so per row". The axis you name disappears, so axis=0 removes the rows and leaves one value per column.
• DON'T take a result's length as proof of the right axis. On a square table both axes return the same length, so decide which dimension the values belong to before writing the call.

keepdims and Standardizing Columns

What. keepdims=True keeps the collapsed axis as a dimension of size 1, so est.mean(axis=1, keepdims=True) is (4, 1) rather than (4,). Standardizing a table, (X - X.mean(axis=0)) / X.std(axis=0), gives every feature a mean of 0 and a standard deviation of 1, except a constant one, whose standard deviation of 0 turns every value into nan with a warning. Why. A summary is often subtracted back from the table it came from, and Lesson 2.2.4's rule decides whether that works: a per-column result, (k,), pairs with (n, k) as it is, and a per-row result needs a column shape, which keepdims keeps for you. Standardizing uses the data's own center and spread to put features with different units on one scale. When. keepdims whenever a per-row summary goes back into an operation with its table, and standardization before anything that compares features across columns, Lesson 2.3.3's distances and Module 3's distance-based models among them.

import numpy as np

est = np.array([[310, 280, 325],
                [245, 231, 250],
                [420, 365, 442],
                [199, 186, 209]])

consensus = est.mean(axis=1, keepdims=True)
print(consensus.shape)
print(est - consensus)

(4, 1)
[[  5. -25.  20.]
 [  3. -11.   8.]
 [ 11. -44.  33.]
 [  1. -12.  11.]]

Each row now says how far each source sits from that listing's consensus, and the tax assessment is below it on all four. Without keepdims, est.mean(axis=1) is (4,), and est - est.mean(axis=1) raises Lesson 2.2.4's ValueError with the shapes (4,3) (4,), or, on a square table, subtracts per column without a word. est.mean(axis=1)[:, np.newaxis] is the Lesson 2.2.4 route to the same (4, 1) column.

Standardizing works down the columns instead, so its (3,) results broadcast as they are:

import numpy as np

# Four listings: area_sqm, rooms, age.
X = np.array([[ 85, 3, 12],
              [ 62, 2, 30],
              [120, 4,  5],
              [ 48, 1, 41]])

mu = X.mean(axis=0)    # one mean per feature, (3,)
sigma = X.std(axis=0)  # one spread per feature, (3,)
Z = (X - mu) / sigma
print(mu)
print(sigma)
print(np.round(Z, 3))
print(Z.mean(axis=0))
print(np.round(Z.mean(axis=0), 3), Z.std(axis=0))

[78.75  2.5  22.  ]
[27.23394022  1.11803399 14.26534262]
[[ 0.229  0.447 -0.701]
 [-0.615 -0.447  0.561]
 [ 1.515  1.342 -1.192]
 [-1.129 -1.342  1.332]]
[ 0.00000000e+00  0.00000000e+00 -5.55111512e-17]
[ 0.  0. -0.] [1. 1. 1.]

A standardized value counts standard deviations from the feature's mean: the third listing's area sits 1.515 above the average area and its age 1.192 below the average age, numbers that compare across features in a way 120 square meters and 5 years never could. The age column's mean came out as -5.55111512e-17 rather than 0, ordinary float noise from a sum that does not cancel exactly, so the check is made on rounded values. Rounding leaves -0., a zero that keeps the sign of the tiny negative value it came from, and -0.0 == 0 is True.

DON'T:

• DON'T subtract a per-row summary without keepdims=True or a new axis. (n, k) against (n,) returns an (n, n) table when k is 1, and otherwise raises when n differs from k and works per column, wrongly, when they are equal.
• DON'T read -5.55e-17 or -0. as a failed standardization. Both are zero up to float noise.
• DON'T standardize with axis=1. That takes each listing's mean across area, rooms and age, three different units.

Missing Values and ddof

What. np.nanmean, np.nansum and np.nanmedian skip a nan, the value Lesson 2.2.3 showed turning a sum or a mean into nan, and np.isnan(x).sum() counts it. std and var divide the summed squared distances by n by default, ddof=0, and ddof=1 divides by n minus 1, ddof being short for delta degrees of freedom, the amount taken off n. Why. Real data arrives with gaps, the housing set of Section 2.4 among it, so a summary either fails loudly with nan or skips the gap, and you need to know which one happened. The divisor matters because the same column gives two different spreads under the two conventions. When. The nan functions when a statistic over the known values is what you want, after counting the missing ones. ddof=1 when the data is a sample and the spread of the larger set it came from is what you estimate, and whenever you compare with a tool that uses it.

import numpy as np

# One listing's price is missing.
price = np.array([305, 248, np.nan, 200, 265, 342])

print(np.isnan(price).sum())                # 1
print(np.nanmean(price), np.nansum(price))  # 272.0 1360.0
print(np.nanmedian(price))                  # 265.0

np.isnan(price) is a boolean mask, Lesson 2.2.2's, and sum counts its True values. The nan functions return the statistics of the five known prices and say nothing about a sixth, which is why the count comes first. Neither changes the array: the nan is still in price, and deciding what to do with a gap is Lesson 2.5.1's subject.

import numpy as np

price = np.array([305, 248, 410, 200])
sq = (price - price.mean()) ** 2        # squared distances

print(price.var(), sq.sum() / 4)        # 6121.6875 6121.6875
print(price.var(ddof=1), sq.sum() / 3)  # 8162.25 8162.25
print(price.std(), price.std(ddof=1))   # 78.24121356420797 90.34517142603693

The default divides the summed squared distances by the 4 prices and ddof=1 by 3, and on four values the two spreads differ by more than 12 thousand.

Dividing by n describes these prices themselves, and dividing by n minus 1 is the usual estimate of the spread in the larger set of listings they were drawn from. The table library installed in Section 2.4 defaults to ddof=1, and Lesson 2.7.1 sets the two side by side.

DON'T:

• DON'T trust a mean before counting the nan values. One nan turns the plain mean into nan, and the nan functions return a number that hides how many were missing.
• DON'T compare standard deviations from two tools without checking ddof. The same four prices give 78.24 with ddof=0 and 90.35 with ddof=1.
• DON'T treat np.nanmean as cleaning. It skips the gap in one calculation and leaves the nan in the array.

Lab Exercises

Lab Exercise 1: Summarize the Prices, Then Each Source and Each Listing

Seven sale prices, in thousands, are [305, 248, 410, 200, 265, 342, 289]. Print their total, mean, median, 25th and 75th percentiles and standard deviation, then the position and the value of the cheapest listing.

Then take five listings' estimates from an agent, the tax office and an online service, in thousands:

import numpy as np

est5 = np.array([[280, 255, 300],
                 [350, 310, 372],
                 [199, 180, 215],
                 [460, 400, 480],
                 [315, 290, 330]])

Before running anything, write in a comment the shape you expect for each source's mean estimate and for each listing's mean estimate. Print both, then use argmax on the right one to print the full row of the listing with the highest mean estimate.

Lab Exercise 2: Standardize, Center, Then Count the Gaps

Standardize these five listings' features, area_sqm, rooms and age, and print the result rounded to 3 places together with its rounded column means and standard deviations:

import numpy as np

X5 = np.array([[110, 4, 18],
               [ 58, 2, 35],
               [205, 6,  3],
               [ 76, 3, 52],
               [ 40, 1, 90]])

Then center each listing's estimates in est5 on that listing's own mean, keeping the collapsed dimension, and print the rounded row means of the result. Last, for the price column [276, 340, np.nan, 450, 318, np.nan, 289], count the missing values, print the plain and the nan-skipping mean, and print the standard deviation of the known prices with ddof=0 and with ddof=1.

import numpy as np

price7 = np.array([305, 248, 410, 200, 265, 342, 289])
print(price7.sum(), price7.mean(), np.median(price7))
print(np.percentile(price7, [25, 75]), price7.std())
print(price7.argmin(), price7[price7.argmin()])

est5 = np.array([[280, 255, 300],
                 [350, 310, 372],
                 [199, 180, 215],
                 [460, 400, 480],
                 [315, 290, 330]])
per_source = est5.mean(axis=0)   # expect (3,)
per_listing = est5.mean(axis=1)  # expect (5,)
print(per_source)
print(per_listing)
top = per_listing.argmax()
print(top, est5[top])

2059 294.14285714285717 289.0
[256.5 323.5] 62.904527594098774
3 200
[320.8 287.  339.4]
[278.33333333 344.         198.         446.66666667 311.66666667]
3 [460 400 480]

import numpy as np

X5 = np.array([[110, 4, 18],
               [ 58, 2, 35],
               [205, 6,  3],
               [ 76, 3, 52],
               [ 40, 1, 90]])
Z5 = (X5 - X5.mean(axis=0)) / X5.std(axis=0)
print(np.round(Z5, 3))
print(np.round(Z5.mean(axis=0), 3), np.round(Z5.std(axis=0), 3))

est5 = np.array([[280, 255, 300],
                 [350, 310, 372],
                 [199, 180, 215],
                 [460, 400, 480],
                 [315, 290, 330]])
centered = est5 - est5.mean(axis=1, keepdims=True)
print(np.round(centered.mean(axis=1), 3))

[[ 0.209  0.465 -0.718]
 [-0.682 -0.697 -0.153]
 [ 1.836  1.627 -1.217]
 [-0.373 -0.116  0.412]
 [-0.99  -1.279  1.676]]
[ 0. -0.  0.] [1. 1. 1.]
[ 0.  0.  0. -0. -0.]

import numpy as np

pm = np.array([276, 340, np.nan, 450, 318, np.nan, 289])
print(np.isnan(pm).sum(), pm.mean(), np.nanmean(pm))
known = pm[~np.isnan(pm)]  # Lesson 2.2.2's mask, negated
print(known.std(), known.std(ddof=1))

2 nan 334.6
61.86307460836392 69.16502006072145

Theory Assessment

Q1: Value or Position, and Which Scale

For price = np.array([305, 248, 410, 200, 265, 342]), say what price.argmax() returns and how you get the most expensive price from it, why price.median() fails and which call works, and why np.percentile(price, 0.25) and np.quantile(price, 0.25) return different numbers.

Hint: One of the calls returns a position. For the last pair, ask what scale each function reads its second argument on.

Q2: The Axis That Disappears

A table X holds 5 listings by 3 features. Give the shapes of X.sum(), X.sum(axis=0), X.sum(axis=1) and X.mean(axis=1, keepdims=True), say which one gives each feature's total, and explain why a wrong axis raises no error and when even the result's shape cannot catch it.

Hint: Cross out the dimension each call names. keepdims leaves a 1 where the crossed-out size was.

Q3: One Gap, Two Divisors

For [305, 248, np.nan, 200, 265, 342], say what mean, np.nanmean and np.isnan(...).sum() return, and why the count belongs before the nan-skipping mean. Then for [305, 248, 410, 200], say which of std() and std(ddof=1) is larger, what each divides by, and why another tool can print the other number for the same prices.

Hint: nan spreads through arithmetic. For the divisor, write the variance as a sum of squared distances divided by something.

Q4: A Standardization That Ran Without Complaint

A teammate standardizes three listings' features, area_sqm, rooms and age, with these lines, and nothing raises:

import numpy as np

X3 = np.array([[85, 3, 12], [62, 2, 30], [120, 4, 5]])
Z = (X3 - X3.mean(axis=1)) / X3.std(axis=1)

Say what the two statistics measure, which dimension of X3 the broadcast then pairs them with, what printing the rounded Z.mean(axis=0) would show, and which line they meant to write.

Hint: Two mistakes stack here. First ask what axis=1 collapses on this table, then use Lesson 2.2.4's rule to ask which dimension a (3,) result lines up with.

Assignment

Assignment: A Weighted Consensus That Survives a Missing Estimate

Three sources value five listings, and the online service failed on the second listing, so that estimate is missing:

import numpy as np

est5 = np.array([[280, 255, 300],
                 [350, 310, np.nan],
                 [199, 180, 215],
                 [460, 400, 480],
                 [315, 290, 330]])
w = np.array([0.5, 0.2, 0.3])  # trust in agent, tax, online

The sources are trusted unequally, so a listing's consensus is a weighted mean. np.average(a, axis=1, weights=w) computes one for every row: each value times its weight, summed, divided by the sum of the weights.

1. Count the missing estimates per source and per listing.
2. Compute each listing's weighted consensus with np.average, and say which listing it fails on and why.
3. Repair it. Give every cell its source's weight, or 0 where the estimate is missing, and divide each listing's weighted sum by that listing's own total weight, so a listing with a missing source is averaged over the sources it has. Print the repaired consensus and confirm that the four complete listings kept their np.average values.
4. Center every estimate on its listing's repaired consensus, leaving the missing cell missing, and print the position and the estimates of the listing the sources value highest.

Hints:

• np.where(np.isnan(est5), 0, w) broadcasts the (3,) weights against the (5, 3) mask, giving every cell its source's weight, or 0 where the estimate is missing.
• nan times 0 is still nan, so the weighted sum needs np.nansum, and each row's divisor is the sum of that row's own weights.

import numpy as np

est5 = np.array([[280, 255, 300],
                 [350, 310, np.nan],
                 [199, 180, 215],
                 [460, 400, 480],
                 [315, 290, 330]])
w = np.array([0.5, 0.2, 0.3])

missing = np.isnan(est5)
print(missing.sum(axis=0), missing.sum(axis=1))

naive = np.average(est5, axis=1, weights=w)
print(naive)

[0 0 1] [0 1 0 0 0]
[281.    nan 200.  454.  314.5]

W = np.where(missing, 0, w)  # (5, 3), each cell's weight
consensus = np.nansum(est5 * W, axis=1) / W.sum(axis=1)
print(W.sum(axis=1))
print(consensus)

complete = missing.sum(axis=1) == 0
print(np.allclose(consensus[complete], naive[complete]))

[1.  0.7 1.  1.  1. ]
[281.         338.57142857 200.         454.         314.5       ]
True

centered = est5 - consensus[:, np.newaxis]
print(np.round(centered, 2))
top = consensus.argmax()
print(top, est5[top])

[[ -1.   -26.    19.  ]
 [ 11.43 -28.57    nan]
 [ -1.   -20.    15.  ]
 [  6.   -54.    26.  ]
 [  0.5  -24.5   15.5 ]]
3 [460. 400. 480.]

Additional Resources

Optional, and nothing later depends on them. Worth opening if you want the reason behind ddof=1, or the rest of NumPy's statistics functions.

• YouTube - "Why Dividing By N Underestimates the Variance" by StatQuest with Josh Starmer. Why a spread measured from a sample's own mean comes out too small, and what dividing by n minus 1 corrects.
• Official docs - NumPy reference, "Statistics", the "Order statistics" section. The range ptp and the nan versions of percentile and quantile, with the averages and variances, average and the nan versions of std and var among them, listed beneath.

Summary

Aggregation is where a table starts to say something. Totals, means, spreads, middles and cut points turn rows of numbers into sentences about them, and argmax turns "the most expensive" into a row you can read. None of it needed a new idea about arrays, only a choice about which dimension to give up.

That choice is the axis, and one sentence carries it: the axis you name is the axis that disappears. Read that way, axis=0 stops sounding like a direction and starts meaning "collapse the listings", and the result's shape confirms it before a single value is read. The shape cannot confirm intent, though, and a square table makes the wrong answer look right, so the habit from the previous lesson returns: decide what every number belongs to before writing the call.

Standardization ties the two lessons together. A mean and a spread collapsed down the rows broadcast straight back across them, and keepdims does the same service for summaries taken across a row. The quieter conventions at the end, a nan spreading through a total and a standard deviation divided by n, decide whether two people computing the same statistic agree, and the next section builds linear algebra on these arrays without explaining any of this again.

Cheat Sheet

Part A - Reference Table

Part B - Compact Code Block

import numpy as np

est = np.array([[310, 280, 325], [245, 231, 250]])
est.mean()                             # every element, one number
est.mean(axis=0)                       # per source, (3,)
est.mean(axis=1)                       # per listing, (2,)
est - est.mean(axis=1, keepdims=True)  # center each row

X = np.array([[85, 3, 12], [62, 2, 30], [120, 4, 5]])
Z = (X - X.mean(axis=0)) / X.std(axis=0)
np.round(Z.mean(axis=0), 3)            # zeros, maybe -0.

p = np.array([305, 248, np.nan, 200])
np.isnan(p).sum(), np.nanmean(p)       # count first, then skip
p[~np.isnan(p)].std(ddof=1)            # divide by n - 1