SciPy: Turning Measurements into Answers

You have some numbers — temperatures you wrote down, sales against advertising spend, a signal from a sensor — and a question about them. SciPy is the Python library that answers that kind of question: fit a curve to what you measured, find where something crosses a value, work out the best setting, add up an area, fill in the gaps between readings, or check whether a difference is real.

This is a lesson, not a tour. SciPy is enormous, and a list of its functions teaches nothing. So each section here starts with a problem, solves it, and shows exactly what comes back.

Where to run these

Put each example in a CodeBook cell or a Python code document. CodeBook is the better fit for this lesson, because it keeps your variables between Python cells — so the curve you fit in the first example is still there in the second, and you can build the work up a piece at a time exactly as it is written here.

Python runs on the device you are holding. One thing to know before you start: the first time you run anything that imports SciPy on a device, it downloads about eleven megabytes, so that first run needs a connection and takes a moment. Every run after that on the same device starts straight away. Nothing to set up, and nothing to clean up.

The version on your device is SciPy 1.11.2, alongside NumPy 1.25.2 on Python 3.11. Every example and every output on this page was run on exactly that.

If you have set a Python step to run on a connected computer instead of on the device in front of you, it uses the Python installed there — so SciPy needs to be installed on that machine too, and its version may differ. See Languages, and Where They Run.

Fitting a curve to measurements

This is the one to learn first, because it is the most common real problem and the others build on it.

Someone poured a cup of coffee and wrote down its temperature for an hour. The readings are a bit noisy, as real readings are. The questions are: how fast is it cooling, and what is the room temperature it is heading towards?

You know the shape of the answer — things cool towards room temperature, quickly at first and slower as they get closer. What you don't know are the three numbers that pin that shape down. curve_fit finds them.

import numpy as np
from scipy.optimize import curve_fit

minutes = np.array([0, 2, 5, 10, 15, 20, 30, 45, 60])
celsius = np.array([93.1, 86.6, 78.1, 67.0, 57.7, 50.4, 39.7, 30.3, 25.6])

def cooling(t, room, start, k):
    return room + start * np.exp(-k * t)

params, covariance = curve_fit(cooling, minutes, celsius, p0=[20, 70, 0.05])
room, start, k = params
spread = np.sqrt(np.diag(covariance))

print(f"room temperature : {room:.1f} C  (give or take {spread[0]:.1f})")
print(f"cooling rate     : {k:.4f} per minute")
print(f"time to halve    : {np.log(2) / k:.1f} minutes")

room temperature : 20.6 C  (give or take 0.3)
cooling rate     : 0.0444 per minute
time to halve    : 15.6 minutes

Nobody measured the room. The fit worked it out from the shape of the cooling — 20.6 °C, which is a believable room.

Three things in that example are worth taking with you:

• You write the shape as an ordinary Python function. cooling takes the thing you varied first (t), then one argument per unknown. That is the whole interface. Any shape you can write, you can fit.
• p0 is your rough guess, and it matters more than beginners expect. Fitting is a search, and it starts from p0. A sensible guess — "room is about 20, it started about 70 above that, it cools a few percent a minute" — finds the answer. A wild guess can wander off and fail. When a fit refuses to converge, a better p0 is the first thing to try.
• The second return value is how much to trust the first. covariance is usually ignored and shouldn't be: the square root of its diagonal gives the uncertainty on each fitted number. "20.6, give or take 0.3" is a result. "20.6" on its own is a number with unknown standing.

Asking the fitted curve a question

Now that cooling has real numbers in it, you have a model of the coffee — not just the nine moments you measured, but every moment in between. So you can ask it things nobody wrote down.

When was it cool enough to drink — say, 60 °C?

That is a root-finding problem: find the time where the temperature minus 60 is zero. brentq does it, given a range you know the answer sits inside.

from scipy.optimize import brentq

drinkable = brentq(lambda t: cooling(t, room, start, k) - 60.0, 0, 60)
print(f"cool enough to drink after {drinkable:.1f} minutes")

cool enough to drink after 13.6 minutes

In a CodeBook this runs as a second cell, with no imports repeated and no data re-entered — cooling, room, start and k are all still there from the first cell.

The 0, 60 is the bracket: the two ends of a range where the answer changes sign. brentq needs that and rewards you for it by never failing to find the answer inside it. If you get an error saying the signs are the same at both ends, the bracket is wrong, not the maths.

Finding the best setting

Root finding asks where does this hit zero. Minimising asks where is this as small as it can be — the cheapest, the fastest, the least wasteful.

A shop sells about 4,000 items a year. Ordering in big batches means fewer orders but more stock sitting in the back room; ordering little and often means the opposite. How many should they order at a time?

from scipy.optimize import minimize_scalar

