`vp.optimize.curve_fit`

Fit a model's parameters to data by least squares and estimate their covariance, as SciPy's curve_fit.

Part of the distribution vipyr; it runs from a checkout of the repository. Planned on PyPI as pip install vipyr.

vp.optimize.curve_fit(f, xdata, ydata, p0) minimizes $\sum_{i} ((f (x_{i}) - y_{i}) / σ_{i})^{2}$ with MINPACK's Levenberg-Marquardt and returns the parameters and their covariance. As in SciPy, popt, pcov = vp.optimize.curve_fit(...) unpacks the two; the result also carries the standard errors perr, the residuals, and the evaluation counts. The defaults are those of scipy.optimize.leastsq, which SciPy's curve_fit calls for method="lm".

The model is evaluated one point at a time: f(x, b1, ..., bk), SciPy's form, or f(x, p) with p an array (a Viv def f(x: Num, p: Arr) -> Num). Its Jacobian is exact: reverse-mode grad of the model at each point, so a Jacobian costs about one model evaluation per point. With a Python model written for SciPy (np.exp, **, arithmetic), py2viv translates it and the call is the SciPy call with vp.optimize in front.

The covariance is SciPy's: the pseudo-inverse of $J^{T} J$ for the Jacobian $J$ of the weighted residuals at the solution, multiplied by the reduced chi-square $rss / (m - n)$ unless absolute_sigma=True. A rank-deficient $J$ gives an infinite covariance and an OptimizeWarning. A fit that stops without converging raises RuntimeError, as SciPy's does.

Signature

vp.optimize.curve_fit(f, xdata, ydata, p0=None, *, sigma=None, absolute_sigma=False, jac=None, ftol=1.49012e-08, xtol=1.49012e-08, gtol=0.0, maxfev=None, epsfcn=None, diag=None, fast=False)

Parameters

f : str | function, required: The model at one point, f(x, b1, ..., bk) or f(x, p) with p an array: 'file.viv:def', a def of a module from vp.load, or a Python function in py2viv's subset.
xdata : ndarray[float64], required: The independent variable, one number per point (length m).
ydata : ndarray[float64], required: The observations, length m.
p0 : ndarray[float64], default None: The starting parameters; None: ones, one per parameter of f.
sigma : ndarray[float64], default None: The standard deviation of each observation; the fit minimizes sum(((f(x) - y) / sigma)^2).
absolute_sigma : bool, default False: True: pcov is (J^T J)^-1 with sigma taken as absolute; False: scaled by the reduced chi-square rss / (m - n), as SciPy.
jac : str, default None: None: exact, reverse-mode grad of the model at each point; 'fd': MINPACK's forward differences (what SciPy's curve_fit uses).
ftol : float, default 1.49012e-08: Leastsq's ftol.
xtol : float, default 1.49012e-08: Leastsq's xtol.
gtol : float, default 0.0: Leastsq's gtol.
maxfev : int, default None: The residual evaluations allowed; None: 100 (n + 1), or 200 (n + 1) with jac='fd', as leastsq.
epsfcn : float, default None: With jac='fd', the relative step squared; None: eps.
diag : ndarray[float64], default None: MINPACK's variable scales; None: the Jacobian's column norms.
fast : bool, default False: Compile with the fast flags (fused multiply-add): faster, last bits may differ.

Returns

A frozen CurveFitResult. Iterating it gives the fields below in order; it also has params, tests, log and check().

popt : ndarray[float64]: The fitted parameters.
pcov : ndarray[float64]: Their covariance, n by n (infinite if the Jacobian is rank deficient).
perr : ndarray[float64]: The standard errors, sqrt(diag(pcov)).
fun : ndarray[float64]: The weighted residuals (f(x) - y) / sigma at popt.
rss : float: Their sum of squares.
nfev : int: Model evaluations over all points, counted as leastsq counts them.
njev : int: Jacobian evaluations.
info : int: MINPACK's info, 1 to 4.
message : str: Info in words.
success : bool: Always True: a fit that does not converge raises RuntimeError.

Raises

ValueError: A tolerance below machine epsilon, fewer residuals than variables, or an array of the wrong shape; vipyr._optimize.OptimizeProblemError (a ValueError) when the function or jac cannot be made into a Viv problem, or py2viv refuses a Python function.
vipyr.VivCompileError: A SyntaxError: the function's .viv file (or py2viv's translation) does not compile.
RuntimeError: Lmder stopped without converging (info 0 or 5 to 8), as SciPy's 'Optimal parameters not found'.
TypeError: An argument of the wrong type, or an unknown keyword.
vipyr.VivBuildError: An ImportError: no C compiler, or the program's C does not build.

