NumPy Arrays and Vectorized Computing

NumPy stores numerical data in arrays so that one expression can work on many values. Import it with import numpy as np; the short name np is the usual way to call its functions.

Read an array's structure

An ndarray is a rectangular grid whose elements share one data type. A dimension is also called an axis. Start with a small table:

import numpy as np
a = np.array([[1, 2, 3], [4, 5, 6]])

Its attributes answer four different questions:

Axis 0 indexes rows; axis 1 indexes columns. A data table with one row per person and columns for height, mass, and age would have shape (100, 3) for 100 people. The product of the shape entries gives the size.

A vector has no row or column axis

np.array([2, 4, 6, 8]) has shape (4,) and one axis. The comma makes this a one-entry tuple. Shape (1, 4) has two axes: one row and four columns. Shape (4, 1) also has two axes, but four rows and one column. These arrays contain four values, yet behave differently when combined with other arrays. np.array(7) is a zero-dimensional array: its shape is () and its ndim is 0.

Arrays change the meaning of operators

A list is a general Python container. An array is designed for calculations on matching positions.

In the array addition, 1 pairs with 3 and 2 pairs with 4. The array stays the same size.

Why numerical arrays are efficient

A newly created numerical array stores same-typed values together in memory. NumPy runs the loop through these values in compiled code, avoiding a separate Python-level operation for every number. A million float64 values occupy 8,000,000 bytes of element storage, available as a.nbytes on an array of that size. A Python list also stores references to Python objects, so it needs substantially more space. Performance gains depend on the calculation and computer.

Array size does not grow through an in-place append method. np.append(a, value) allocates a new array. For data collected one value at a time, build a Python list and convert once. Nested input must be rectangular: np.array([[1, 2], [3]]) raises ValueError because its row lengths differ.

Recognize displayed results

print(np.array([1, 2, 3])) shows [1 2 3]. An interactive display uses array([1, 2, 3]). For np.array([5, 6])[0] on a 64-bit computer, NumPy 2 displays np.int64(5); it participates in arithmetic as the number 5.

Choose a creation function by asking what you know: the actual values, a shape and fill value, a step size, or a required number of points. That choice makes the code express how the data should be built.

Start from values or a shape

np.array([3, 1, 4]) creates a vector. A list of equal-length lists creates a table. Pass the values as one sequence; np.array(1, 2, 3) does not supply three elements correctly.

For a multidimensional shape, the parentheses matter. np.zeros((2, 3)) passes one tuple. In np.zeros(2, 3), Python sends 3 as the dtype argument, which causes a type error. Zeros and ones use float64 by default; full infers its dtype from the fill value. Add dtype=int when you need whole-number storage.

empty does not promise zeros or meaningful starting values. Read from it only after assigning every element you will use. zeros is a simpler choice when an initial zero has meaning.

Choose a step or choose a point count

np.arange(start, stop, step) follows the stepping idea of Python's range. np.arange(2, 11, 3) gives [2 5 8]: start at 2 and add 3 while staying below 11. With one argument, np.arange(5) starts at 0 and gives [0 1 2 3 4].

Floating-point steps can give surprising endpoints. For example, np.arange(1, 1.3, 0.1) includes the displayed value 1.3 because binary rounding affects the calculation of the sequence. For a decimal grid with a fixed point count, use linspace:

grid = np.linspace(0, 1, 5)
# array([0.  , 0.25, 0.5 , 0.75, 1.  ])

Both endpoints are included by default. Five points make four gaps. Ten points between 0 and 1 make nine gaps, so their spacing is about 0.111; steps of 0.1 need 11 points.

Reuse an existing structure

np.zeros_like(a), np.ones_like(a), and np.full_like(a, value) create new arrays with a's shape and dtype. The inherited dtype matters: filling an integer template with a decimal can lose the decimal part unless you request a floating dtype. np.eye(3) creates a square identity array with ones on the diagonal and zeros elsewhere.

Allocate a data table once

For rainfall in millimeters at three stations over four weeks:

rain = np.zeros((3, 4))
rain[0] = [12.0, 0.0, 5.5, 8.0]

The first assignment fills one row and preserves shape (3, 4). If the number of readings is known, reserving the full table avoids repeatedly allocating larger arrays. If the amount of input is unknown, accumulate values in a list and call np.array when collection is finished.

An array's dtype controls what it can store and how much memory each value needs. Choosing the right type also prevents calculations from quietly losing information.

Inspect the storage type

int64 holds 64-bit whole numbers, float64 holds binary floating-point numbers, and bool holds true/false values. Unicode dtypes such as <U5 store text with a maximum character length. a.itemsize is the bytes per element; a.nbytes is a.size * a.itemsize. A float64 uses 8 bytes and has about 15 significant decimal digits of precision.

