If NumPy is the foundation of scientific Python and matplotlib draws the pictures, SciPy is the toolbox in between. It’s been there since 2001, it’s enormous, and most working Python developers only touch one or two of its submodules before forgetting it exists. That’s a shame, because the moment you need a t-test, a curve fit, a sparse matrix, or an FFT, the answer is almost always already in scipy.something and the docs are excellent.
This lesson is a tour. Not a reference — that would be a book — but enough of a map that you know which submodule to open when you hit a problem that NumPy alone doesn’t cover. We’ll go deep on the three submodules that come up most often (stats, optimize, sparse) and quickly name the rest.
How SciPy is organized
SciPy is a collection of submodules, each focused on a domain. You almost never import scipy and use it directly; you import the submodule you need:
from scipy import stats
from scipy import optimize
from scipy import sparse
from scipy import signal
from scipy import interpolate
from scipy import spatialThe library is on version 1.x in 2026 — it’s been on the 1.x line since 2017 and the API is rock-solid. Things rarely break between versions. Install with uv add scipy; on most systems it pulls in NumPy automatically.
scipy.stats: probability and statistical tests
This is the submodule I open most often. It has two main pieces: probability distribution objects, and statistical tests.
Distributions are objects that bundle the PDF, CDF, sampling, and parameter estimation for a named distribution:
from scipy import stats
import numpy as np
# Standard normal
stats.norm.pdf(0) # 0.3989... (peak of the bell curve)
stats.norm.cdf(1.96) # 0.975 — the famous 95% boundary
stats.norm.ppf(0.975) # 1.959... — inverse CDF, gives you the quantile
stats.norm.rvs(size=1000) # 1000 random samples
# Non-standard? Pass loc and scale
stats.norm(loc=100, scale=15).rvs(size=1000) # IQ-like distribution
# Other distributions follow the same interface
stats.binom(n=10, p=0.3).pmf(3) # P(X=3) for Binomial(10, 0.3)
stats.poisson(mu=4.5).cdf(6) # P(X<=6) for Poisson(4.5)
stats.t(df=20).ppf(0.975) # critical value for a t with 20 dfThe fact that every distribution has the same .pdf / .cdf / .ppf / .rvs / .fit interface is the under-celebrated nice thing about this module. Once you know the pattern you can use any of the dozens of distributions without learning a new API.
Statistical tests is where the practical wins are. The ones that earn their keep in real data work:
# Two-sample t-test — A/B test analysis, classic
control = np.array([2.1, 2.3, 1.9, 2.5, 2.0, 2.2])
treatment = np.array([2.4, 2.6, 2.5, 2.7, 2.3, 2.5])
result = stats.ttest_ind(control, treatment)
print(result.statistic, result.pvalue)
# Mann-Whitney U — non-parametric alternative when data isn't normal
stats.mannwhitneyu(control, treatment)
# Chi-squared test of independence — for contingency tables
table = np.array([[50, 30], [20, 40]])
chi2, p, dof, expected = stats.chi2_contingency(table)
# Pearson and Spearman correlation
stats.pearsonr(control, treatment)
stats.spearmanr(control, treatment)When this matters in day-to-day data work: A/B test analysis (t-test or Mann-Whitney depending on your data shape), distribution fitting to figure out what kind of variable you’re dealing with, quantile-based outlier detection (stats.norm.ppf(0.99) to find the 1% upper tail of a fitted distribution), and the occasional chi-squared when you’re comparing categorical groups.
There’s also stats.describe(arr) which gives you mean, variance, skewness, kurtosis, and the min/max in one call — handy for a quick look at a column’s shape.
scipy.optimize: minimization, root-finding, curve fitting
Three jobs: find the minimum of a function, find where a function equals zero, and fit parameters of a model to data.
Minimization is the most general:
from scipy import optimize
# Minimize a simple parabola
def f(x):
return (x - 3) ** 2 + 1
result = optimize.minimize(f, x0=0.0)
print(result.x) # [3.] — the minimum
print(result.fun) # 1.0 — the value at the minimum
print(result.success) # Truex0 is the starting guess. For convex functions the starting point barely matters; for non-convex ones it matters a lot. optimize.minimize accepts an N-dimensional function — x0 can be an array of any size — and picks a sensible default solver (BFGS for unconstrained, L-BFGS-B if you pass bounds, SLSQP for constraints).
Curve fitting is the one most data people actually use, even if they’ve forgotten it lives here:
# Fit a logistic curve to growth data
def logistic(t, L, k, t0):
return L / (1 + np.exp(-k * (t - t0)))
t = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
y = np.array([0.1, 0.15, 0.3, 0.6, 1.1, 1.8, 2.5, 2.9, 3.1, 3.15, 3.18])
popt, pcov = optimize.curve_fit(logistic, t, y, p0=[3.0, 1.0, 4.0])
L_fit, k_fit, t0_fit = popt
print(f"L={L_fit:.2f}, k={k_fit:.2f}, t0={t0_fit:.2f}")curve_fit takes your model function (first arg is the independent variable, the rest are parameters), the x and y data, and an optional p0 initial guess. It returns the best-fit parameters and the covariance matrix (the diagonal of which gives you parameter standard errors). This is how you fit logistic curves to product growth, exponentials to decay, Gompertz to disease spread, anything where you have a known functional form and want the coefficients.
Root finding for when you want to solve f(x) = 0:
optimize.brentq(lambda x: x**3 - 2*x - 5, 1, 3) # root inside [1, 3]brentq is robust and almost always the right choice when you have a bracketing interval where the function changes sign.
scipy.sparse: matrices that are mostly zero
The third submodule that earns its keep. A sparse matrix stores only the non-zero entries. For a million-by-million matrix where 99.99% of entries are zero — extremely common in NLP, recommendation systems, graph problems — sparse storage is the difference between “fits in 2 GB of RAM” and “would need 8 TB.”
from scipy import sparse
# Build a sparse matrix from coordinate triplets
rows = np.array([0, 1, 2, 0])
cols = np.array([0, 1, 2, 2])
vals = np.array([1.0, 2.0, 3.0, 4.0])
M = sparse.csr_matrix((vals, (rows, cols)), shape=(3, 3))
print(M.toarray())
# [[1. 0. 4.]
# [0. 2. 0.]
# [0. 0. 3.]]The four formats you’ll see: csr_matrix (compressed sparse row — fast row slicing, fast matrix-vector products), csc_matrix (column version of the same thing), coo_matrix (the triplet format, good for building a sparse matrix incrementally then converting), and lil_matrix (the one to use when you need to mutate entries). The pattern: build in COO or LIL, convert to CSR for computation.
The classic use: term-document matrices. A corpus with 100k documents and a vocabulary of 50k words is, in dense form, a 5-billion-cell matrix. As a sparse csr_matrix it’s a few hundred megabytes. scikit-learn’s TfidfVectorizer returns a csr_matrix directly; you can pass it straight to a LogisticRegression and the whole pipeline stays sparse end to end.
The rest of the toolbox
The submodules I haven’t gone into, but you should know exist by name:
scipy.interpolate— when you have data points and want a continuous estimate between them.interp1dfor 1-D,RegularGridInterpolatorfor n-D on a grid, splines viaUnivariateSpline. The “I have lookup-table data and need to evaluate it at arbitrary points” submodule.scipy.signal— signal processing. Designing and applying digital filters (butter,filtfilt), computing FFTs (the olderscipy.fftis the modern home for this), correlations, convolutions, peak finding (find_peaks). If you’re working with sensor data, audio, or any time series where frequency content matters, this is your submodule.scipy.spatial— distance metrics and spatial data structures.cKDTreefor fast nearest-neighbor queries (much faster than brute force on more than a few thousand points),distance.cdistfor all-pairs distances between two sets of points.scipy.cluster— hierarchical clustering (linkage,dendrogram,fcluster) and a basic k-means. For most clustering work in 2026 people reach for scikit-learn instead, but the hierarchical functions in scipy are still the standard for that specific algorithm.scipy.integrate— numerical integration (quadfor 1-D integrals) and ODE solvers (solve_ivp).scipy.linalg— linear algebra. Mostly a superset ofnumpy.linalgwith extras: matrix factorizations, special solvers, thelstsqeveryone uses for least-squares.scipy.ndimage— n-dimensional image processing. Filters, morphology, labeling. The thing scikit-image is built on.
The pattern with all of these: you don’t memorize them. You go through this lesson once, you remember “right, that’s a SciPy thing,” and when you hit the problem in two months you open scipy.org/doc/scipy/reference/ and find the function in five minutes.
A real example: A/B test, end to end
Putting stats to work on something concrete. You ran an A/B test on a checkout button. Group A saw the old design, group B saw the new one. You logged conversion times in seconds for users who completed checkout. Did the new design make people faster?
import numpy as np
from scipy import stats
rng = np.random.default_rng(42)
group_a = rng.lognormal(mean=2.0, sigma=0.5, size=500)
group_b = rng.lognormal(mean=1.85, sigma=0.5, size=500)
# Quick descriptives
print(f"A: median={np.median(group_a):.2f}, mean={group_a.mean():.2f}")
print(f"B: median={np.median(group_b):.2f}, mean={group_b.mean():.2f}")
# T-test on raw values — but wait, lognormal data is skewed
t_result = stats.ttest_ind(group_a, group_b)
print(f"t-test: p={t_result.pvalue:.4f}")
# Better: log-transform first, since underlying data is lognormal
log_a, log_b = np.log(group_a), np.log(group_b)
t_log = stats.ttest_ind(log_a, log_b)
print(f"log t-test: p={t_log.pvalue:.4f}")
# Or skip the parametric assumption entirely
mw = stats.mannwhitneyu(group_a, group_b, alternative="greater")
print(f"Mann-Whitney: p={mw.pvalue:.4f}")Three answers to roughly the same question, each with different assumptions. This is what scipy.stats is for — having the tools to ask “is this real?” in two or three ways and seeing whether they agree.
A second example: interpolating sensor readings
A second pattern I hit constantly. You have a sensor that logs at irregular intervals. You want a clean reading every minute for downstream processing. scipy.interpolate.interp1d does the work:
from scipy import interpolate
# Irregular timestamps (in seconds since start) and the sensor values at each
t_irregular = np.array([0, 7, 23, 41, 58, 79, 102, 130, 165])
values = np.array([20.1, 20.3, 21.0, 22.4, 23.1, 22.8, 21.9, 20.5, 19.8])
# Build the interpolator (linear by default; "cubic" for smooth)
f = interpolate.interp1d(t_irregular, values, kind="cubic")
# Sample on a clean grid
t_regular = np.arange(0, 165, 5)
values_regular = f(t_regular)f is now a callable — give it any time inside the original range and it returns an interpolated value. kind="linear" is the default and the right choice when you have noisy data; kind="cubic" is smoother but can ring around outliers. For data outside the original range, you’ll want fill_value="extrapolate" or bounds_error=False — just be aware that extrapolation past your training data is a confidence trap.
For 2-D and higher, RegularGridInterpolator is the modern API that replaces the older, deprecated griddata for grid-structured data.
Wrapping the module
That’s Module 8 done. Lesson 43 gave you NumPy’s array model and broadcasting. Lesson 44 covered the three plotting libraries. Lesson 45 — this one — pointed at the rest of the scientific toolbox.
The pattern across all three lessons is the same: a small core of operations covers the vast majority of real work, and the libraries are big enough that nobody memorizes everything. You learn the shape of the toolkit once, you bookmark the docs, and you go look things up when you need them. That habit — “I know SciPy has a thing for this, let me find it” — is more valuable than memorizing function signatures, which the IDE and the LLM can both supply on demand anyway.
Module 9 picks up the thread and starts using all of this for real ML work: scikit-learn for classical models, then a pivot toward the modern stack (PyTorch basics, hugging face for NLP, the bits of MLOps that aren’t terrible). The numerical foundation we just built is what makes any of that possible.
Reference: SciPy documentation, retrieved 2026-05-01.