Notes

Against NIST. The 26 single-predictor problems of NIST's Statistical Reference Datasets for nonlinear regression, from both certified starting points, are fitted in lib/minpack/bench_optimize.py; the table in lib/minpack/bench_optimize.json gives each problem's certified digits matched (the log relative error of the worst parameter, and of the worst standard deviation), its evaluation counts and SciPy's. The gallery's NIST page draws one of them.

Limits. One independent variable per point; for several (NIST's Nelson, xdata of shape (k, m) in SciPy) write the residuals for vp.optimize.least_squares and pass the columns in args. sigma is a vector of standard deviations; a full covariance matrix of the data is not taken. No bounds (as SciPy's method="lm").

References

NIST Statistical Reference Datasets, nonlinear regression: 27 problems with certified parameter values and standard deviations, in lib/minpack/data. https://www.itl.nist.gov/div898/strd/nls/nls_main.shtml
J. J. Moré, The Levenberg-Marquardt algorithm: implementation and theory, in Numerical Analysis (Dundee 1977), Lecture Notes in Mathematics 630, 105-116, Springer (1978). doi:10.1007/BFb0067700
P. Virtanen et al., SciPy 1.0: fundamental algorithms for scientific computing in Python, Nat. Methods 17, 261-272 (2020): scipy.optimize.curve_fit, whose interface and covariance this follows. doi:10.1038/s41592-019-0686-2

Examples

import numpy as np
import vipyr as vp

# NIST StRD Misra1a, from its first certified start
x = [77.6, 114.9, 141.1, 190.8, 239.9, 289.0, 332.8, 378.4, 434.8, 477.3, 536.8, 593.1, 689.1, 760.0]
y = [10.07, 14.73, 17.94, 23.93, 29.61, 35.18, 40.02, 44.82, 50.76, 55.05, 61.01, 66.40, 75.47, 81.78]

def misra1a(x, b1, b2):
    return b1 * (1.0 - np.exp(-b2 * x))

popt, pcov = vp.optimize.curve_fit(misra1a, x, y, [500.0, 1e-4])
print(popt, np.sqrt(np.diag(pcov)))      # NIST: 238.94212918, 5.5015643181e-4; sd 2.7070075241, 7.2668688436e-6
r = vp.optimize.curve_fit(misra1a, x, y, [500.0, 1e-4])
assert r.check()

[2.38942129e+02 5.50156432e-04] [2.70700752e+00 7.26686884e-06]

A model written as a Viv def taking its parameters as one array, with measurement errors.

import tempfile
from pathlib import Path
import numpy as np
import vipyr as vp

src = Path(tempfile.mkdtemp()) / "peak.viv"
src.write_text("""
def peak(x: Num, p: Arr) -> Num {
    return p[0] * exp(-0.5 * ((x - p[1]) / p[2]) ^ 2) + p[3]
}
""")
x = np.linspace(-3.0, 3.0, 61)
y = 2.0 * np.exp(-0.5 * ((x - 0.4) / 0.7) ** 2) + 0.1 + 0.02 * np.cos(11.0 * x)
r = vp.optimize.curve_fit("%s:peak" % src, x, y, [1.0, 0.0, 1.0, 0.0], sigma=np.full(61, 0.02),
                          absolute_sigma=True)
print(np.round(r.popt, 5), np.round(r.perr, 5), r.nfev)

[1.99903 0.4     0.69935 0.10134] [0.00754 0.00281 0.00346 0.00409] 7

Each example ran as written when these pages were generated, and the output shown is what it printed.

Checked against

Reference: NIST's 26 single-predictor StRD nonlinear regression problems from both certified starts: the certified parameters and standard deviations, and SciPy 1.12's curve_fit counts (lib/minpack/bench_optimize.py, lib/minpack/bench_optimize.json); Misra1a in lib/minpack/fit.viv.

Credits

A port of MINPACK (Argonne National Laboratory), under MINPACK's BSD-style copyright notice (1999, University of Chicago), which the port and the Fortran in lib/minpack/ref both carry. This product includes software developed by the University of Chicago, as Operator of Argonne National Laboratory.