SciPy curve_fit in Python: maxfev, bounds and p0

SciPy curve fit takes a model function and some data, then finds the parameter values that make the model match. It lives in SciPy‘s optimize module:

from scipy.optimize import curve_fit

popt, pcov = curve_fit(model, xdata, ydata, p0=[1, 1, 1])

popt comes back with the fitted parameters and pcov with their covariance. The one thing that goes wrong more than anything else is RuntimeError: Optimal parameters not found, which has its own section below.

Everything shown here was produced on Python 3.12.5, SciPy 1.18.1, NumPy 2.5.3.

A basic curve_fit example

The model function takes x first and the parameters after it:

import numpy as np
from scipy.optimize import curve_fit

# the model: its FIRST argument is x, the rest are the parameters to find
def decay(x, a, b, c):
    return a * np.exp(-b * x) + c

rng = np.random.default_rng(0)
x = np.linspace(0, 4, 50)
y = decay(x, 2.5, 1.3, 0.5) + rng.normal(0, 0.05, x.size)

popt, pcov = curve_fit(decay, x, y)

print("fitted a, b, c ->", np.round(popt, 4))
print("true values    -> [2.5    1.3    0.5   ]")
print()

errors = np.sqrt(np.diag(pcov))
print("standard errors ->", np.round(errors, 4))
print()
print("popt holds the parameters; pcov is the covariance matrix,")
print("and the square root of its diagonal gives the one-sigma errors.")

Output:

fitted a, b, c -> [2.5038 1.3844 0.5336]
true values    -> [2.5    1.3    0.5   ]

standard errors -> [0.026  0.0303 0.0107]

popt holds the parameters; pcov is the covariance matrix,
and the square root of its diagonal gives the one-sigma errors.

SciPy works out how many parameters to fit by counting the arguments of your function, so nothing needs declaring. Recovering 2.5, 1.3 and 0.5 from noisy data is the whole job.

The covariance matrix is the part people skip. Its diagonal holds the variances, so the square root gives the one-sigma uncertainty on each parameter.

RuntimeError: Optimal parameters not found

This is the error nearly everyone meets:

import numpy as np
from scipy.optimize import curve_fit

def peak(x, a, mu, sigma, offset):
    return a * np.exp(-(x - mu) ** 2 / (2 * sigma ** 2)) + offset

rng = np.random.default_rng(1)
x = np.linspace(0, 20, 300)
y = peak(x, 5, 12.0, 1.2, 0.4) + rng.normal(0, 0.05, x.size)

for budget in (5, 50):
    try:
        curve_fit(peak, x, y, p0=[1, 1, 1, 1], maxfev=budget)
    except RuntimeError as err:
        print(f"maxfev={budget:<4} -> RuntimeError: {err}")

print()
print("the number in the message is whatever budget was in force.")
print("when you did not set one, it is the default described below.")

Output:

maxfev=5    -> RuntimeError: Optimal parameters not found: Number of calls to function has reached maxfev = 5.
maxfev=50   -> RuntimeError: Optimal parameters not found: Number of calls to function has reached maxfev = 50.

the number in the message is whatever budget was in force.
when you did not set one, it is the default described below.

maxfev is a budget: the maximum number of times SciPy may call your model function while searching. When the search uses it all without settling, you get this error.

The message is not saying your data is unfittable. It is saying the solver ran out of attempts before it got there.

Why the message says maxfev = 600, 800 or 1000

Those numbers are not arbitrary. curve_fit passes the work to leastsq, which sets the budget itself when you do not:

print("curve_fit hands the work to leastsq, which sets the budget like this:")
print()
print("    if Dfun is None:")
print("        maxfev = 200 * (n + 1)")
print()
print("n is the number of PARAMETERS, and curve_fit never supplies Dfun,")
print("so the default budget is:")
print()
for n in (2, 3, 4, 5):
    print(f"  {n} parameters -> 200 * ({n} + 1) = {200 * (n + 1)}")

print()
print("so the number in your error message tells you how many parameters")
print("your model has: 600 means two, 800 means three, 1000 means four.")

Output:

curve_fit hands the work to leastsq, which sets the budget like this:

    if Dfun is None:
        maxfev = 200 * (n + 1)