def yearly_cost(batch):
    return 4000.0 / batch + 3.0 * batch

best = minimize_scalar(yearly_cost, bounds=(1, 200), method='bounded')
print(f"order {best.x:.0f} at a time")
print(f"yearly cost {best.fun:.2f}")

order 37 at a time
yearly cost 219.09

best.x is where the minimum is; best.fun is the value there. Giving bounds with method='bounded' keeps the search inside a sensible range — you cannot order half an item or a million of them, and saying so stops the search wandering somewhere meaningless.

For more than one unknown at a time, scipy.optimize.minimize takes the same shape of problem with a list of starting values instead of a single bound.

Adding up an area

quad computes a definite integral — the area under a curve between two points. You hand it a function and the two ends.

from scipy.integrate import quad
import numpy as np

area, error = quad(lambda x: np.exp(-x ** 2), 0, np.inf)
print(f"area        {area:.10f}")
print(f"uncertainty {error:.1e}")

area        0.8862269255
uncertainty 7.1e-09

np.inf is allowed as an endpoint, which is genuinely useful — this is a curve with no end, and the area under all of it is still finite.

Note the second return value again: quad tells you how far off it might be. Here, nine decimal places in. As with curve_fit, SciPy is consistently willing to tell you how much to trust it, and it is worth getting into the habit of looking.

Change over time

Some problems don't give you a formula for the answer — they give you a rule for how fast things change, and you want to know what happens. That is a differential equation, and solve_ivp steps it forward for you.

A thousand people. One of them is ill. Ill people infect well people at some rate, and recover at another. What happens over the next four months?

from scipy.integrate import solve_ivp
import numpy as np

def outbreak(day, state, infect_rate, recover_rate):
    well, ill, better = state
    people = well + ill + better
    caught = infect_rate * well * ill / people
    recovered = recover_rate * ill
    return [-caught, caught - recovered, recovered]

days = solve_ivp(outbreak, [0, 120], [999, 1, 0], args=(0.3, 0.1), max_step=1.0)

ill = days.y[1]
print(f"worst day    : day {days.t[ill.argmax()]:.0f}")
print(f"ill that day : {ill.max():.0f} people")
print(f"never ill    : {days.y[0][-1]:.0f} people")

worst day    : day 38
ill that day : 301 people
never ill    : 60 people

You never wrote a formula for "how many are ill on day 38". You wrote the three rates of change, and solve_ivp did the rest.

The pieces: your function returns the rate of change of each quantity; [0, 120] is the time range; [999, 1, 0] is how things start; args passes your constants through. Results come back as days.t (the times) and days.y (one row per quantity). max_step keeps the steps small enough that "the worst day" is a real day rather than an artefact of a long stride.

Filling in between measurements

You measured at some points and want values in between. Interpolation draws a smooth curve through what you have.

from scipy.interpolate import CubicSpline
import numpy as np

hours = np.array([0, 3, 6, 9, 12, 15, 18, 21])
outside = np.array([11.2, 10.1, 12.8, 17.4, 21.0, 22.3, 18.9, 14.2])

curve = CubicSpline(hours, outside)
print(f"at 07:30 it was about {curve(7.5):.1f} C")

minute = np.linspace(0, 21, 1000)
print(f"warmest {curve(minute).max():.1f} C at {minute[curve(minute).argmax()]:.1f} hours")

at 07:30 it was about 15.1 C
warmest 22.4 C at 14.5 hours

The warmest point it reports, 14.5 hours, is between two readings — the thermometer was never read then.

This is the section to be most careful with. Interpolation is not fitting. Fitting says I believe the world behaves like this shape, find the numbers; interpolation says join these dots smoothly and ask no questions. A spline goes exactly through every point you gave it, including the wrong ones, and beyond the ends of your data it is guessing with confidence. Use it to fill small gaps between trustworthy readings. Use curve_fit when you have a reason to believe a shape, when your readings are noisy, or when you want to know what is happening outside the range you measured.

Is this difference real?

scipy.stats is where the statistics live. Three jobs cover most of what people actually need.

A straight line through data, with the numbers that say whether to believe it:

from scipy import stats
import numpy as np

spend = np.array([1.2, 2.4, 3.1, 4.8, 5.5, 6.9, 8.2, 9.0])
sales = np.array([15.1, 19.8, 22.0, 29.4, 31.8, 38.2, 43.9, 46.1])

line = stats.linregress(spend, sales)
print(f"every extra 1 of spend brings {line.slope:.2f} more sales")
print(f"explains {line.rvalue ** 2:.1%} of the variation")
print(f"p-value {line.pvalue:.1e}")

every extra 1 of spend brings 4.07 more sales
explains 99.9% of the variation
p-value 8.0e-10

