Instead of examining every single value in an array, a single summarizing number — the total, the average, the range — is often what you actually need. numpy statistical functions provide optimized, vectorized ways to compute exactly these summaries. This article covers numpy mean median std and the other core aggregates, with particular care given to the axis parameter — the single concept most competing guides explain poorly.
Summarizing Data Instead of Reading Every Value
The core need these functions serve
Raw data, by itself, is often too much to reason about directly — a thousand individual sensor readings, a spreadsheet of a hundred students' scores. What you usually need instead is a summary: the average, the middle value, how spread out the values are. NumPy provides optimized, vectorized functions specifically for these aggregations, built on exactly the same compiled-C-speed foundation covered in the earlier vectorization article.
What this article covers
Measures of central tendency (mean, median), measures of spread (variance, standard deviation), other common aggregates (sum, min/max, argmin/argmax), and — given genuine, careful attention — the axis parameter, which controls how every one of these functions behaves once your data has more than one dimension.
Measures of Central Tendency: Mean, Median, and Weighted Average
np.mean(): the arithmetic average
import numpy as np
scores = np.array([85, 90, 78, 92, 88])
print(np.mean(scores)) # 86.6np.mean() calculates the arithmetic mean — the sum of every element divided by the count of elements — along a specified axis, or, if no axis is given, across the entire flattened array, regardless of its original shape.
np.median(): the middle value
print(np.median(scores)) # 88.0The median is the value separating the higher half of sorted data from the lower half: the middle element, once the data is sorted, for an odd-length array — or the average of the two middle elements, for an even-length array.
even_scores = np.array([70, 75, 80, 85])
print(np.median(even_scores)) # 77.5 — average of 75 and 80, the two middle valuesThe median is genuinely useful as a complement to the mean, specifically because it's far less sensitive to outliers — a single extreme value can dramatically shift a mean, while barely affecting a median at all.
np.average(): a weighted mean
grades = np.array([85, 92, 78])
weights = np.array([0.5, 0.3, 0.2]) # exam worth 50%, midterm 30%, homework 20%
weighted_avg = np.average(grades, weights=weights)
print(weighted_avg) # 85.5np.average() supports an optional weights argument, letting different values contribute unequally to the result — genuinely important whenever a plain, unweighted mean would misrepresent the actual relative importance of different values, exactly as shown in this grade-weighting example.
A note on mode
NumPy itself has no built-in mode function (the most frequently occurring value) — this is a genuine gap worth knowing about, since it's easy to assume every basic statistic lives directly in NumPy. scipy.stats.mode(), from the separate SciPy library, fills exactly this gap:
from scipy import stats
data = np.array([1, 2, 2, 3, 3, 3, 4])
result = stats.mode(data)
print(result.mode) # 3Measures of Spread: Variance, Standard Deviation, and Range
np.var(): variance
scores = np.array([85, 90, 78, 92, 88])
print(np.var(scores)) # 24.24Variance is the average of the squared deviations from the mean — it captures how spread out the data is, but in squared units, which makes it somewhat awkward to interpret directly (variance of "scores" comes out in "squared score points," which isn't a naturally meaningful unit).
np.std(): standard deviation
print(np.std(scores)) # 4.923...Standard deviation is simply the square root of the variance, bringing the measure of spread back into the same units as the original data — genuinely more interpretable than raw variance, and the far more commonly reported statistic in practice. A larger standard deviation means the data points are, on average, further from the mean; a smaller one means they cluster more tightly around it.
np.ptp(): the range (peak-to-peak)
print(np.ptp(scores)) # 14 — 92 (max) minus 78 (min)np.ptp() (peak-to-peak) calculates the simplest possible spread measure: the maximum value minus the minimum. It's genuinely useful as a fast, rough sanity check on a dataset's overall scale — is this data ranging over a handful of units, or thousands? — before reaching for the more nuanced variance and standard deviation.
Understanding the axis Parameter
This is the concept that trips up more learners than anything else in this article, so it's worth working through concretely, with a genuinely representative example, rather than stated abstractly.
A grade matrix: the concrete scenario
Imagine a 2D array where each row is a student, and each column is an assignment:
grades = np.array([
[85, 90, 78], # Student 1's scores on Assignment 1, 2, 3
[92, 88, 95], # Student 2's scores
[70, 75, 80], # Student 3's scores
])Two genuinely different, equally valid questions could be asked of this exact same data: "what's the average score per assignment?" (aggregating down each column) versus "what's the average score per student?" (aggregating across each row). The axis parameter is precisely what tells NumPy which of these two directions you actually mean.
axis=0: aggregating down columns
column_means = grades.mean(axis=0)
print(column_means) # [82.33 84.33 84.33] — one average PER ASSIGNMENTaxis=0 means "collapse the rows, keeping the columns separate" — producing one output value per column: the average score on Assignment 1 across all three students, the average on Assignment 2, and the average on Assignment 3.
A column-sum walkthrough, step by step
column_sums = grades.sum(axis=0)
print(column_sums) # [247 253 253]Walking through exactly what happened for the first output value, 247: NumPy took the first column — [85, 92, 70], one value from each of the three rows — and added those three values together: 85 + 92 + 70 = 247. The same happens independently for each remaining column, producing one summed value per column, with the row dimension having been "collapsed" entirely.
axis=1: aggregating across rows
row_means = grades.mean(axis=1)
print(row_means) # [84.33 91.67 75. ] — one average PER STUDENTaxis=1 means "collapse the columns, keeping the rows separate" — producing one output value per row: each individual student's average score across their three assignments. The first value, 84.33, is the average of [85, 90, 78] — Student 1's own three scores.
The general rule worth remembering
The genuinely useful mental shortcut: axis=N means "the dimension at position N disappears from the result, having been collapsed down into a single aggregated value." For a standard 2D array, axis=0 collapses rows (aggregating down each column), and axis=1 collapses columns (aggregating across each row) — the exact opposite of what the number might intuitively suggest to some people at first, which is precisely why working through a concrete example like the grade matrix above is worth the effort.
The same logic applies to every aggregate function
print(grades.std(axis=0)) # standard deviation per assignment (down each column)
print(np.median(grades, axis=1)) # median score per student (across each row)Every aggregate function covered in this article — mean, median, std, var, sum, and the rest — accepts this exact same axis parameter, with exactly this same behavior. Once the grade-matrix mental model above genuinely clicks, it transfers directly to every one of these functions without needing to relearn the concept each time.
Beyond the Basics: argmin/argmax, Percentiles, and NaN-Aware Functions
argmin() and argmax(): finding the position, not the value
scores = np.array([85, 90, 78, 92, 88])
print(np.argmax(scores)) # 3 — the INDEX where the maximum value (92) occurs
print(np.argmin(scores)) # 2 — the INDEX where the minimum value (78) occursThis is worth being precise about: argmax()/argmin() don't return the maximum or minimum value itself (that's what plain .max()/.min() are for) — they return the position where that value occurs. This is genuinely useful whenever you need to know which record holds an extreme value, not just what that extreme value is:
grades = np.array([
[85, 90, 78],
[92, 88, 95],
[70, 75, 80],
])
best_student_per_assignment = np.argmax(grades, axis=0)
print(best_student_per_assignment) # [1 1 1] — student at index 1 scored highest on every assignmentnp.percentile(): arbitrary percentiles
data = np.array([10, 20, 30, 40, 50, 60, 70, 80, 90, 100])
print(np.percentile(data, 50)) # 55.0 — the 50th percentile, equivalent to the median
print(np.percentile(data, 90)) # 91.0 — the 90th percentilenp.percentile() generalizes the median (which is simply the 50th percentile) to any percentile you specify, and it supports the same axis parameter as every other aggregate function covered in this article, behaving exactly according to the grade-matrix mental model from Section 4.
np.corrcoef(): measuring the relationship between variables
hours_studied = np.array([1, 2, 3, 4, 5])
exam_scores = np.array([50, 60, 65, 80, 85])
correlati np.corrcoef(hours_studied, exam_scores)
print(correlation_matrix)
# [[1. 0.98...]
# [0.98... 1. ]]np.corrcoef() computes the correlation coefficient between two variables — a value close to 1 indicates a strong positive relationship (as one variable increases, so does the other), close to -1 indicates a strong negative relationship, and close to 0 indicates little to no linear relationship. The result is a full correlation matrix (since it also includes each variable's correlation with itself, always 1), with the actual correlation between your two variables sitting at position [0, 1] (or equivalently [1, 0]).
NaN-aware functions: handling missing data correctly
Real datasets frequently have gaps — missing sensor readings, incomplete survey responses — commonly represented as NaN (Not a Number) in NumPy. This creates a genuine, practical problem: a single NaN in an array poisons an ordinary aggregate calculation entirely.
data_with_gaps = np.array([85, 90, np.nan, 92, 88])
print(np.mean(data_with_gaps)) # nan — the ENTIRE result becomes NaN, even though only one value was missing!
print(np.nanmean(data_with_gaps)) # 88.75 — correctly ignores the NaN and averages the restNumPy's NaN-aware functions — np.nanmean(), np.nanstd(), np.nanmedian(), and others, each mirroring a standard aggregate function with a nan prefix — ignore NaN values entirely when computing along an axis, rather than letting a single missing value silently turn an entire row's, column's, or array's result into a useless NaN. This is genuinely essential for working with any real-world dataset that has gaps — reaching for the nan-prefixed version of whatever aggregate you need is a small habit that avoids a surprisingly common, easy-to-miss source of silently broken results.