On typical 64-bit computers, whole-number input uses 64-bit integers. NumPy 2 also uses that default on 64-bit Windows. Decimal input normally gives float64.

Inference chooses one shared type

np.array([1, 2.5]) becomes [1. 2.5]: the integer is represented as a float. np.array([True, 2]) becomes [1 2]. Mixing numbers and text can turn the entire array into text, so a label such as 'n/a' can make later arithmetic fail. Inspect dtype when apparently numerical data will not add or multiply.

Convert deliberately

astype returns a converted copy:

x = np.array([1.7, -1.7, 2.5])
whole = x.astype(int)
# array([ 1, -1,  2])

This conversion cuts toward zero; it does not round. np.round(x).astype(int) gives [2 -2 2] because exact halves round to the nearest even value. Text such as '1.5' converts with astype(float), but 'n/a' raises ValueError. Clean placeholders before conversion.

Assignments use the destination's dtype. After a = np.array([1, 2, 3]), assigning a[0] = 9.9 stores 9 without warning. Converting the entire array to float first lets later assignments retain decimals.

Fixed-width integers can overflow

An int8 holds values from −128 to 127; a uint8 holds 0 to 255. Adding two uint8 arrays containing 250 and 10 produces 4 through wraparound. Multiplying an int8 array containing 100 by 2 produces −56. A Python integer that cannot fit the array type, such as 300 added to an int8 array in NumPy 2, raises OverflowError instead. Widen small integer types before calculations that can exceed their range.

Division and approximate decimals

Floor division rounds down toward negative infinity, which differs from truncation for negative results. True division returns floating-point results for integer inputs.

Many decimal fractions are approximations in binary. 0.1 + 0.2 == 0.3 is false; np.isclose compares corresponding values with a tolerance, and np.allclose gives one whole-array answer.

The missing marker np.nan is floating-point. Integer arrays cannot store it. np.array([1, 2, np.nan]) therefore becomes floating-point; converting missing values to integers does not repair them and can yield meaningless results.

Indexing locates values by position. Slicing selects a regular run of positions. In a table, write one selector for each axis so the code says both which rows and which columns you need.

Read a vector

For a = np.array([10, 20, 30, 40, 50]), position 0 contains 10 and position −1 contains 50. A slice has the form start:stop:step, with the stop excluded:

Leaving out start or stop extends to the relevant end. A slice can run beyond the array and be shortened automatically. A single invalid position, such as a[5], raises IndexError.

Give each table axis its own selector

m = np.array([[1, 2, 3],
              [4, 5, 6],
              [7, 8, 9]])

m[1, 2] reads row 1, column 2, giving 6. The comma separates axes. A colon alone selects the full axis: m[0, :] gives the first row [1 2 3], and m[:, 1] gives the second column [2 5 8]. For a block, m[:2, 1:] selects rows 0 and 1 and columns 1 and 2, giving [[2 3] [5 6]].

Avoid treating repeated square-bracket indexing as interchangeable with comma indexing. m[1][2] happens to locate the same scalar as m[1, 2]. But m[:][0] first takes all rows and then selects the first row. It does not select the first column, which is m[:, 0].

Predict which axes remain

A single integer selector removes its axis. A slice retains its axis even when only one position is selected. Thus m[:, 1] has shape (3,), while m[:, 1:2] has shape (3, 1). Both select 2, 5, and 8, but only the second keeps a column axis. This difference becomes important when subtracting or multiplying arrays with other shapes.

The same rule extends beyond tables. If t.shape is (2, 3, 4), t[0] has shape (3, 4), t[0, 1] has shape (4,), and t[:, :, 0] has shape (2, 3). Any trailing axes omitted from the index are selected in full.

Change the selected positions

Assignment through an index writes into the existing array. a[1:3] = 0 replaces both selected values with zero; the scalar is applied throughout the selection. m[0] = [0, 0, 0] replaces a complete row. The assigned values must fit the selection's shape and are converted to its dtype.

For four students and three tests:

scores = np.array([[78, 85, 92], [64, 70, 58],
                   [88, 91, 95], [72, 68, 80]])

scores[2] is the third student's [88 91 95]. scores[:, -1] gives the final-test results [92 58 95 80]. scores[1:3, :2] gives [[64 70] [88 91]]: two students, each retaining two test columns.

A selection can share its data with the original array or own independent data. Knowing which one you have prevents edits to a small selection from unexpectedly changing a larger table.

Distinguish an alias, a view, and a copy

For a = np.arange(6) and b = a[2:5], assigning b[0] = 99 makes a equal to [0 1 99 3 4 5]. The first position of the view refers to position 2 in the original. This differs from slicing a Python list, which makes a new list.

