Mathematical and statistical judgment with NumPy

The linear-algebra and random-number lessons provide the starting point. Here we connect their calculations to three questions: what does a summary represent, how can we check a numerical solution, and what does a random experiment establish? This is preparation for later statistical modelling, not a claim that one short lesson replaces a full mathematics or statistics text.

Read a summary as a claim about observations

Suppose four invented reading times are 10, 12, 14 and 40 minutes. The mean uses every value; the median describes the middle of the ordered values. Neither tells the entire story. Before removing the longest reading time, investigate whether it is an error, a different kind of reading session, or a genuine observation relevant to your question.

import numpy as np
times = np.array([10.0, 12.0, 14.0, 40.0])
print(times.mean(), np.median(times))
print(times.std(ddof=0), times.std(ddof=1))
print(np.quantile(times, [0.25, 0.75]))
assert times.mean() == 19.0 and np.median(times) == 13.0

19.0 13.0
12.206555615733702 14.094916341243982
[11.5 20.5]

For observations $x_1,\ldots,x_n$, the mean is $\bar{x}=\sum_i x_i/n$. The descriptive population variance divides the squared deviations by $n$; the usual sample variance divides by $n-1$. In NumPy, ddof=1 selects the latter denominator. The standard-deviation reference explains this parameter. State which quantity you report and why. The standard deviation measures spread among observations; it is not automatically the uncertainty of their mean.

If observations are a time sequence, repeated measurements of the same person, or clustered within institutions, treating them as independent observations requires justification. Ten visits from one person do not provide the same information about variation between people as ten independently selected people.

Separate association from explanation

A positive association means larger values of one variable tend to occur with larger values of another in the observed data. It does not identify what would happen after intervening. Group composition, time, selection and other variables can change the apparent relationship. The NIST correlation discussion provides further background. A correlation also describes a particular kind of relationship; a small linear correlation need not imply the absence of structure.

hours = np.array([1., 2., 3., 4.])
items = np.array([2., 4., 6., 8.])
print(np.corrcoef(hours, items)[0, 1])
assert np.isclose(np.corrcoef(hours, items)[0, 1], 1.0)

1.0

This constructed example has perfect linear association because we defined items to be twice hours. It is an arithmetic demonstration, not empirical evidence about real reading behaviour. In an actual investigation, inspect the units, observation process, sample size, subgroups and plausible alternative explanations before writing a conclusion.

Solve a system and check the residual

The matrix chapter introduces determinants and inverses. An inverse is a mathematical object worth understanding, but you do not need to form it explicitly to solve $A\mathbf{x}=\mathbf{b}$. A residual, $A\hat{\mathbf{x}}-\mathbf{b}$, checks how closely the computed answer satisfies the supplied equations.

A = np.array([[2., 1.], [1., 3.]])
b = np.array([8., 9.])
solution = np.linalg.solve(A, b)
print(solution)
print(A @ solution - b)
assert np.allclose(solution, [3., 2.])
assert np.allclose(A @ solution, b)

[3. 2.]
[0. 0.]

For an overdetermined system, a least-squares calculation can find coefficients that minimize the sum of squared residuals. The NumPy least-squares reference defines the returned coefficients, residual information, rank and singular values. Use a solver rather than presenting explicit inversion of the normal equations as the general computational recipe.

X = np.column_stack([np.ones(4), hours])
coefficients, residuals, rank, singular_values = np.linalg.lstsq(X, items, rcond=None)
print(coefficients, rank)
assert np.allclose(coefficients, [0., 2.])

[4.85700758e-16 2.00000000e+00] 2

A small residual does not prove that a real-world model is suitable. It checks the fit to these values under this representation. Nearly dependent columns can make coefficient interpretation unstable. Rank asks how many independent directions the columns contain. Eigenvectors and eigenvalues, introduced in later mathematical work, describe invariant directions of square transformations; they prepare you for PCA, but estimating PCA components and choosing their number belong in the modelling sequence.

Use a local random generator

Use np.random.default_rng(seed) for a local generator rather than depending on unrelated cells changing shared random state. The Generator documentation describes this interface. Save package versions as well as the seed; do not promise that every future implementation will produce an identical sequence.

rng = np.random.default_rng(23)
repeated = rng.choice(times, size=(2000, len(times)), replace=True).mean(axis=1)
print(np.quantile(repeated, [0.025, 0.975]))
assert repeated.min() >= times.min() and repeated.max() <= times.max()

[11. 33.]

Each row resamples the four observations with replacement, then calculates a mean. The quantiles describe this resampling distribution. Interpreting them as an uncertainty interval needs assumptions about the sampling process; four convenient observations do not become representative because you resample them thousands of times. Dependence may require resampling groups or blocks, and a biased measurement remains biased. See SciPy’s bootstrap reference for a fuller method and its limitations. Do not interpret the two quantiles as a guarantee that a specific unobserved population mean lies between them.

Practice and independent check

Replace 40 with 16 and predict which summaries change most. Solve a different two-equation system by hand and verify it with a residual. Create two identical columns in a design matrix and inspect its rank; explain why two separately interpretable coefficients cannot be identified from those columns alone. Finally, write down the observation-selection assumptions required before using a resampling interval in your own project.

Transition

Numerical foundations and interpretation belong together. Keep the original NumPy cultural-production challenge as a checkpoint, then move to labelled tables, where alignment and observation identity become just as important as arithmetic.

Chapter exercisesWorked answers after an attempt