Reductions, axes, NaN and masked data

Five students sat four papers, and two of the twenty marks are missing. marks.sum(axis=1) gives NaN for both students who missed a paper, np.nansum ranks one of them below a weaker student, and np.nanmean ranks the same student first: three answers from one table, each correct for its own question.

You will reduce along the axis you mean, keep dimensions so a result lines up with its table, choose the divisor of a standard deviation, handle missing values with the nan-functions and masked arrays, and pick a percentile method on purpose.

An axis argument collapses that axis and keeps the others

The table has one row per student and one column per paper: maths, physics, chemistry and English. A missed paper is NaN. A reductionreduction: An operation that combines many values into fewer, such as sum, mean or max over an axis. The reduced axis disappears from the shape, unless keepdims=True keeps it with length 1. such as sum combines values along an axis and removes that axis from the result. Step through summing along each axis:

The run computes the same sums, the means and the best student each way:

"""Five students, four papers, two missed papers: reductions along each axis."""
import numpy as np

# Rows: students 0 to 4. Columns: maths, physics, chemistry, English. NaN: a paper the student missed.
marks = np.array([
    [78, 64, 81, 70],
    [92, 88, np.nan, 85],
    [55, 71, 60, 66],
    [83, 79, 90, 94],
    [61, 58, 72, np.nan],
])
print("shape", marks.shape)
print("marks.sum(axis=0), one per paper:  ", marks.sum(axis=0))
print("marks.sum(axis=1), one per student:", marks.sum(axis=1))
print("np.nansum(marks, axis=1):          ", np.nansum(marks, axis=1))
print("np.nanmean(marks, axis=1):         ", np.nanmean(marks, axis=1).round(2))
print("np.nanmean(marks, axis=0):         ", np.nanmean(marks, axis=0).round(2))
print("best student by total:", np.argmax(np.nansum(marks, axis=1)), "| by mean of the papers sat:", np.argmax(np.nanmean(marks, axis=1)))
print("best mark in each paper:", np.nanmax(marks, axis=0), "| by students", np.nanargmax(marks, axis=0))

shape (5, 4)
marks.sum(axis=0), one per paper:   [369. 360.  nan  nan]
marks.sum(axis=1), one per student: [293.  nan 252. 346.  nan]
np.nansum(marks, axis=1):           [293. 265. 252. 346. 191.]
np.nanmean(marks, axis=1):          [73.25 88.33 63.   86.5  63.67]
np.nanmean(marks, axis=0):          [73.8  72.   75.75 78.75]
best student by total: 3 | by mean of the papers sat: 1
best mark in each paper: [92. 88. 90. 94.] | by students [1 1 3 3]

marks.sum(axis=0) removed the students and left one value per paper, shape (4,); marks.sum(axis=1) removed the papers and left one per student, (5,); with no axis at all, every value would become one number. Because any sum with NaNNaN: Not a Number: the floating-point value for a missing or undefined result, such as 0 / 0. Arithmetic with NaN gives NaN, and NaN is not even equal to itself, so np.isnan is how to find it. in it is NaN, both students who missed a paper got NaN totals, and so did the two papers someone missed. np.nansum skips the NaN and np.nanmean divides by the number of marks present, so student 1 has a total of 265 but a mean of 88.33. np.argmax of the totals picked student 3, the one with 346; np.argmax of the means picked student 1. np.nanmax and np.nanargmax along axis 0 gave each paper’s best mark and who scored it.

keepdims keeps the reduced axis, so the result broadcasts back

A reduced result is often combined with the table it came from: each mark as a share of its student’s total, or each paper’s marks centred on that paper’s mean. keepdims=True keeps the reduced axis with length 1, and NumPy’s reference describes what that buys: the result “will broadcast correctly against the original array”. The run uses it both ways, and shows what goes wrong without it:

"""keepdims=True keeps the reduced axis with length 1, so the result lines up with the table again."""
import numpy as np

marks = np.array([[78, 64, 81, 70], [92, 88, np.nan, 85], [55, 71, 60, 66], [83, 79, 90, 94], [61, 58, 72, np.nan]])
totals = np.nansum(marks, axis=1, keepdims=True)
print("totals with keepdims=True:", totals.shape, "| without:", np.nansum(marks, axis=1).shape)
print("each mark as a share of its student's total:")
print((marks / totals).round(3))
try:
    marks / np.nansum(marks, axis=1)