Comparing two groups. A process was changed; times were recorded before and after. Is it genuinely faster, or did it just look that way?

from scipy import stats
import numpy as np

before = np.array([12.1, 11.8, 13.0, 12.4, 11.9, 12.7, 12.2, 11.6])
after  = np.array([11.2, 10.9, 11.8, 11.1, 10.7, 11.5, 11.0, 10.8])

t, p = stats.ttest_ind(before, after)
print(f"before {before.mean():.2f}s, after {after.mean():.2f}s")
print(f"p-value {p:.5f}")

before 12.21s, after 11.12s
p-value 0.00015

A small p-value means a difference this large would rarely happen by chance alone if the change had done nothing. It does not tell you the change caused it, and it does not tell you the difference matters — a second saved may be worth nothing. It answers one narrow question, and reading more into it than that is the most common mistake in applied statistics.

Working with a distribution. Build one, then ask it anything:

from scipy import stats

scores = stats.norm(loc=100, scale=15)
print(f"above 130   : {1 - scores.cdf(130):.2%}")
print(f"top 5% from : {scores.ppf(0.95):.0f}")

above 130   : 2.28%
top 5% from : 125

cdf gives the share below a value, and ppf is its reverse — the value at a given share. Every distribution in scipy.stats answers the same methods, so swapping norm for poisson, expon or binom changes the model without changing how you ask.

Finding a rhythm hidden in a signal

Some data has a repeating pattern buried under noise you cannot see by eye. A Fourier transform reports which frequencies are present and how strongly.

Here is a reading made of a slow 7 Hz rhythm, a weaker fast 33 Hz one, and a pile of noise on top — and the transform finding both without being told they are there.

import numpy as np
from scipy.fft import rfft, rfftfreq

rate = 500
seconds = np.arange(0, 4, 1 / rate)
reading = (2.0 * np.sin(2 * np.pi * 7 * seconds)
           + 0.8 * np.sin(2 * np.pi * 33 * seconds)
           + 0.4 * np.random.default_rng(0).standard_normal(seconds.size))

strength = np.abs(rfft(reading))
hz = rfftfreq(seconds.size, 1 / rate)

for i in sorted(np.argsort(strength)[-2:]):
    print(f"{hz[i]:4.1f} Hz  strength {strength[i]:.0f}")

 7.0 Hz  strength 1993
33.0 Hz  strength 812

Both recovered exactly, and their relative strengths match how they were built — 2.0 against 0.8.

The pairing to remember is rfft with rfftfreq: the first gives strengths, the second tells you which frequency each one belongs to, and it needs to know your sampling rate to do it. Without the second you have a list of numbers with no labels. Use rfft for ordinary real-world readings; plain fft is for the complex-valued case.

When SymPy is the better tool

This is the most useful distinction on this page, and it is worth more than any function above.

SciPy computes with numbers. SymPy computes with symbols.

quad above returned 0.8862269255 — a number, correct to about nine decimal places, arrived at by clever approximation. SymPy would tell you that the same area is exactly half the square root of pi. Neither answer is better; they are answers to different questions.

Reach for SymPy when you want:

• an exact answer — a fraction that is really a third, not 0.3333333333
• a formula rather than a value, including one you can rearrange or read a meaning out of
• to do algebra: solve for x, differentiate, expand, simplify, factorise
• to check what a result is before you ever put numbers in it

Reach for SciPy when you want:

• an answer from data you measured, which has no formula behind it
• a problem with no exact answer at all — most real integrals and most equations genuinely have none, and approximating well is the only option available
• speed over many values, because it works on whole arrays at once
• anything statistical, or a signal, or a fit

A good habit combines them: work the problem out exactly in SymPy while it is still small, then hand the result to SciPy to run over your actual numbers. The SymPy lesson covers that side properly.

What else is in here

Everything below is present on your device and imports without anything extra. If a problem in one of these areas comes up, the module is already there:

Two notes on the neighbours. scipy.linalg overlaps with NumPy's own numpy.linalg and is generally the more capable of the two — it has decompositions and solvers NumPy does not. And SciPy has no drawing of its own: to see any of this as a picture, hand the numbers to Matplotlib.

Where to go next

Try changing one number in an example and re-running the cell. Give curve_fit a deliberately bad p0 and watch it struggle. Move a single reading in the spline example and see how far the curve moves. Drop the 33 Hz term out of the signal and check that the transform stops reporting it. Each of those takes a few seconds in a CodeBook cell and teaches more than reading another section.

Related

• SymPy — exact answers and algebra, the other half of this pairing
• NumPy — the arrays everything here is built on
• Matplotlib — turning these results into pictures
• CodeBook — the notebook these examples are written for
• Code Editor — running Python as a plain document
• Languages, and Where They Run — where a Python step actually executes