Copy before editing when the original must remain intact. A small retained view can also keep a large original data buffer alive; copying that small piece lets it stand independently.

Select arbitrary positions

Fancy indexing uses an integer list or array. It allows reordered and repeated positions:

a = np.array([10, 20, 30, 40, 50])
a[[4, 0, 0]]
# array([50, 10, 10])

For m = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]), m[[0, 2]] selects the first and third rows. m[:, [2, 0]] selects the last and first columns in that order, giving [[3 1] [6 4] [9 7]].

Two integer lists in one index select paired coordinates. m[[0, 2], [1, 2]] pairs row 0 with column 1 and row 2 with column 2, producing [2 9]. It does not take all combinations. To get the block formed by both rows and both columns, use:

m[np.ix_([0, 2], [1, 2])]
# array([[2, 3], [8, 9]])

Selecting the rows first and then the columns, m[[0, 2]][:, [1, 2]], gives the same block.

Reading and assigning behave differently

Fancy indexing and boolean filtering return copies. If c = a[[0, 1]], changing c[0] leaves a unchanged. But a[[0, 1]] = -1 writes directly into those original positions. A saved selection and an assignment through an index are different operations.

This also explains why chained edits can fail to reach the source. In a[[0, 1]][0] = -1, the first selection is already a temporary copy. The final write changes that temporary array, not a.

Check sharing directly

np.shares_memory(x, y) checks whether two arrays overlap in memory. It returns true for an ordinary nonempty slice and false for a fancy-indexed copy. The base attribute can identify a backing owner, but it need not identify the array you most recently sliced. If a = np.arange(6), b = a[2:5], and c = b[0:2], then c.base is a is true and c.base is b is false. Use shares_memory when the question is specifically whether two arrays share data.

A boolean mask records which values satisfy a condition. Once a mask exists, the same true/false decisions can select values, count events, or choose replacements.

Build and apply a mask

t = np.array([18.5, 22.0, 15.2, 25.8, 21.1])
warm = t > 20
t[warm]
# array([22. , 25.8, 21.1])

The comparison produces [False True False True True]. Boolean indexing keeps the true positions and returns a copy. A 1-D mask used on an axis must have the same length as that axis.

To combine conditions, use & for and, | for or, and ~ for not. Put parentheses around each comparison:

comfortable = (t > 18) & (t < 23)
t[comfortable]
# array([18.5, 22. , 21.1])

Python's and and or try to get one truth value from an entire array and fail for multi-element arrays. Omitting the parentheses also changes how Python groups the expression because bitwise operators bind before comparisons.

Count, measure, or test conditions

True contributes 1 and False contributes 0. warm.sum() is 3 and warm.mean() is 0.6, the fraction of days above 20 °C. np.count_nonzero(warm) also counts true entries. warm.any() asks whether at least one day is warm; warm.all() asks whether every day is warm.

These operations answer different questions from t[warm].mean(), which calculates the average temperature among warm days. The mean of a mask measures a share; the mean of selected measurements retains the measurement's unit.

Replace values or choose between branches

Assignment through a mask updates the selected positions. With t2 = t.copy(), the line t2[t2 < 16] = 16 raises low readings to 16. np.clip(t, 16, 25) returns a new array bounded at both ends, [18.5 22. 16. 25. 21.1].

np.where(condition, x, y) builds a new result taking x at true positions and y at false positions. np.where(t > 20, 'warm', 'cool') supplies a label at every original position. In contrast, filtering shortens the result. With a condition alone, np.where(t > 20) returns a tuple of index arrays, here (array([1, 3, 4]),).

Keep values or keep complete records

A same-shape mask on a 2-D array extracts matching elements into one flat vector. In m[m > 5], the row structure is lost. To retain whole records, derive one boolean per row from a column:

scores = np.array([[78, 85, 92], [64, 70, 58],
                   [88, 91, 95], [72, 68, 80]])
scores[scores[:, 2] >= 90]
# array([[78, 85, 92], [88, 91, 95]])

Find missing values explicitly

NaN does not compare equal even to itself. Equality and ordering comparisons with NaN are false; inequality is true. Thus neither x == np.nan nor x != np.nan identifies present and missing entries correctly. np.isnan(x) marks missing positions, and x[~np.isnan(x)] keeps present values.

Changing shape changes how values are grouped into axes. Reshaping, transposing, and adding an axis do different jobs, even when some of them produce the same final shape.

Reshape without changing the element count

a = np.arange(6)
r = a.reshape(2, 3)
# array([[0, 1, 2], [3, 4, 5]])

The default order fills one row at a time: the last axis changes fastest. Reshaping the same values to (3, 2) gives [[0 1] [2 3] [4 5]]. The product of the new dimensions must equal the original size. Twelve elements fit (3, 4) but not (5, 3).