except ValueError as error:
    print("without keepdims: ValueError:", error)
centred = marks - np.nanmean(marks, axis=0, keepdims=True)  # each paper's mean subtracted from its column
print("centred on each paper's mean, student 0:", centred[0].round(2))
print("np.nansum(marks, axis=(0, 1)):", np.nansum(marks, axis=(0, 1)), "| np.nansum(marks):", np.nansum(marks))
square = np.array([[1.0, 3.0], [2.0, 6.0]])  # a square table: the shapes fit either way
print("square / row totals, no keepdims:", (square / square.sum(axis=1)).tolist(), "(each column divided by a row's total)")
print("square / row totals, keepdims:   ", (square / square.sum(axis=1, keepdims=True)).tolist())

totals with keepdims=True: (5, 1) | without: (5,)
each mark as a share of its student's total:
[[0.266 0.218 0.276 0.239]
 [0.347 0.332   nan 0.321]
 [0.218 0.282 0.238 0.262]
 [0.24  0.228 0.26  0.272]
 [0.319 0.304 0.377   nan]]
without keepdims: ValueError: operands could not be broadcast together with shapes (5,4) (5,) 
centred on each paper's mean, student 0: [ 4.2  -8.    5.25 -8.75]
np.nansum(marks, axis=(0, 1)): 1347.0 | np.nansum(marks): 1347.0
square / row totals, no keepdims: [[0.25, 0.375], [0.5, 0.75]] (each column divided by a row's total)
square / row totals, keepdims:    [[0.25, 0.75], [0.25, 0.75]]

With keepdims=True the totals have shape (5, 1), so marks / totals divides each row by its own total: 78 / 293 = 0.266 for student 0’s maths. Without it the totals have shape (5,), NumPy lines them up with the last axis, the 4 papers, and the division raised a ValueError. np.nanmean(marks, axis=0, keepdims=True) has shape (1, 4), which lines up with the columns: student 0’s maths mark is 4.2 above the maths mean. axis may also be a tuple: np.nansum with axis=(0, 1) reduced both axes at once and gave 1,347, the sum of all 18 marks.

std divides by n minus ddof

The standard deviation measures spread: how far the values lie from their mean, on average, in the units of the data. NumPy adds up the squared distances from the mean, divides by N − ddof and takes the square root:

In words: the square root of the squared deviations added up and divided by N minus ddofddof: Delta degrees of freedom: std and var divide the sum of squared deviations by N minus ddof. ddof=0, NumPy's default, gives the population value; ddof=1 gives the sample value.. For 1, 2, 3 and 4 the mean is 2.5 and the squared deviations add up to 5, so dividing by 4 gives 1.25 and by 3 gives 1.667, whose square roots are 1.118 and 1.291:

"""std divides by n - ddof: ddof=0 treats the data as the whole population, ddof=1 as a sample of it."""
import numpy as np

x = np.array([1, 2, 3, 4])
print("np.std(x):", np.std(x), "| np.std(x, ddof=1):", np.std(x, ddof=1))
print("np.var(x):", np.var(x), "| np.var(x, ddof=1):", np.var(x, ddof=1))
squares = ((x - x.mean()) ** 2).sum()  # the squared distances from the mean, added up
print("sum of squared deviations:", squares, "| divided by 4:", squares / 4, "| by 3:", squares / 3)
print("correction=1 is the same as ddof=1:", np.std(x, correction=1) == np.std(x, ddof=1))
maths = np.array([78, 92, 55, 83, 61])
print("maths marks: np.std", np.std(maths).round(3), "| ddof=1", np.std(maths, ddof=1).round(3))

np.std(x): 1.118033988749895 | np.std(x, ddof=1): 1.2909944487358056
np.var(x): 1.25 | np.var(x, ddof=1): 1.6666666666666667
sum of squared deviations: 5.0 | divided by 4: 1.25 | by 3: 1.6666666666666667
correction=1 is the same as ddof=1: True
maths marks: np.std 13.79 | ddof=1 15.418

The default, ddof=0, treats the data as the whole population, and NumPy’s reference notes that statistics calls the result the “population standard deviation”. ddof=1 gives what it calls the “sample standard deviation”, for data that is a sample of something larger, such as the reaction times of 30 volunteers standing in for everyone. NumPy 2.0 added correction=, the array API name for the same argument. With few values the choice is visible: the five maths marks have a spread of 13.79 with ddof=0 and 15.418 with ddof=1.

NaN spreads through a reduction unless you skip it

Ordinary reductions let NaN spread, the nan-functions skip it, and data that is all NaN forces a decision. The run shows each case, with the warnings it caught:

"""NaN spreads through ordinary reductions, the nan-functions skip it, and all-NaN data needs a decision."""
import warnings

import numpy as np

v = np.array([1.0, np.nan, 3.0])
print("np.max:", np.max(v), "| np.argmax:", np.argmax(v), "| np.nanmax:", np.nanmax(v), "| np.nanargmax:", np.nanargmax(v), "| np.fmax.reduce:", np.fmax.reduce(v))
print("np.mean:", np.mean(v), "| np.nanmean:", np.nanmean(v), "| np.nansum:", np.nansum(v))


def report(label, fn):
    with warnings.catch_warnings(record=True) as caught:
        warnings.simplefilter("always")
        try:
            value = fn()
        except ValueError as error:
            value = f"ValueError: {error}"
    seen = "; ".join(f"{w.category.__name__}: {w.message}" for w in caught) or "no warning"
    print(f"{label:<26} {value!s:<38} {seen}")


all_nan = np.array([np.nan, np.nan])
report("np.nanmean(all_nan)", lambda: np.nanmean(all_nan))
report("np.nanmax(all_nan)", lambda: np.nanmax(all_nan))
report("np.nanmedian(all_nan)", lambda: np.nanmedian(all_nan))
report("np.nansum(all_nan)", lambda: np.nansum(all_nan))
report("np.nanargmax(all_nan)", lambda: np.nanargmax(all_nan))
report("np.nanstd([5.0], ddof=1)", lambda: np.nanstd(np.array([5.0]), ddof=1))

np.max: nan | np.argmax: 1 | np.nanmax: 3.0 | np.nanargmax: 2 | np.fmax.reduce: 3.0
np.mean: nan | np.nanmean: 2.0 | np.nansum: 4.0
np.nanmean(all_nan)        nan                                    RuntimeWarning: Mean of empty slice
np.nanmax(all_nan)         nan                                    RuntimeWarning: All-NaN slice encountered
np.nanmedian(all_nan)      nan                                    RuntimeWarning: All-NaN slice encountered
np.nansum(all_nan)         0.0                                    no warning
np.nanargmax(all_nan)      ValueError: All-NaN slice encountered  no warning
np.nanstd([5.0], ddof=1)   nan                                    RuntimeWarning: Degrees of freedom <= 0 for slice.

On [1, nan, 3], np.max returned nan and np.argmax pointed at the NaN, position 1, while np.nanmax gave 3.0 and np.nanargmax position 2; np.fmax.reduce, a ufunc method as in the previous lesson, skips NaN too. On data with nothing but NaN there is nothing to skip to, and the functions disagree about what to do:

NumPy’s reference for nanmean states the first row: for all-NaN slices, “NaN is returned and a RuntimeWarning is raised”. A sample standard deviation of one value hits a similar wall: np.nanstd([5.0], ddof=1) divides by 1 − 1 = 0 and returns nan with the warning Degrees of freedom <= 0 for slice.

Percentiles depend on the method you choose

A percentilepercentile: The value below which a given percentage of the data lies: the 25th percentile has a quarter of the values below it. With few values the answer depends on the method used between neighbouring values. is the value below which a given share of the data lies, and with few values there is more than one sensible answer. np.percentile offers 13 methods: the nine of the paper by Hyndman and Fan that its reference cites, and four NumPy variants of the default, "linear". Step through the default on four values:

The run computes six of the methods, works the default out by hand, and applies it per paper:

"""The 25th percentile of 1, 2, 3 and 4 by six of NumPy's methods, and the default one worked by hand."""
import numpy as np

x = np.array([1, 2, 3, 4])
for method in ["linear", "lower", "higher", "nearest", "midpoint", "median_unbiased"]:
    print(f"{method:<16} {float(np.percentile(x, 25, method=method)):.4f}")
position = (len(x) - 1) * 25 / 100  # where the linear method reads the sorted values: 0.75
below = int(np.floor(position))
print("by hand: position", position, "->", x[below] + (position - below) * (x[below + 1] - x[below]))
print("np.quantile(x, 0.25):", np.quantile(x, 0.25), "| np.median(x):", np.median(x))
try:
    np.percentile(x, 25, interpolation="lower")
except TypeError as error:
    print("interpolation=:", error)
marks = np.array([[78, 64, 81, 70], [92, 88, np.nan, 85], [55, 71, 60, 66], [83, 79, 90, 94], [61, 58, 72, np.nan]])
print("np.nanpercentile(marks, 50, axis=0):", np.nanpercentile(marks, 50, axis=0))

linear           1.7500
lower            1.0000
higher           2.0000
nearest          2.0000
midpoint         1.5000
median_unbiased  1.4167
by hand: position 0.75 -> 1.75
np.quantile(x, 0.25): 1.75 | np.median(x): 2.5
interpolation=: percentile() got an unexpected keyword argument 'interpolation'
np.nanpercentile(marks, 50, axis=0): [78.  71.  76.5 77.5]

The default reads the sorted values at position (N − 1) × q, here 3 × 0.25 = 0.75, three quarters of the way from 1 to 2: 1.75, which the hand calculation matched. lower and higher take the neighbour below or above, 1 and 2, nearest the closer one, 2, midpoint the half-way value, 1.5, and median_unbiased uses a different position rule and gives 1.4167. np.quantile(x, 0.25) is the same as np.percentile(x, 25) with q as a fraction. The old keyword interpolation= raised a TypeError: NumPy 2.4 removed it, and method= replaces it. np.nanpercentile skips NaN, so the median chemistry mark comes from the four marks present, 76.5.

Masked arrays carry the missing positions with the data

A masked arraymasked array: An array from numpy.ma that carries a mask of True and False beside its data. Masked positions, True in the mask, are left out of computations such as mean and count. from numpy.ma keeps a mask of True and False beside the data, and NumPy’s guide states the rule it follows: masked entries are not used in computations. np.ma.masked_invalid masks every NaN and inf:

"""A masked array carries the missing positions along with the data."""
import numpy as np

marks = np.array([[78, 64, 81, 70], [92, 88, np.nan, 85], [55, 71, 60, 66], [83, 79, 90, 94], [61, 58, 72, np.nan]])
m = np.ma.masked_invalid(marks)  # NaN and inf become masked positions
print(m)
print("mask of student 1:", m.mask[1])
print("mean per paper:", m.mean(axis=0))
print("marks counted per paper:", m.count(axis=0))
print("student 1 filled with 0:", m.filled(0)[1], "| fill_value:", m.fill_value)
print("inf is masked as well:", np.ma.masked_invalid(np.array([1.0, np.inf, np.nan, 4.0])).mean())

[[78.0 64.0 81.0 70.0]
 [92.0 88.0 -- 85.0]
 [55.0 71.0 60.0 66.0]
 [83.0 79.0 90.0 94.0]
 [61.0 58.0 72.0 --]]
mask of student 1: [False False  True False]
mean per paper: [73.8 72.0 75.75 78.75]
marks counted per paper: [5 5 4 4]
student 1 filled with 0: [92. 88.  0. 85.] | fill_value: 1e+20
inf is masked as well: 2.5

The printout shows the masked marks as --. m.mean(axis=0) skipped them and gave the same means as np.nanmean, m.count(axis=0) counted the marks present per paper, 5, 5, 4 and 4, and m.filled(0) made an ordinary array with 0 in the masked places, for code that does not understand masks; the default fill_value of a float array is 1e+20. Unlike NaN, a mask can mark any value as missing, an inf or a placeholder such as −1, without changing it. For tables with named rows and columns, a labelled table library such as pandas is usually the better tool.

Practise: percentiles and column summaries by hand

Two functions in plain Python make you decide every case this lesson raised: the default percentile, worked out from the sorted values, and a per-column summary that skips missing values and returns NaN, without any warning, where there is too little data. Their sample tests run here:

def percentile_linear(values, q):
    """Return the q-th percentile of values (q from 0 to 100) as np.percentile's default "linear" method gives it."""
    ordered = sorted(values)
    return float(ordered[round((len(ordered) - 1) * q / 100)])

import math

from percentile import percentile_linear


def test_the_lessons_quartile():
    """gives 1.75 for the 25th percentile of 1, 2, 3 and 4, in any input order"""
    assert percentile_linear([1, 2, 3, 4], 25) == 1.75
    assert percentile_linear([4, 1, 3, 2], 25) == 1.75


def test_the_ends_and_the_median():
    """gives the smallest value at 0, the largest at 100 and the middle at 50"""
    assert percentile_linear([1, 2, 3, 4], 0) == 1.0
    assert percentile_linear([1, 2, 3, 4], 100) == 4.0
    assert percentile_linear([1, 2, 3, 4], 50) == 2.5
    assert percentile_linear([7], 90) == 7.0


def test_the_maths_marks():
    """matches np.percentile on the lesson's maths marks"""
    marks = [78, 92, 55, 83, 61]
    assert math.isclose(percentile_linear(marks, 10), 57.4)
    assert percentile_linear(marks, 50) == 78.0
    assert percentile_linear(marks, 75) == 83.0


def test_bad_input():
    """refuses an empty list and a q outside 0 to 100"""
    for values, q in (([], 50), ([1, 2], 101), ([1, 2], -1)):
        try:
            percentile_linear(values, q)
        except ValueError:
            continue
        assert False, f"{values!r} with q = {q} must raise ValueError"

import math


def column_summary(rows):
    """Return (count, mean, std with ddof=1, median) of each column, skipping NaN, as the nan-functions give them."""
    result = []
    for column in zip(*rows):
        mean = sum(column) / len(column)
        result.append((len(column), mean, 0.0, sorted(column)[len(column) // 2]))
    return result

import math

from summary import column_summary

nan = math.nan
MARKS = [
    [78, 64, 81, 70],
    [92, 88, nan, 85],
    [55, 71, 60, 66],
    [83, 79, 90, 94],
    [61, 58, 72, nan],
]


def close(got, expected):
    """True when two numbers agree to 1e-9, counting NaN as equal to NaN."""
    if math.isnan(expected):
        return math.isnan(got)
    return math.isclose(got, expected, rel_tol=1e-9, abs_tol=1e-12)


def test_the_lessons_marks():
    """summarises maths and chemistry as np.nanmean, np.nanstd(ddof=1) and np.nanmedian do"""
    maths, _, chemistry, english = column_summary(MARKS)
    assert maths[0] == 5 and close(maths[1], 73.8) and close(maths[2], 15.417522498767434) and close(maths[3], 78.0)
    assert chemistry[0] == 4 and close(chemistry[1], 75.75) and close(chemistry[3], 76.5)
    assert close(chemistry[2], 12.816005617976296)
    assert english[0] == 4 and close(english[3], 77.5)


def test_one_value_and_no_value():
    """gives NaN for the spread of one value, and NaN everywhere for a column with no value"""
    only = column_summary([[5.0], [nan]])[0]
    assert only[0] == 1 and close(only[1], 5.0) and math.isnan(only[2]) and close(only[3], 5.0)
    empty = column_summary([[nan], [nan]])[0]
    assert empty[0] == 0 and math.isnan(empty[1]) and math.isnan(empty[2]) and math.isnan(empty[3])

The self-check revisits all six ideas, beginning where your answers so far are least sure:

Interview questions

Warm-up (fresher): For an array of shape (5, 4), what do sum(axis=0) and sum(axis=1) return, and what does keepdims=True change?

A reduction removes the axis you name. sum(axis=0) adds down the first axis, the 5 rows, and returns one value per column, shape (4,); sum(axis=1) adds across each row and returns one value per row, shape (5,). With no axis, the whole array becomes one number. keepdims=True keeps the reduced axis with length 1, so the results have shapes (1, 4) and (5, 1), and they broadcast against the original array: a / a.sum(axis=1, keepdims=True) divides each row by its own total. Without keepdims, a (5,) result lines up with the last axis instead and raises an error here, or, on a square array, silently pairs the wrong values.

Two interview questions take these decisions to real reports, a ranking with missing papers and two percentiles that disagree; the module’s question page says where each answer starts:

Key takeaways

• A reduction removes the axis you name: axis=0 combines down the rows and axis=1 across them, and NaN anywhere in a slice makes an ordinary sum or mean NaN.
• keepdims=True keeps the reduced axis with length 1 so the result broadcasts back against its array; without it, a square table pairs the wrong values without an error.
• std and var divide by N − ddof: ddof=0, NumPy’s default, for a whole population and ddof=1 for a sample, while pandas defaults to ddof=1.
• The nan-functions skip NaN, but on all-NaN data nanmean and nanmax return nan with a warning, nansum returns 0.0 silently and nanargmax raises.
• np.percentile has 13 methods, and the default linear one interpolates at position (N − 1) × q, 1.75 for the 25th percentile of 1 to 4.
• A masked array carries the missing positions with the data, skips them in computations and fills them on request.