PYTHON / SCIPY
Optimisation and curve fitting
Fit parameterised models to data with scipy.optimize.curve_fit and minimise your own objective functions, then judge whether the result is trustworthy.
What you will learn
- Write a model as f(x, *params) and fit it with curve_fit(f, x, y, p0=...)
- Read parameter uncertainties from sqrt(np.diag(pcov)) and check residuals
- Constrain parameters with bounds=(lo, hi) so fits stay physically meaningful
- Use minimize / minimize_scalar when the objective is not a sum of residuals
Understanding Optimisation and curve fitting
Curve fitting asks a different question from interpolation. Interpolation forces a curve through every point; fitting assumes you already know the shape of the law, say y = a*exp(-b*x) + c, and searches for the parameter values that make the model as close as possible to the data in the least-squares sense. scipy.optimize.curve_fit does exactly that: it minimises the sum of squared residuals sum((y_i - f(x_i, params))**2) over the parameter vector, and hands back both the best parameters and their estimated covariance.
The mental model to hold is a landscape over parameter space whose height is that residual sum. For a linear model the landscape is a single bowl and there is one answer. For a nonlinear model like an exponential or a damped sine it can have several valleys, and curve_fit only walks downhill from where you put it. That is why p0 is not a formality: with no p0 SciPy starts every parameter at 1.0, which for a decay constant that should be 0.001, or an amplitude that should be 5000, can start you in a flat region where the solver stalls and raises RuntimeError.
The second return value, pcov, comes from the curvature of that landscape at the minimum. A narrow, steep valley in one direction means the data pins that parameter down, so sqrt(diag(pcov)) is small; a long flat trough means two parameters trade off against each other and the errors are large and strongly correlated. If you do not pass sigma, curve_fit rescales pcov using the observed residual spread, so the numbers only mean something once you have confirmed the model actually describes the data. When you are not minimising residuals at all, drop down to minimize for vector arguments or minimize_scalar for one variable, and supply the same objective you would compute by hand.
import numpy as np
from scipy.optimize import curve_fit, minimize_scalar
def decay(x, a, b, c):
return a * np.exp(-b * x) + c
x = np.linspace(0.0, 4.0, 25)
y = decay(x, 2.5, 1.3, 0.5)
popt, pcov = curve_fit(decay, x, y, p0=[1.0, 1.0, 0.0])
print("fitted a, b, c:", np.round(popt, 4))
worst = np.max(np.abs(y - decay(x, *popt)))
print("largest residual below 1e-6:", bool(worst < 1e-6))
hit = minimize_scalar(lambda t: (decay(t, *popt) - 1.0) ** 2,
bounds=(0.0, 4.0), method="bounded",
options={"xatol": 1e-10})
print("x where the curve reaches 1.0:", round(hit.x, 3))curve_fit is a downhill search over the sum of squared residuals, so the starting guess decides which minimum you land in and pcov describes how sharply that minimum is defined.
Worked examples
Minimising an objective directly
Shows minimize solving a general two-variable problem that is not a residual sum.
import numpy as np
from scipy.optimize import minimize
def cost(v):
x, y = v
return (x - 3.0) ** 2 + (y + 1.0) ** 2 + 2.0
res = minimize(cost, x0=[0.0, 0.0], method="BFGS")
print(res.success)
print(f"{res.x[0]:.3f} {res.x[1]:.3f}")
print(f"{res.fun:.4f}")Example explained
Line 1cost takes a single array argument, not separate scalars, because minimize always passes one parameter vector.
Line 2x0 is the starting point; BFGS estimates the gradient by finite differences when you do not supply jac.
Line 3res.fun is the objective value at the minimum, 2.0 here, not the distance to the optimum.
Line 4Always check res.success before using res.x, since a failed run still returns whatever point it stopped at.
Keeping parameters positive with bounds
Fits a saturating Michaelis-Menten curve where both parameters must stay above zero.
import numpy as np
from scipy.optimize import curve_fit
def mm(s, vmax, km):
return vmax * s / (km + s)
s = np.array([0.5, 1.0, 2.0, 4.0, 8.0, 16.0])
v = mm(s, 4.0, 2.0)
popt, _ = curve_fit(mm, s, v, p0=[1.0, 1.0], bounds=(0.0, np.inf))
print(f"vmax={popt[0]:.3f} km={popt[1]:.3f}")Example explained
Line 1bounds=(0.0, np.inf) applies the same lower and upper limit to every parameter; pass sequences for per-parameter limits.
Line 2Supplying bounds switches curve_fit from the Levenberg-Marquardt method to trf, so it can never step to a negative km.
Line 3A negative km would put a pole at s = -km inside the data range and produce nonsense, which the bound prevents outright.
Line 4The recovered values match the generating parameters because the data here is noise-free.
Seeing that the fit really is the minimum
Computes the residual sum at the fitted parameters and at a slightly worse slope.
import numpy as np
from scipy.optimize import curve_fit
def line(x, m, c):
return m * x + c
x = np.array([0.0, 1.0, 2.0, 3.0])
y = np.array([1.0, 3.0, 5.0, 8.0])
popt, _ = curve_fit(line, x, y, p0=[0.0, 0.0])
ssr = np.sum((y - line(x, *popt)) ** 2)
worse = np.sum((y - line(x, popt[0] + 0.1, popt[1])) ** 2)
print(f"m={popt[0]:.3f} c={popt[1]:.3f}")
print(f"ssr={ssr:.4f}")
print(f"slope +0.1 -> ssr={worse:.4f}")Example explained
Line 1The data is not exactly collinear, so the fitted line misses every point and ssr is 0.3 rather than 0.
Line 2Nudging the slope by 0.1 raises the residual sum to 0.44, which is what 'best fit' means numerically.
Line 3line(x, *popt) unpacks the fitted parameters into the model, the same call pattern curve_fit uses internally.
Line 4For this linear model the landscape is a single bowl, so p0=[0, 0] reaches the same answer as any other start.
Important notes
curve_fit is a thin wrapper over least_squares; it uses method 'lm' when unbounded and 'trf' when you pass bounds, so error messages and convergence behaviour change once bounds appear.
If parameters differ by many orders of magnitude, rescale them (fit b in units of 1/ms, not 1/s) or pass x_scale, because the shared step tolerances otherwise favour the large parameter.
Common mistakes
Omitting p0 for a nonlinear model: every parameter starts at 1.0, and for a rate like 1e-4 or an amplitude like 1e4 the solver either stalls or raises 'Optimal parameters not found: Number of calls to function has reached maxfev'.
Writing the model as f(params, x) instead of f(x, a, b, c): curve_fit inspects the signature to count parameters, so you get a TypeError or a fit over the wrong number of unknowns.
Quoting sqrt(diag(pcov)) as an error bar without looking at the residuals: with no sigma given, pcov is scaled by the residual spread, so a badly shaped model can still yield small, entirely meaningless uncertainties.
Try it yourself
Change, predict, then run
Build x = np.linspace(0, 5, 60) and y = 3*np.exp(-0.7*x)*np.cos(2*x) plus noise from np.random.default_rng(0).normal(0, 0.05, 60), then fit the model f(x, A, k, w) = A*np.exp(-k*x)*np.cos(w*x) with p0=[1, 1, 1] and print the parameters together with sqrt(np.diag(pcov)).
Open the Python workspaceCheck your understanding
You run curve_fit twice on identical data, changing only p0, and get two noticeably different parameter sets. What is the most likely explanation?
- The sum of squared residuals has more than one local minimum, and the solver settles in whichever one lies downhill from the start
- curve_fit perturbs the input data slightly on each call, so repeated runs never agree
- p0 acts as a prior and is blended into the reported parameters
- pcov is recomputed from p0, and that recomputation shifts the parameters
Show answer
Nonlinear least squares is a downhill search over a residual landscape that can have several valleys, so the starting point selects the valley. Option 3 is tempting because p0 clearly influences the result, but p0 is only a starting point and is discarded once the solver converges; it is never averaged into the answer, and curve_fit is fully deterministic for fixed inputs.