NumPy

Lesson 4 — Python

Lesson 4 of 20 Intermediate ~75 min

Learning objectives

Create and inspect NumPy arrays
Use vectorised operations instead of loops
Apply broadcasting rules deliberately
Index and slice arrays, including boolean masks
Handle missing values with nan
Understand where views differ from copies

Why NumPy

Python lists are flexible and slow. NumPy arrays are homogeneous, contiguous in memory, and operations on them run in compiled C.

import numpy as np

# Pure Python
values = list(range(1_000_000))
squares = [x**2 for x in values]        # ~60 ms

# NumPy
arr = np.arange(1_000_000)
squares = arr**2                        # ~1 ms

NumPy is also the foundation of pandas, scikit-learn, SciPy and every other scientific Python library. Understanding arrays is prerequisite to understanding pandas.

This is exactly R’s vectorisation argument. An R vector is essentially a NumPy array — contiguous, homogeneous, operated on element-wise in compiled code.

Creating arrays

np.array([1, 2, 3])                     # from a list
np.array([[1, 2], [3, 4]])              # 2-D
np.zeros(5)                             # [0. 0. 0. 0. 0.]
np.ones((2, 3))                         # 2x3 of ones
np.full(5, 7)                           # [7 7 7 7 7]
np.empty(5)                             # uninitialised — contains garbage
np.arange(0, 10, 2)                     # [0 2 4 6 8]
np.linspace(0, 1, 5)                    # [0. 0.25 0.5 0.75 1.]
np.eye(3)                               # identity matrix
np.random.default_rng(42).normal(0, 1, 100)   # modern RNG API

Inspect:

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

a.shape         # (2, 3)
a.ndim          # 2
a.size          # 6
a.dtype         # dtype('int64')
a.nbytes        # 48
a.T             # transpose

dtypes

np.array([1, 2, 3]).dtype               # int64
np.array([1.0, 2, 3]).dtype             # float64
np.array([1, 2, 3], dtype=np.float32)
np.array(["a", "b"]).dtype              # <U1  (unicode, length 1)
np.array([True, False]).dtype           # bool

a.astype(np.float64)

A NumPy array is homogeneous — all elements share one dtype. Mixing types promotes, exactly as R’s coercion rules do:

np.array([1, 2.5])          # array([1. , 2.5])       -> float64
np.array([1, "a"])          # array(['1', 'a'])       -> <U21

Integer overflow is silent

np.array([2**62], dtype=np.int64) * 4
# array([0])          — wrapped around, no warning

np.array([2**62], dtype=np.float64) * 4
# array([1.8446744e+19])

R’s integers are 32-bit and produce NA with a warning on overflow. NumPy’s default is 64-bit and wraps silently. For counts and IDs this rarely matters; for accumulating large sums, use float64 or np.int64 with awareness.

Vectorised operations

a = np.array([1, 2, 3, 4])
b = np.array([10, 20, 30, 40])

a + b               # [11 22 33 44]
a * b               # [10 40 90 160]
a ** 2              # [1 4 9 16]
np.sqrt(a)
np.log(a)
np.exp(a)
np.abs(-a)
np.round(a / 3, 2)

a > 2               # [False False True True]
(a > 2) & (b < 35)  # element-wise AND — note the parentheses
~(a > 2)            # NOT

Use &, |, ~ — not and, or, not

(a > 2) and (b < 35)
# ValueError: The truth value of an array with more than one element is ambiguous

(a > 2) & (b < 35)          # correct

and/or try to evaluate the whole array as a single boolean. NumPy and pandas require the bitwise operators, and the parentheses are mandatory because & binds more tightly than >.

This is identical to R’s & versus && distinction.

Aggregations:

a.sum(); a.mean(); a.std(); a.min(); a.max()
a.cumsum(); a.cumprod()
np.median(a); np.percentile(a, [25, 50, 75])
np.corrcoef(a, b)

m = np.array([[1, 2, 3], [4, 5, 6]])
m.sum()             # 21     everything
m.sum(axis=0)       # [5 7 9]     down the columns
m.sum(axis=1)       # [6 15]      across the rows

axis=0 collapses rows (giving one value per column); axis=1 collapses columns. Remembering which is which: the axis you name is the one that disappears.

Broadcasting

Operations between arrays of different shapes, where NumPy stretches the smaller one.

a = np.array([1, 2, 3])
a * 2                       # [2 4 6]           scalar broadcast

m = np.array([[1, 2, 3],
              [4, 5, 6]])
