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 fitting 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 matrix. 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 freedom. 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-square χ²/ν = Σ((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.