n is the number of PARAMETERS, and curve_fit never supplies Dfun,
so the default budget is:

  2 parameters -> 200 * (2 + 1) = 600
  3 parameters -> 200 * (3 + 1) = 800
  4 parameters -> 200 * (4 + 1) = 1000
  5 parameters -> 200 * (5 + 1) = 1200

so the number in your error message tells you how many parameters
your model has: 600 means two, 800 means three, 1000 means four.

So the number in your error message tells you something useful straight away. A message saying 600 means your model has two parameters, 800 means three, and 1000 means four.

If that number does not match the model you thought you were fitting, the real problem is the function signature rather than the budget.

Why raising maxfev is usually the wrong fix

The common advice is to set maxfev to something huge. Here is what that actually buys you:

import numpy as np
from scipy.optimize import curve_fit
import warnings
warnings.simplefilter("ignore")

def wave(x, amplitude, frequency, phase):
    return amplitude * np.sin(frequency * x + phase)

rng = np.random.default_rng(3)
x = np.linspace(0, 10, 500)
y = wave(x, 4.0, 9.7, 0.6) + rng.normal(0, 0.2, x.size)

print("true parameters      -> [4.0   9.7   0.6 ]")
print()

lazy, _ = curve_fit(wave, x, y, p0=[1.0, 1.0, 0.0], maxfev=100000)
print("p0=[1, 1, 0], maxfev raised to 100000:")
print("  result ->", np.round(lazy, 3))
print("  no error was raised, and the answer is nonsense")
print()

good, _ = curve_fit(wave, x, y, p0=[4.0, 10.0, 0.5])
print("p0=[4, 10, 0.5], default maxfev:")
print("  result ->", np.round(good, 3))
print("  correct, and it needed no extra budget at all")

Output:

true parameters      -> [4.0   9.7   0.6 ]

p0=[1, 1, 0], maxfev raised to 100000:
  result -> [-0.15   1.237 -1.324]
  no error was raised, and the answer is nonsense

p0=[4, 10, 0.5], default maxfev:
  result -> [3.999 9.7   0.595]
  correct, and it needed no extra budget at all

With the budget raised, the fit finished without complaint and returned values nowhere near the truth. The error had been the honest outcome, and silencing it produced a worse one.

With a sensible starting guess, the same data fitted correctly inside the default budget. Running out of evaluations is nearly always a symptom of a bad p0, not of a budget that is too small.

Raise maxfev when the fit is genuinely converging but slowly, for instance with many parameters or an expensive model. Fix p0 in every other case.

How to pick a good initial guess

You can usually read the starting values off the data rather than guess them:

import numpy as np
from scipy.optimize import curve_fit

def wave(x, amplitude, frequency, phase):
    return amplitude * np.sin(frequency * x + phase)

rng = np.random.default_rng(3)
x = np.linspace(0, 10, 500)
y = wave(x, 4.0, 9.7, 0.6) + rng.normal(0, 0.2, x.size)

# read the starting values off the data instead of guessing
amplitude_guess = (y.max() - y.min()) / 2
spacing = x[1] - x[0]
spectrum = np.fft.rfft(y - y.mean())
peak_bin = np.argmax(np.abs(spectrum))
frequency_guess = 2 * np.pi * np.fft.rfftfreq(y.size, spacing)[peak_bin]

print("amplitude from the range   ->", round(amplitude_guess, 3))
print("frequency from an FFT peak ->", round(frequency_guess, 3))
print()

popt, _ = curve_fit(wave, x, y, p0=[amplitude_guess, frequency_guess, 0.0])
print("fitted ->", np.round(popt, 3))
print("true   -> [4.0   9.7   0.6 ]")

Output:

amplitude from the range   -> 4.41
frequency from an FFT peak -> 9.406

fitted -> [-3.999  9.7   -2.547]
true   -> [4.0   9.7   0.6 ]

Amplitude comes from the range of the data, and an FFT peak gives the dominant frequency. Both are a few lines, and both beat trying numbers by hand.