One dimension can be -1, letting NumPy infer its length. np.arange(12).reshape(-1, 3) gives four rows of three; reshape(2, -1) gives two rows of six. Two unknown dimensions are not allowed because their individual lengths cannot be determined from the total alone.

Transpose swaps the axes

For r above, r.T is [[0 3] [1 4] [2 5]], with shape (3, 2). Each old column is now a row. That arrangement differs from r.reshape(3, 2), even though the resulting shapes match.

A transpose is useful when the next operation needs the opposite axis order. A reshape is useful when a sequence already has a meaningful grouping, such as months grouped into quarters.

A 1-D transpose does not make a column

v = np.arange(3) has only one axis, so v.T still has shape (3,). Add an axis explicitly:

column = v[:, np.newaxis]  # shape (3, 1)
row = v[np.newaxis, :]     # shape (1, 3)

v.reshape(-1, 1) and np.expand_dims(v, 1) also create a column. np.squeeze(column) removes its length-one axis and returns shape (3,). Adding an axis changes no values but tells later arithmetic how they should line up.

Flattening and memory sharing

r.ravel() and r.flatten() both list values in a 1-D result. ravel returns a view when possible; flatten always copies. reshape(-1) is another way to request a 1-D arrangement, but whether it shares memory depends on the input's layout. Reshape may need a copy; transpose shares data in its array view. Use .copy() when independent storage is required.

Group monthly values into quarters

For monthly = np.array([120, 135, 150, 160, 155, 170, 180, 175, 165, 150, 140, 130]), monthly.reshape(4, 3) makes four quarter rows, each with three months. Summing each row with .sum(axis=1) gives [405 485 520 420]. The first total is 120 + 135 + 150 = 405. The reshape works because the chronological input already puts each quarter's months together.

Joining combines values from several arrays. Splitting cuts one array into pieces. First decide whether the operation should extend an existing axis or introduce a new one.

Concatenate along an existing axis

np.concatenate([a, b], axis=0) joins vectors end to end. Joining [1 2] with [3 4 5] gives [1 2 3 4 5].

For two arrays each shaped (2, 3), axis 0 adds rows and gives (4, 3). Axis 1 adds columns and gives (2, 6). Only the joining axis can differ in length. Thus shapes (2, 3) and (2, 2) can join along axis 1 into (2, 5), but cannot join along axis 0 because their column counts disagree.

The arrays must be passed together as a sequence, such as [a, b]. Joining creates a new combined array and copies the inputs' data.

Stack along a new axis

For p = np.array([1, 2, 3]) and q = np.array([4, 5, 6]):

stack requires equal input shapes and adds an axis. concatenate preserves the number of axes. The helper functions express common arrangements: vstack([p, q]) makes rows, hstack([p, q]) joins these 1-D inputs end to end, and column_stack([p, q]) makes columns. The 1-D behavior of hstack is why its name alone is not enough to predict a table.

Cut by boundaries or by part count

np.split(np.arange(10), [3, 7]) cuts before positions 3 and 7, yielding [0 1 2], [3 4 5 6], and [7 8 9]. The boundary positions are not part lengths.

An integer second argument instead requests equal parts: np.split(x, 3) asks for three equal pieces. Ten elements cannot be divided equally into three integer-sized pieces, so this raises ValueError. np.array_split(np.arange(10), 3) allows unequal lengths and produces 4, 3, and 3. These splitting calls return a collection of arrays, not one array with a new stacking axis.

Build a measurement table

h = np.array([1.62, 1.75, 1.80])  # meters
w = np.array([55.0, 72.5, 80.0])  # kilograms
data = np.column_stack([h, w])

Shape (3, 2) means three people with height and mass in columns. Selecting data[:, 1] / data[:, 0] ** 2 calculates mass divided by squared height for each row, giving approximately 21.0, 23.7, and 24.7 kg/m². Keeping the two measurement arrays in the same person order is essential.

If input arrives in batches, collect the batch arrays in a Python list and call concatenate once. Rejoining a growing array during every loop pass repeatedly copies its existing data.

Vectorization expresses a calculation for entire arrays. NumPy handles the repeated element-level work in compiled code, so the Python expression can stay close to the calculation you mean.

Translate a loop into matching-position arithmetic

prices = np.array([2.50, 4.00, 1.20])
qty = np.array([3, 1, 5])
line_totals = prices * qty
order_total = line_totals.sum()

The products are [7.5 4. 6.] and their sum is 17.5. Each price pairs with its own quantity. A Python loop updating a total produces the same result, but makes Python handle each pair separately.

All the arithmetic operators act element by element: +, -, *, /, //, %, and **. A scalar applies to every entry. 10 - np.array([1, 2, 3]) gives [9 8 7]; operand order still matters. Two arrays need equal shapes or compatible broadcasting shapes.