m + np.array([10, 20, 30])  # adds the row vector to EVERY row
#> [[11 22 33]
#>  [14 25 36]]

m + np.array([[10], [20]])  # adds the column vector to every column
#> [[11 12 13]
#>  [24 25 26]]

The rules, applied from the trailing dimension backwards:

Dimensions are compatible if they are equal, or one of them is 1
Missing dimensions are treated as 1
Otherwise it is an error

(3, 4) + (4,)      -> (3, 4)    OK
(3, 4) + (3, 1)    -> (3, 4)    OK
(3, 4) + (3,)      -> ERROR     4 vs 3

Centring a matrix by column:

data = np.random.default_rng(0).normal(50, 10, (100, 5))
centred = data - data.mean(axis=0)         # (100,5) - (5,) -> broadcast
standardised = centred / data.std(axis=0)

Broadcasting can silently do the wrong thing

Like R’s recycling, broadcasting is powerful and occasionally produces a plausible result from a mistake. If two arrays happen to have compatible shapes for the wrong reason, you get numbers rather than an error.

Assert shapes when it matters:

assert data.shape[1] == len(weights), (
    f"weights has {len(weights)} elements, data has {data.shape[1]} columns"
)

Indexing

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

a[0]                # 10
a[-1]               # 50
a[1:4]              # [20 30 40]
a[::2]              # [10 30 50]
a[::-1]             # reversed

m = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
m[0, 1]             # 2          row 0, column 1
m[0]                # [1 2 3]    row 0
m[:, 1]             # [2 5 8]    column 1
m[0:2, 1:3]         # [[2 3], [5 6]]
m[[0, 2], :]        # rows 0 and 2 — fancy indexing

Boolean masks

The workhorse of data selection:

ages = np.array([45, 72, 38, 66, 29])

mask = ages >= 65
mask                # [False True False True False]
ages[mask]          # [72 66]
ages[ages >= 65]    # same, inline

ages[(ages >= 40) & (ages < 70)]        # [45 66]

np.where(ages >= 65, "elderly", "adult")
#> array(['adult', 'elderly', 'adult', 'elderly', 'adult'], dtype='<U7')

np.where(ages >= 65)                     # indices where True
#> (array([1, 3]),)

np.select(
    [ages < 18, ages < 65],
    ["child", "adult"],
    default="elderly",
)

np.select is the NumPy equivalent of dplyr::case_when().

Modification through a mask:

ages[ages > 100] = np.nan       # ValueError on an int array — nan is a float
ages = ages.astype(float)
ages[ages > 100] = np.nan       # now fine

Missing values

NumPy has no NA. Missingness is np.nan, which is a float.

a = np.array([1, 2, np.nan, 4])
a.dtype             # float64 — the nan forced it

a.sum()             # nan     — propagates
np.nansum(a)        # 7.0
np.nanmean(a)       # 2.333
np.nanstd(a)
np.nanmax(a)

np.isnan(a)         # [False False True False]
a[~np.isnan(a)]     # [1. 2. 4.]

np.nan == np.nan    # False!  — always use np.isnan()

nan is not equal to itself

np.nan == np.nan            # False
np.nan is np.nan            # True (same object, but do not rely on this)
np.isnan(np.nan)            # True — the correct test

This is IEEE 754 behaviour, shared by R’s NaN. The difference is that R has a separate NA for statistical missingness, typed per column; NumPy has only the float nan.

The practical consequence: an integer NumPy array cannot hold a missing value. pandas addressed this with nullable dtypes — see lesson 5.

Views versus copies

Slicing returns a view — a window onto the same memory.

a = np.array([1, 2, 3, 4, 5])
b = a[1:4]          # a VIEW
b[0] = 99
a                   # [1 99 3 4 5]   — a changed!

c = a[1:4].copy()   # an explicit copy
c[0] = 0
a                   # unchanged

b.base is a         # True  — b is a view of a
c.base is None      # True  — c owns its data

Fancy indexing and boolean masking return copies:

d = a[[1, 2, 3]]    # fancy indexing -> copy
e = a[a > 2]        # boolean mask   -> copy

Views are a performance feature — slicing a 10 GB array costs nothing. They are also a source of surprising mutation, in exactly the way R’s copy-on-modify semantics protect against. When in doubt, .copy().

Random numbers

rng = np.random.default_rng(42)         # the modern API

rng.normal(50, 10, 100)
rng.uniform(0, 1, 10)
rng.integers(1, 7, 10)                  # dice
rng.choice(["A", "B", "C"], 20)
rng.choice(subjects, 50, replace=False) # sample without replacement
rng.permutation(arr)
rng.shuffle(arr)                        # in place

