Curve fitting

Why this matters

Measurements become useful when they pin down the numbers in a model: the spring constant of a sample, the rate constant of a reaction, the time constant of a sensor, the Michaelis constant of an enzyme, the gravitational acceleration from a falling object. Curve fittingFinding the parameter values that make a model match the data best, along with their uncertainties. finds the parameter values that make the model match the data best, and, just as important, tells you how uncertain those values are. A fitted value without an uncertainty is only half an answer.

Core idea: least squares

For a model f(x; p) with parameters p and data points (x_i, y_i), least squares chooses p to minimise the sum of squared residuals:

S(p) = Σ (y_i - f(x_i; p))^2

If the model is linear in its parameters (a straight line, a polynomial, a sin t + b cos t), there is a direct solution with a design matrix, as you saw with np.linalg.lstsq. If not (A e^(-t/τ), Vmax S / (Km + S)), an iterative optimiser starts from a guess and walks downhill. SciPy's curve_fit handles both.

scipy.optimize.curve_fit

import numpy as np
from scipy.optimize import curve_fit

def decay(t, A, tau, C):
    return A * np.exp(-t / tau) + C

t = np.linspace(0, 10, 30)                      # seconds
rng = np.random.default_rng(0)
y = decay(t, 5.0, 2.5, 0.8) + rng.normal(0, 0.05, t.size)

p0 = [4.0, 2.0, 0.5]                            # starting guess
popt, pcov = curve_fit(decay, t, y, p0=p0)
perr = np.sqrt(np.diag(pcov))                   # 1-sigma uncertainties

for name, val, err in zip(["A", "tau", "C"], popt, perr):
    print(f"{name} = {val:.3f} ± {err:.3f}")

The pieces:

• The model function takes x first, then each parameter as a separate argument.
• p0 is the starting guess. For non-linear models a sensible guess matters: estimate A from the first point, C from the last, tau from where the curve has fallen about two thirds of the way.
• popt holds the best-fit parameters.
• pcov is the covariance matrixThe matrix pcov from a fit; its diagonal holds parameter variances and the off-diagonal terms show how parameters trade off.. Its diagonal holds the variances, so sqrt(diag(pcov)) gives standard errors. Off-diagonal terms show how parameters trade off against each other.

If you know each point's measurement uncertainty, pass it: curve_fit(f, x, y, sigma=y_err, absolute_sigma=True). Without absolute_sigma=True, SciPy rescales the uncertainties to match the scatter of the residuals.

Where the uncertainties come from

For a linear model with design matrix X (one column per parameter), the covariance of the fitted parameters is

cov(p) = s^2 (X^T X)^(-1),      s^2 = S_min / (n - p)

where s^2 is the residual variance and n - p the degrees of freedomThe number of data points minus the number of fitted parameters, n - p.. This is exactly what curve_fit reports for linear models and approximates for non-linear ones. You can compute it yourself with NumPy, which is what the exercise does.

Worked example: g from a dropped ball

A video tracker gives the height of a tossed ball every 0.05 s. The model is h(t) = h0 + v0 t - g t^2 / 2: linear in h0, v0 and g, so the design matrix has columns [1, t, -t^2/2].

X = np.column_stack([np.ones_like(t), t, -0.5 * t**2])
p, *_ = np.linalg.lstsq(X, h, rcond=None)
res = h - X @ p
s2 = res @ res / (len(t) - 3)
cov = s2 * np.linalg.inv(X.T @ X)
h0, v0, g = p
g_err = np.sqrt(cov[2, 2])

With 20 points and 1 cm of tracking noise, a typical result is g = 9.78 ± 0.08 m/s²: consistent with 9.81, and the ± tells you that a second toss could easily give 9.86.

Judging the fit

The optimiser always returns numbers. Check them:

1. Plot the residuals against x. They should look like random scatter around zero. A curve or trend means the model is wrong or incomplete.
2. Compare the residual size with your measurement uncertainty. With known errors, the reduced chi-squareThe sum of squared, error-scaled residuals divided by the degrees of freedom; near 1 means the fit matches the stated errors. χ²/ν = Σ((y - f)/σ)² / (n - p) should be near 1. Much larger: the model misses something or errors are underestimated. Much smaller: errors are overestimated.
3. Look at the uncertainties. A parameter with error larger than its value is not determined by your data.
4. Check correlations. Strongly correlated parameters (correlation near ±1) mean the data cannot separate them; consider fixing one or collecting data over a wider range.

Common mistakes

• No starting guess for a non-linear model: curve_fit defaults every parameter to 1, which may converge to nonsense or fail with "Optimal parameters not found".
• Reporting pcov diagonals as uncertainties: they are variances; take the square root.
• Fitting too many parameters for the data, which gives huge, correlated uncertainties.
• Extrapolating the fitted model far outside the measured range.
• Ignoring units: the units of a parameter follow from the model (if t is in seconds, tau is in seconds).

Quick recap

• popt, pcov = curve_fit(model, x, y, p0=guess); uncertainties are np.sqrt(np.diag(pcov)).
• For linear-in-parameter models, cov = s² (XᵀX)⁻¹ with s² = SSR / (n - p).
• Always look at residuals, the uncertainties and parameter correlations before trusting a fit.

Try it

In the Do tab you will fit the dropped-ball model with NumPy, compute the uncertainty in g from the covariance matrix, and test it on seeded noisy data.