Keep element-wise and matrix products distinct

* multiplies corresponding entries. For A = np.array([[1, 2], [3, 4]]) and B = np.array([[5, 6], [7, 8]]), A * B gives [[5 12] [21 32]]. Matrix multiplication uses @ and would give [[19 22] [43 50]]. Choose * for paired measurements such as price and quantity.

Comparisons also return one answer per position

a == b produces a boolean array. It does not ask whether two arrays are equal as whole objects. An if statement needs one truth value, so if a == b: fails for multi-element results. np.array_equal(a, b) tests exact whole-array equality. For approximate decimals, np.allclose(a, b) tests agreement within tolerance.

Understand division warnings

With default NumPy error settings, floating-point division by zero produces a warning and special values:

np.array([1.0, -1.0, 0.0]) / 0
# array([ inf, -inf, nan])

np.isfinite identifies results that are neither infinite nor NaN. Integer floor division by zero instead returns integer zeros with a warning. Those zeros look finite, so inspect integer denominators with (d == 0).any() before floor division. A warning deserves attention even when the program continues.

Decide whether to update in place

a += value modifies existing storage and retains its dtype. If a is an integer array, a += 0.5 raises a casting error because the floating-point result cannot be stored under that operation's casting rule. a = a + 0.5 creates a floating-point result and binds a to it. This differs from assigning one decimal into an integer element, which can silently truncate.

Build a multi-step expression

For old = np.array([2.50, 4.00, 8.00]) and new = np.array([3.00, 3.00, 8.00]), (new - old) / old * 100 gives [20. -25. 0.]. First find each change, then divide by its original value, then express the ratio as a percentage. The signs show a rise, a fall, and no change. The original values must be nonzero for this percentage calculation.

A universal function, or ufunc, applies a calculation separately at each array position. These functions extend vectorization beyond ordinary arithmetic to roots, logarithms, rounding, comparisons, and other numerical operations.

Choose an array-aware function

Unary ufuncs take one input. np.sqrt(np.array([1, 4, 9])) gives [1. 2. 3.], and np.abs(np.array([-2, 3])) gives [2 3]. np.exp is the exponential function, np.log is the natural logarithm, and np.log10 is the base-10 logarithm. For example, np.exp(0) is 1.0 and np.log10(1000) is 3.0.

Python's math module works on individual numbers. math.sqrt(np.array([1, 4])) cannot calculate both roots and raises a type error. Use NumPy's version for numerical arrays. Operators also call ufuncs: a + b corresponds to np.add(a, b) and a * b to np.multiply(a, b).

Compare corresponding values or reduce a whole array

Binary ufuncs take two inputs. np.maximum([1, 5, 3], [4, 2, 6]) gives [4 5 6], choosing the larger value at every position. np.minimum chooses the smaller, and np.power raises corresponding bases to powers.

np.maximum(x, 0) sets negative entries to zero in a new result. It is different from np.max(x), an aggregation that returns one largest value. For [1, 5, 3], max gives 5; maximum needs another input and retains the element-wise structure.

Round with the intended rule

np.round(x, 2) keeps two decimal places. Exact halves follow the even rule: np.round([0.5, 1.5, 2.5]) gives [0. 2. 2.]. Binary approximation can affect decimal-looking ties, so interpret displayed decimals with their floating-point representation in mind.

Inspect values outside a function's domain

For real floating-point inputs, np.sqrt(-1.0) returns NaN with a warning, while np.log(0.0) returns negative infinity with a warning under default error settings. Unlike math.sqrt(-1), these calls do not normally halt with a domain exception. Check the inputs when the mathematical operation requires a restricted range. np.isnan marks NaN; np.isfinite marks values that are neither NaN nor infinite; np.isclose compares values with a tolerance.

Combine ufuncs into a useful calculation

For coordinate arrays x = np.array([3, 5, 8]) and y = np.array([4, 12, 15]), calculate distance from the origin with np.sqrt(x ** 2 + y ** 2). Squaring gives [9 25 64] and [16 144 225]; adding gives [25 169 289]; the roots are [5. 13. 17.]. np.hypot(x, y) calculates those distances directly.

Broadcasting lets arrays with different shapes participate in element-wise calculations. The aligned axes determine which values pair and which repeat.

Compare dimensions from the right

Line up shapes at their right-hand ends. Each pair of axis lengths must be equal or have one length equal to 1. Missing leading axes count as length 1. A length-one axis expands conceptually to match the other length; NumPy need not copy the repeated input values.

For (4, 3) and (3,), the last lengths match at 3. The missing leading dimension acts as 1 and expands to 4. The result is (4, 3): the vector applies to every row. (4, 3) and (4,) fail because the trailing lengths 3 and 4 disagree and neither is 1. An incompatible pair raises ValueError.