Always seed for reproducibility. Use default_rng rather than the legacy np.random.seed() / np.random.normal() functions — the generator object is explicit, and independent generators do not interfere.

Linear algebra

A = np.array([[1, 2], [3, 4]])
b = np.array([5, 6])

A @ A                       # matrix multiplication
A * A                       # ELEMENT-WISE — not matrix multiplication
np.linalg.inv(A)
np.linalg.solve(A, b)       # solve Ax = b — better than inv(A) @ b
np.linalg.det(A)
np.linalg.eig(A)
np.linalg.cholesky(A)

@ versus * is the equivalent of R’s %*% versus *, and the mistake goes the same way.

R and Python side by side

Task	R	NumPy
Create	`c(1, 2, 3)`	`np.array([1, 2, 3])`
Sequence	`1:10`	`np.arange(1, 11)`
Repeat	`rep(0, 5)`	`np.zeros(5)`
Even spacing	`seq(0, 1, length.out = 5)`	`np.linspace(0, 1, 5)`
Length	`length(x)`	`x.size`
Dimensions	`dim(x)`	`x.shape`
Type	`typeof(x)`	`x.dtype`
Index	`x[1]` (1-based)	`x[0]` (0-based)
Slice	`x[2:4]` (inclusive)	`x[1:4]` (exclusive stop)
Boolean filter	`x[x > 5]`	`x[x > 5]`
Conditional	`ifelse(c, a, b)`	`np.where(c, a, b)`
Multi-branch	`case_when()`	`np.select()`
Missing	`NA`	`np.nan` (float only)
Ignore missing	`mean(x, na.rm = TRUE)`	`np.nanmean(x)`
Test missing	`is.na(x)`	`np.isnan(x)`
Matrix multiply	`A %*% B`	`A @ B`
Column sums	`colSums(m)`	`m.sum(axis=0)`
Assignment	Copies	View for slices

Common mistakes

Mistake	Consequence	Fix
`and`/`or` on arrays	`ValueError`	`&`, `\|`, with parentheses
Missing parentheses around comparisons	Wrong precedence	`(a > 2) & (b < 5)`
`x == np.nan`	Always `False`	`np.isnan(x)`
Expecting a slice to copy	Original mutated	`.copy()`
`*` for matrix multiplication	Element-wise result	`@`
Wrong `axis`	Aggregated the wrong way	The named axis disappears
Integer array assigned `nan`	`ValueError`	Cast to float first
Looping over an array	Slow	Vectorise

Exercise 4.1 — Vectorise a loop

Rewrite this without a loop, and measure the difference.

def normalise(values, lower, upper):
    result = []
    for v in values:
        if v is None:
            result.append(None)
        elif v < lower:
            result.append(0.0)
        elif v > upper:
            result.append(1.0)
        else:
            result.append((v - lower) / (upper - lower))
    return result

Show solution

import numpy as np


def normalise(values, lower, upper):
    """Scale values to [0, 1], clipping outside the range.

    Args:
        values: Array-like of numbers; NaN propagates.
        lower: Value mapped to 0.0.
        upper: Value mapped to 1.0. Must be greater than `lower`.
    """
    if upper <= lower:
        raise ValueError(f"upper ({upper}) must exceed lower ({lower})")

    arr = np.asarray(values, dtype=float)      # None becomes nan
    scaled = (arr - lower) / (upper - lower)
    return np.clip(scaled, 0.0, 1.0)           # nan survives clip

values = [None, 5.0, 15.0, 25.0, 35.0]
normalise(values, 10, 30)
#> array([ nan, 0. , 0.25, 0.75, 1.  ])

The three pieces:

np.asarray(values, dtype=float) converts a list to an array and turns None into nan in one step. asarray rather than array avoids a copy when the input is already an array.
The arithmetic is vectorised — one expression for the whole array.
np.clip replaces the two conditional branches, and propagates nan rather than clipping it.

Timing:

import timeit

data = list(np.random.default_rng(0).normal(20, 10, 1_000_000))

# Loop version
timeit.timeit(lambda: normalise_loop(data, 10, 30), number=1)
#> 0.42 s

# Vectorised
timeit.timeit(lambda: normalise(data, 10, 30), number=1)
#> 0.006 s

# Already an array — no conversion cost
arr = np.asarray(data)
timeit.timeit(lambda: normalise(arr, 10, 30), number=1)
#> 0.003 s