The same idea works generally: an exponential rate from the time to halve, an offset from the mean, a midpoint from where the data crosses halfway. For harder problems, scipy.optimize.minimize lets you drive the search yourself.

Setting bounds on the parameters

Bounds keep the search inside values that make physical sense:

import numpy as np
from scipy.optimize import curve_fit

def decay(x, a, b, c):
    return a * np.exp(-b * x) + c

rng = np.random.default_rng(0)
x = np.linspace(0, 4, 50)
y = decay(x, 2.5, 1.3, 0.5) + rng.normal(0, 0.05, x.size)

# keep every parameter positive, and cap them
popt, _ = curve_fit(decay, x, y, bounds=([0, 0, 0], [10, 5, 2]))
print("bounded fit ->", np.round(popt, 4))
print()

print("the guess has to sit inside the bounds:")
try:
    curve_fit(decay, x, y, p0=[100, 100, 100], bounds=([0, 0, 0], [10, 5, 2]))
except ValueError as err:
    print("  ValueError:", err)

print()
print("bounds also switch the algorithm, which changes the error wording:")
try:
    curve_fit(decay, x, y, bounds=([0, 0, 0], [10, 5, 2]), maxfev=5)
except RuntimeError as err:
    print("  RuntimeError:", err)

Output:

bounded fit -> [2.5038 1.3844 0.5336]

the guess has to sit inside the bounds:
  ValueError: Initial guess is outside of provided bounds

bounds also switch the algorithm, which changes the error wording:
  RuntimeError: Optimal parameters not found: The maximum number of function evaluations is exceeded.

Two details matter here. The starting guess has to lie inside the bounds, otherwise you get ValueError: Initial guess is outside of provided bounds before any fitting happens.

And adding bounds switches SciPy from Levenberg-Marquardt to a trust region method, which words its failure differently. Seeing “The maximum number of function evaluations is exceeded” instead of the maxfev message simply means bounds were in play.

Weighting points with sigma and absolute_sigma

sigma tells the fit how much to trust each point, and the flag beside it changes what the errors mean:

import numpy as np
from scipy.optimize import curve_fit

def decay(x, a, b, c):
    return a * np.exp(-b * x) + c

rng = np.random.default_rng(0)
x = np.linspace(0, 4, 50)
y = decay(x, 2.5, 1.3, 0.5) + rng.normal(0, 0.05, x.size)
measurement_error = np.full(x.size, 0.05)

_, pcov_rel = curve_fit(decay, x, y, sigma=measurement_error)
_, pcov_abs = curve_fit(decay, x, y, sigma=measurement_error, absolute_sigma=True)

print("errors with sigma given, absolute_sigma left off:")
print(" ", np.round(np.sqrt(np.diag(pcov_rel)), 4))
print("errors with absolute_sigma=True:")
print(" ", np.round(np.sqrt(np.diag(pcov_abs)), 4))
print()
print("Left off, SciPy treats sigma as RELATIVE weights and rescales the")
print("covariance by the residuals. Set to True, it trusts your error bars.")
print("Use True when sigma holds real measurement uncertainties.")

Output:

errors with sigma given, absolute_sigma left off:
  [0.026  0.0303 0.0107]
errors with absolute_sigma=True:
  [0.0303 0.0353 0.0124]

Left off, SciPy treats sigma as RELATIVE weights and rescales the
covariance by the residuals. Set to True, it trusts your error bars.
Use True when sigma holds real measurement uncertainties.

Left alone, SciPy treats sigma as relative weights and rescales the covariance using the residuals it found. That is right when your error bars express relative confidence only.

With absolute_sigma=True the values are taken as real measurement uncertainties and the covariance is not rescaled. If your sigma came from instrument specifications, that is the one you want.

When pcov comes back as infinity

An infinite covariance is not a failed fit, it is an unconstrained one:

import numpy as np
from scipy.optimize import curve_fit
import warnings

def line(x, slope, intercept):
    return slope * x + intercept

print("two data points and two parameters leaves nothing to estimate from:")
with warnings.catch_warnings(record=True) as caught:
    warnings.simplefilter("always")
    popt, pcov = curve_fit(line, np.array([1.0, 2.0]), np.array([1.0, 2.0]))
    print("  popt ->", popt)
    print("  pcov all infinite? ->", bool(np.all(np.isinf(pcov))))
    for w in caught:
        print(f"  {w.category.__name__}: {w.message}")