Center or standardize columns

X = np.array([[1, 10], [2, 20], [3, 30]])
centered = X - X.mean(axis=0)

The column means are [2. 20.], shape (2,). They line up with the last axis of (3, 2), giving [[-1. -10.] [0. 0.] [1. 10.]]. Dividing the centered result by X.std(axis=0) standardizes each column, producing approximately −1.22, 0, and 1.22 down both columns. Each column's spread must be nonzero for this division.

Keep an axis for row calculations

Use a score table with four students and three tests:

scores = np.array([[78, 85, 92], [64, 70, 58],
                   [88, 91, 95], [72, 68, 80]])

scores.mean(axis=1) has shape (4,). It cannot broadcast as one mean per row against (4, 3). Preserve the removed axis:

row_mean = scores.mean(axis=1, keepdims=True)
centered = scores - row_mean

row_mean now has shape (4, 1). Its leading 4 matches the rows; its trailing 1 expands across three columns. scores.mean(axis=1)[:, np.newaxis] creates the same shape. For the row [78, 85, 92], subtracting its mean of 85 gives [-7. 0. 7.].

Use the same rule to build a table

np.arange(1, 4)[:, np.newaxis] * np.arange(1, 5) combines (3, 1) and (4,). Both arrays expand to the result shape (3, 4), producing [[1 2 3 4] [2 4 6 8] [3 6 9 12]].

Expansion can also be accidental. A (3, 1) column minus a (3,) vector creates all nine pairwise differences, shape (3, 3). If you intended three corresponding differences, make both inputs shape (3,), for example with column.ravel().

Weight each column before adding

For the same score table and w = np.array([0.2, 0.3, 0.5]), scores * w applies the three weights to every row. Summing with .sum(axis=1) gives [87.1 62.8 92.4 74.8]. The first result is 15.6 + 25.5 + 46.0 = 87.1.

An aggregation combines values into fewer results: a total, mean, extreme, or measure of spread. Choosing the axis is part of choosing the question, because it decides which groups of values are combined.

Collapse the axis you name

For a (4, 3) table with four students and three tests, an aggregation without axis uses all 12 values. scores.sum() gives one total. scores.mean(axis=0) collapses the student axis and leaves three test means. scores.mean(axis=1) collapses the test axis and leaves four student means.

scores = np.array([[78, 85, 92], [64, 70, 58],
                   [88, 91, 95], [72, 68, 80]])

Its overall sum is 941. The column means are [75.5 78.5 81.25]; the row means are approximately [85. 64. 91.33 73.33]. Ask which axis should disappear. Saying “axis 0 means each row” reverses the operation: axis 0 combines rows.

Choose the quantity to return

Methods such as a.sum() and functions such as np.sum(a) express the same basic operation. NumPy's methods avoid the different behavior of Python's built-ins. For example, Python's sum(scores) adds row arrays to produce column totals, while scores.sum() combines all entries. Python's max(scores) tries to compare rows and raises an ambiguous-truth error.

Positions differ from values

scores.mean(axis=1).argmax() is 2: the third student has the highest mean. argmax returns the first position when maxima tie. Without an axis on a table, it returns a row-major flattened index; scores.argmax() is 8. np.unravel_index(8, scores.shape) recovers coordinate (2, 2).

Use keepdims=True when a summary must align with the original table. scores.sum(axis=1, keepdims=True) has shape (4, 1), making it suitable for one normalization value per row. Without it, the result has shape (4,).

Keep intermediate totals when order matters

np.cumsum([3, 1, 4, 1, 5]) gives [3 4 8 9 14]. Each position is the total through that point. np.cumprod([1, 2, 3, 4]) similarly gives [1 2 6 24]. np.diff([3, 1, 4, 1, 5]) gives the four neighbor changes [-2 3 -3 4]. On a table, specify axis=1 for separate running totals across rows.

State the spread convention

std and var divide the squared-deviation total by the number of values by default: ddof=0. Use ddof=1 for the sample formula, whose divisor is one smaller. For [2, 4, 4, 4, 5, 5, 7, 9], the mean is 5, the default standard deviation is 2.0, and the sample standard deviation is about 2.14. A report should make clear which convention it uses.

Order-based summaries describe where values sit in a sorted distribution. Sorting also supports ranking records, but a table needs careful treatment to keep each record's fields together.

Locate the middle and the quartiles

For x = np.array([12, 15, 11, 20, 18, 14]), the sorted order is [11 12 14 15 18 20]. With six values, np.median(x) averages the middle two, 14 and 15, giving 14.5. For an odd count, the median is the one middle value.