Roughly 70× faster on a list, 140× on an array. Most of the remaining cost in the list case is the list-to-array conversion, which is why keeping data in arrays (or a DataFrame) throughout a pipeline matters more than optimising any single operation.

A subtlety worth checking. Does np.clip handle nan correctly?

np.clip(np.array([np.nan, 0.5, 2.0]), 0, 1)
#> array([nan, 0.5, 1. ])

Yes — nan passes through. That is the behaviour we want, and it is worth verifying rather than assuming, because np.minimum/np.maximum and np.fmin/np.fmax differ from each other precisely on nan handling.

Exercise 4.2 — Reference range flags

Given arrays of lab values, lower limits and upper limits, compute a normal-range indicator ("LOW", "NORMAL", "HIGH", "MISSING") and the count in each category — without any loops. Handle missing limits.

Show solution

import numpy as np


def anrind(values, lower, upper):
    """Derive a normal range indicator.

    Args:
        values: Lab results. NaN gives "MISSING".
        lower: Lower limits of normal. NaN means no lower limit.
        upper: Upper limits of normal. NaN means no upper limit.

    Returns:
        An array of "LOW", "NORMAL", "HIGH" or "MISSING".
    """
    values = np.asarray(values, dtype=float)
    lower  = np.asarray(lower,  dtype=float)
    upper  = np.asarray(upper,  dtype=float)

    if not (values.shape == lower.shape == upper.shape):
        raise ValueError(
            f"Shape mismatch: values {values.shape}, "
            f"lower {lower.shape}, upper {upper.shape}"
        )

    # Comparisons with nan are False, which is what we want: a missing limit
    # means the value cannot be low/high on that side.
    with np.errstate(invalid="ignore"):
        is_low  = values < lower
        is_high = values > upper

    return np.select(
        [np.isnan(values), is_low, is_high],
        ["MISSING",        "LOW",  "HIGH"],
        default="NORMAL",
    )

values = np.array([15.0, 45.0, 80.0, np.nan, 30.0])
lower  = np.array([20.0, 20.0, 20.0, 20.0,   np.nan])
upper  = np.array([60.0, 60.0, 60.0, 60.0,   60.0])

result = anrind(values, lower, upper)
result
#> array(['LOW', 'NORMAL', 'HIGH', 'MISSING', 'NORMAL'], dtype='<U7')

labels, counts = np.unique(result, return_counts=True)
dict(zip(labels, counts))
#> {'HIGH': 1, 'LOW': 1, 'MISSING': 1, 'NORMAL': 2}

Three details that make this correct rather than nearly correct:

1. np.isnan(values) must be the first condition in np.select. np.select takes the first matching condition. A nan value compared with a limit gives False for both is_low and is_high, so without the explicit missing check it would fall through to the default and be labelled "NORMAL" — which is badly wrong for a lab dataset.

2. Missing limits are handled by the comparison semantics. The last row has no lower limit (nan), so 30.0 < nan is False and the value cannot be flagged "LOW". It is correctly "NORMAL" because it is below the upper limit. That falls out of IEEE 754 rather than needing special handling — but it is worth stating explicitly in a comment, because a reader will otherwise wonder whether it was considered.

3. np.errstate(invalid="ignore") suppresses the RuntimeWarning NumPy emits when comparing with nan. Without it, a 200,000-row lab dataset produces a wall of warnings. Scoping the suppression to just these two lines is important — a blanket warnings.filterwarnings("ignore") would hide genuine problems elsewhere.

Verification

# A value with NO limits at all must be NORMAL, not MISSING
anrind([30.0], [np.nan], [np.nan])
#> array(['NORMAL'], dtype='<U7')

# A value exactly on a limit is NORMAL (limits are inclusive)
anrind([20.0, 60.0], [20.0, 20.0], [60.0, 60.0])
#> array(['NORMAL', 'NORMAL'], dtype='<U7')

That second check encodes a convention: are the reference range limits inclusive? CDISC ADaM treats ANRLO <= AVAL <= ANRHI as normal, so strict < and > for the abnormal flags is correct. It is exactly the kind of boundary decision that double programming exists to catch, and it belongs in a test.

Recap

Arrays are homogeneous and contiguous; operations run in compiled C
&, |, ~ with parentheses — never and/or on arrays
Broadcasting stretches smaller arrays; assert shapes when correctness depends on it
axis=0 collapses rows, axis=1 collapses columns — the named axis disappears
nan is a float and is never equal to itself; use np.isnan() and the nan* functions
Slices are views; fancy indexing and masks are copies
np.where for two branches, np.select for many
np.random.default_rng(seed), not the legacy API

Next: pandas.