print()
print("an infinite pcov means the fit is unconstrained, not that it failed.")
print("Add data points, remove a parameter, or apply bounds.")
print()

def decay(x, a, b, c):
    return a * np.exp(-b * x) + c

rng = np.random.default_rng(0)
x = np.linspace(0, 4, 50)
y = decay(x, 2.5, 1.3, 0.5) + rng.normal(0, 0.05, x.size)
y[5] = np.nan

try:
    curve_fit(decay, x, y)
except ValueError as err:
    print("a NaN anywhere in ydata -> ValueError:", err)

print("nan_policy='omit' ->", np.round(curve_fit(decay, x, y, nan_policy="omit")[0], 4))

Output:

two data points and two parameters leaves nothing to estimate from:
  popt -> [1. 0.]
  pcov all infinite? -> True
  OptimizeWarning: Covariance of the parameters could not be estimated

an infinite pcov means the fit is unconstrained, not that it failed.
Add data points, remove a parameter, or apply bounds.

a NaN anywhere in ydata -> ValueError: array must not contain infs or NaNs
nan_policy='omit' -> [2.5004 1.3867 0.5343]

Two points and two parameters gives a perfect fit with nothing left over, so there is no information to estimate uncertainty from. SciPy warns and fills pcov with infinity.

The fix is more data points, fewer parameters, or bounds. A NaN anywhere in the data raises instead, and nan_policy="omit" drops those points for you.

Checking whether the fit is any good

popt always returns something, so judge it yourself:

import numpy as np
from scipy.optimize import curve_fit

def decay(x, a, b, c):
    return a * np.exp(-b * x) + c

rng = np.random.default_rng(0)
x = np.linspace(0, 4, 50)
y = decay(x, 2.5, 1.3, 0.5) + rng.normal(0, 0.05, x.size)

popt, pcov = curve_fit(decay, x, y)
residuals = y - decay(x, *popt)

ss_res = np.sum(residuals ** 2)
ss_tot = np.sum((y - np.mean(y)) ** 2)
r_squared = 1 - ss_res / ss_tot

print("R^2                 ->", round(r_squared, 6))
print("residual std dev    ->", round(residuals.std(), 4))
print("noise actually added-> 0.05")
print()
print("parameter +/- one standard error:")
for name, value, error in zip("abc", popt, np.sqrt(np.diag(pcov))):
    print(f"  {name} = {value:.4f} +/- {error:.4f}")

Output:

R^2                 -> 0.995718
residual std dev    -> 0.0416
noise actually added-> 0.05

parameter +/- one standard error:
  a = 2.5038 +/- 0.0260
  b = 1.3844 +/- 0.0303
  c = 0.5336 +/- 0.0107

R squared near 1 and a residual spread close to the known noise level are both good signs. The plot above shows the other half of the check: residuals scattered around zero with no pattern in them.

Curved or fanning residuals mean the model shape is wrong, and no amount of refitting will help. Comparing against a straight line through simple NumPy arithmetic is often a useful baseline.

Common curve_fit errors at a glance

Message	Cause	Fix
Number of calls to function has reached maxfev = 800	Search ran out of evaluations, almost always from a poor p0	Fix p0 first; raise maxfev only if it was converging
The maximum number of function evaluations is exceeded	Same thing, but on a bounded fit	As above; bounds use a different solver
Initial guess is outside of provided bounds	p0 falls outside bounds	Move p0 inside, or widen the bounds
array must not contain infs or NaNs	Missing values in the data	nan_policy="omit", or clean the arrays first
Covariance of the parameters could not be estimated	Too few points for the number of parameters	Add data, drop a parameter, or add bounds
pcov full of inf	The fit is unconstrained	Same as above

Other SciPy guides that go with this one:

• SciPy in Python: what it is and how to use it
• scipy.optimize.minimize
• scipy.optimize.leastsq
• SciPy interpolate
• scipy.stats.norm
• scipy.optimize.linprog

Frequently asked questions