np.percentile(x, q) takes a percentage from 0 to 100. np.quantile(x, q) takes the equivalent fraction from 0 to 1. Their default linear interpolation can give a value between observations. For six values, the 25th percentile's zero-based position is 0.25 × 5 = 1.25. That lies one quarter of the way from 12 to 14, giving 12.5. The 50th and 75th percentiles are 14.5 and 17.25, so the interquartile range is 4.75.

The median resists an extreme value more than the mean does. Adding 200 to these data moves the mean from 15.0 to about 41.4; the median moves from 14.5 to 15.0. Neither measure alone describes the whole distribution.

Choose sorted data or sorting positions

np.sort(x) returns a sorted copy. x.sort() changes x and returns None, so y = x.sort() does not put sorted values into y. For descending values, use np.sort(x)[::-1].

On a 2-D array, sorting defaults to the last axis, sorting each row independently. axis=0 sorts each column independently. Neither operation is a record sort. For tab = np.array([[3, 1], [1, 2], [2, 0]]), np.sort(tab, axis=0) gives [[1 0] [2 1] [3 2]], separating values that originally belonged together.

To preserve records, get an ordering from the key column and apply it to whole rows:

order = tab[:, 0].argsort()
ranked = tab[order]
# array([[1, 2], [2, 0], [3, 1]])

The same index order can reorder related vectors. If avg = np.array([85.0, 64.0, 91.33, 73.33]) and names = np.array(['Ana', 'Ben', 'Chen', 'Dara']), order = np.argsort(avg)[::-1] ranks the means in descending order. names[order] gives ['Chen' 'Ana' 'Dara' 'Ben'].

Count distinct values and test membership

np.unique returns sorted distinct values. With return_counts=True, it also returns corresponding counts:

colors = np.array(['red', 'blue', 'red', 'green', 'blue', 'red'])
values, counts = np.unique(colors, return_counts=True)

The outputs are ['blue' 'green' 'red'] and [2 1 3]. values[counts.argmax()] finds the most frequent value, 'red'. The count positions match the unique-value positions.

np.isin(a, allowed) returns a mask for membership in a collection. For a = np.arange(1, 6), np.isin(a, [2, 4]) gives [False True False True False], ready to select those values or their associated records.

Missing data and filtered subsets need an explicit summary policy. A numerical result can be wrong even when the array operation runs successfully: a placeholder can be counted as a measurement, or an empty subset can be mistaken for a meaningful average.

Find missing values before summarizing

NaN propagates through ordinary sums and means. For x = np.array([3.0, np.nan, 5.0]), both x.sum() and x.mean() give NaN. np.nansum(x) gives 8.0 and np.nanmean(x) gives 4.0, skipping the missing position. Related functions include nanmax, nanmin, nanstd, and nanmedian. If every value is missing, nanmean returns NaN with a warning because no observed values remain.

For sensor readings in °C:

r = np.array([21.5, np.nan, 22.1, 23.0, np.nan, 21.9])
missing = np.isnan(r)

missing.sum() is 2, missing.mean() is about 0.33, and r[~missing] keeps four readings. np.nanmean(r) is 22.125, using only those four values. Comparing r == np.nan would fail to mark the missing readings.

Convert placeholder codes first

A number such as −999 may represent missing data in an input system. NumPy treats it as a real measurement until you replace it:

v = np.array([12, -999, 15, -999, 9]).astype(float)
v[v == -999] = np.nan
np.nanmean(v)  # 12.0

Convert to float first because an integer array cannot store NaN. This replacement is appropriate only when −999 is defined as the missing code for those data.

Filling changes the distribution

Replacing missing temperatures with zero makes the example's mean 14.75, even though no zero-temperature measurement was observed. Filling them with the observed mean, 22.125, preserves that mean but reduces the apparent spread. A filled value is an introduced estimate, so its effect belongs in the analysis decision.

For a 2-D table D, np.nanmean(D, axis=0) summarizes each column's observed values. Each column can have a different count. To drop rows containing any missing value:

complete = ~np.isnan(D).any(axis=1)
clean_rows = D[complete]

The reduction creates one decision per row; negation retains rows with no missing entries. With D = np.array([[1., 2.], [np.nan, 4.], [5., np.nan], [7., 8.]]), this keeps [[1. 2.] [7. 8.]].

Summarize a condition or a group

For temperatures t, t[t > 20].mean() averages only days above 20 °C. In grouped data, use a matching label mask: with group = np.array(['A', 'B', 'A', 'B', 'A']) and value = np.array([10, 20, 14, 26, 12]), value[group == 'A'].mean() is 12.0 and the corresponding B mean is 23.0. Iterate over np.unique(group) to repeat this for each group.

Before summarizing a filtered subset, check mask.any(). An empty selection's mean is NaN with a warning, while its sum is 0. That zero can be a valid empty total, but it does not establish that any measurement was present.

A random number generator supports sampling, shuffling, and simulation. A seed makes the starting state repeatable, so you can rerun an analysis while investigating its logic.

Create one Generator and reuse it

rng = np.random.default_rng(42)

NumPy recommends this Generator interface for new code. Its default bit generator is PCG64. Without a seed, default_rng() obtains fresh entropy from the operating system. Older calls such as np.random.seed and np.random.rand use shared global state; Generator objects let an analysis manage its own state directly.

Pseudorandom values come from a deterministic algorithm. The same starting state and same sequence of calls reproduce the stream in the same environment. Every draw advances the state. Creating a new Generator with seed 42 inside every loop pass restarts the stream, repeating the same initial draws instead of continuing it. Record the NumPy version and the code as well as the seed, because stream compatibility is not guaranteed across versions or all environments. For passwords or security tokens, use Python's secrets module.

Pick a distribution and an output shape

size can be an element count or a shape tuple. For normal, loc is the distribution mean and scale is its standard deviation. A variance of 100 therefore needs scale=10. A small sample's actual mean and standard deviation need not match the distribution parameters exactly.

For integer draws, the upper bound is excluded by default. A die needs integers(1, 7), since integers(1, 6) cannot draw 6. integers(1, 6, endpoint=True) includes the upper bound and also models a die.

Sample values with or without replacement

rng.choice(['heads', 'tails'], size=5) makes five draws with replacement. A value can appear several times. For distinct sampled positions, use replace=False; requesting more positions than the population has raises ValueError. Probabilities can be supplied as p=[0.7, 0.3], in the same order as the values, and must sum to 1.

Without replacement prevents selecting an input position twice. If the input itself contains duplicate values, equal values can still appear in the sample.

Shuffle records without separating their fields

rng.shuffle(a) rearranges the existing array and returns None. rng.permutation(a) returns a shuffled copy. rng.permutation(10) gives integers 0 through 9 in random order, useful as row indexes.

For a 2-D array, the default axis is 0: complete rows move together. That keeps each observation's fields aligned. To shuffle two related arrays consistently, generate one permutation of their shared record indexes and apply it to both.

A vectorized simulation represents repeated trials as arrays, applies the event rule to those arrays, and summarizes the results. The array shape should show what one trial contains.

Turn an event into a proportion

rng = np.random.default_rng(42)
rolls = rng.integers(1, 7, size=100_000)
estimate = (rolls == 6).mean()

The mask is true once for every observed 6. Its mean estimates the probability of that event, close to 1/6, or about 0.1667. rolls.mean() answers a different question: the average die face, close to 3.5. A fixed seed repeats an estimate; it does not make that estimate exact.

Vectorization avoids Python-level work for each draw, often making large simulations much faster than a loop. The speedup depends on the calculation, array size, and computer, so it is not a fixed property of the event.

Use one row per multi-value trial

To simulate two dice in each trial:

dice = rng.integers(1, 7, size=(100_000, 2))
totals = dice.sum(axis=1)
p_seven = (totals == 7).mean()

The table has 100,000 trial rows and two die columns. Reducing axis 1 produces one total per trial. Six of the 36 equally likely ordered pairs sum to 7, so the estimated probability should be near 6/36. Reducing axis 0 would instead give two totals across all trials and lose the trial-level event.

Larger simulations still vary

For independent trials, the standard error of an estimated proportion is:

Here $p$ is the event probability and $n$ the trial count. For a die's 6 with 100,000 trials, it is about 0.0012. A normal approximation puts about 95% of estimates within 0.0023 of 1/6. Increasing the trial count by a factor of 100 reduces standard error by a factor of 10. More draws do not correct an incorrectly coded event.

Build and inspect a random walk

steps = rng.choice([-1, 1], size=1000) draws left or right steps. walk = steps.cumsum() gives the position after each step. Its maximum and minimum locate the extreme positions visited. To find the first point at least 10 units from the start:

hit = np.abs(walk) >= 10
if hit.any():
    step_count = hit.argmax() + 1

argmax locates the first true entry, but also returns 0 when every entry is false. The any() check distinguishes a real hit from no hit. Add 1 because array position 0 is the position after step 1.

For many walks, rng.choice([-1, 1], size=(5000, 1000)).cumsum(axis=1) accumulates separately along each row. walks[:, -1] selects all final positions; their mean is near 0 and standard deviation near the square root of 1,000, about 31.6.

Use permutations and exact checks

A random split of 100 records can use idx = rng.permutation(100), then data[idx[:80]] and data[idx[80:]]. The two disjoint slices of one permutation assign every record exactly once. For a simulation, compare a simple case with an exact answer, such as a die mean of 3.5, before trusting more complicated output.