This site is version 0.0.1. Most of the code on this site is drafted by an agent and needs review. I apologize for the AI-speak, we are working on rewriting it!
viv-lang.org
Python

Thurber's electron mobility, a NIST certified fit

NIST's Statistical Reference Datasets give nonlinear regression problems with parameters certified to eleven digits. Thurber (higher difficulty) is a rational function of degree three over three, seven parameters, fitted to 37 measurements of electron mobility against the log of density. vp.optimize.curve_fit fits it from both of NIST's starting points with MINPACK's Levenberg-Marquardt and an exact Jacobian, and the table compares every parameter and standard deviation with the certified value.

The model and the data

The model is NIST's, written for SciPy in lib/minpack/nist_strd.py:

def rational3(x, b1, b2, b3, b4, b5, b6, b7):
    return (b1 + b2 * x + b3 * x ** 2 + b4 * x ** 3) / (1.0 + b5 * x + b6 * x ** 2 + b7 * x ** 3)

py2viv translates it to Viv once; the driver differentiates it with reverse-mode grad at each point, so the Jacobian has no step-size error.

import json
import sys

import numpy as np
import vipyr as vp
from vipyr.figure import figure

sys.path.insert(0, "lib/minpack")
import nist_strd as N  # noqa: E402

p = N.load("Thurber")
print(
    "Thurber: %d observations, %d parameters, %s difficulty"
    % (len(p.y), len(p.certified), p.difficulty)
)
fits = [
    vp.optimize.curve_fit(N.rational3, p.x, p.y, start) for start in p.starts
]
for k, r in enumerate(fits, 1):
    print(
        "start %d: %3d evaluations, %3d Jacobians, MINPACK info %d, residual "
        "sum of squares %.10e" % (k, r.nfev, r.njev, r.info, r.rss)
    )
print("certified residual sum of squares          %.10e" % p.rss)
Thurber: 37 observations, 7 parameters, higher difficulty
start 1:  33 evaluations,  27 Jacobians, MINPACK info 1, residual sum of squares 5.6427082614e+03
start 2:  18 evaluations,  17 Jacobians, MINPACK info 1, residual sum of squares 5.6427082568e+03
certified residual sum of squares          5.6427082397e+03

Against the certified values

The digits that agree are NIST's log relative error, LRE = -log10(|value - certified| / |certified|). With leastsq's defaults lmder stops when the sum of squares falls by less than 1.5 × 10−8 relative (MINPACK info 1). Thurber's parameters are strongly correlated, so a sum of squares right to nine digits leaves them right to about five. NIST counts four or more digits as a successful fit.

r = fits[0]
lre, lre_sd = N.lre(r.popt, p.certified), N.lre(r.perr, p.cert_sd)
print(
    "     %22s %22s %5s   %16s %16s %5s"
    % ("certified", "vp.optimize", "LRE", "certified sd", "vp sd", "LRE")
)
for i in range(len(p.certified)):
    print(
        "b%d  %22.10e %22.10e %5.1f   %16.10e %16.10e %5.1f"
        % (
            i + 1,
            p.certified[i],
            r.popt[i],
            lre[i],
            p.cert_sd[i],
            r.perr[i],
            lre_sd[i],
        )
    )
                  certified            vp.optimize   LRE       certified sd            vp sd   LRE
b1        1.2881396800e+03       1.2881397251e+03   7.5   4.6647963344e+00 4.6648024297e+00   5.9
b2        1.4910792535e+03       1.4910709812e+03   5.3   3.9571156086e+01 3.9576874779e+01   3.8
b3        5.8323836877e+02       5.8323259612e+02   5.0   2.8698696102e+01 2.8702835321e+01   3.8
b4        7.5416644291e+01       7.5415515842e+01   4.8   5.5675370270e+00 5.5683609863e+00   3.8
b5        9.6629502864e-01       9.6628945360e-01   5.2   3.1333340687e-02 3.1337226644e-02   3.9
b6        3.9797285797e-01       3.9797031629e-01   5.2   1.4984928198e-02 1.4986795863e-02   3.9
b7        4.9727297349e-02       4.9725635312e-02   4.5   6.5842344623e-03 6.5852093944e-03   3.8

Tighter tolerances, more digits

ftol and xtol are leastsq's. Each two decades smaller buys about one certified digit, for a handful of extra evaluations.

for tol in (1e-10, 1e-12, 1e-14):
    for k, start in enumerate(p.starts, 1):
        t = vp.optimize.curve_fit(
            N.rational3, p.x, p.y, start, ftol=tol, xtol=tol
        )
        print(
            "ftol = xtol = %.0e, start %d: %2d evaluations, LRE of "
            "parameters %.1f, of standard deviations %.1f"
            % (
                tol,
                k,
                t.nfev,
                N.lre(t.popt, p.certified).min(),
                N.lre(t.perr, p.cert_sd).min(),
            )
        )
ftol = xtol = 1e-10, start 1: 40 evaluations, LRE of parameters 5.7, of standard deviations 5.0
ftol = xtol = 1e-10, start 2: 24 evaluations, LRE of parameters 5.6, of standard deviations 4.9
ftol = xtol = 1e-12, start 1: 46 evaluations, LRE of parameters 6.7, of standard deviations 6.1
ftol = xtol = 1e-12, start 2: 30 evaluations, LRE of parameters 6.6, of standard deviations 5.9
ftol = xtol = 1e-14, start 1: 51 evaluations, LRE of parameters 7.6, of standard deviations 6.9
ftol = xtol = 1e-14, start 2: 38 evaluations, LRE of parameters 7.9, of standard deviations 7.3

Against SciPy

SciPy's curve_fit calls the same Fortran (lmdif, with differences), so with jac="fd" the two take the same steps. The counts and the certified digits of each are from lib/minpack/bench_optimize.py, over all 27 NIST problems; the times are for one warm fit and depend on the machine and its load, recorded with them.

bench = json.load(open("lib/minpack/bench_optimize.json"))
print(
    "measured on",
    bench["machine"],
    "| load average",
    bench["load_average_start"],
    "| SciPy",
    bench["scipy"],
)
print(
    "%-9s %5s  %22s  %22s  %22s"
    % (
        "problem",
        "start",
        "vp (grad): LRE nfev ms",
        "vp (fd): LRE nfev ms",
        "SciPy (fd): LRE nfev ms",
    )
)


def cell(c):
    return (
        "%4.1f %5d %7.2f" % (c["lre"], c["nfev"], 1e3 * c["s"])
        if c["ok"]
        else "%22s" % "did not converge"
    )


for row in bench["nist"]:
    if row["difficulty"] == "higher":
        print(
            "%-9s %5d  %22s  %22s  %22s"
            % (
                row["problem"],
                row["start"],
                cell(row["vp"]),
                cell(row["vp_fd"]),
                cell(row["scipy"]),
            )
        )
measured on Apple M2 Pro, Darwin | load average 12.06 | SciPy 1.12.0
problem   start  vp (grad): LRE nfev ms    vp (fd): LRE nfev ms  SciPy (fd): LRE nfev ms
MGH09         1       4.1   494    1.32       4.0  1757    0.98        did not converge
MGH09         2       4.0    18    0.28       4.0    82    0.27       4.0    82    0.39
Thurber       1       4.5    33    0.80       4.6   230    1.35       4.5   222    2.04
Thurber       2       4.5    18    0.59       4.5   137    0.89       4.5   137    1.28
BoxBOD        1       0.0    10    0.46       0.0    10    0.24       0.0    10    0.07
BoxBOD        2       5.1     9    0.24       5.1    25    0.23       5.1    25    0.11
Rat42         1       6.8    10    0.27       6.8    34    0.29       6.8    34    0.37
Rat42         2       7.1     6    0.79       7.1    21    0.83       7.1    21    0.11
MGH10         1       9.8   286    1.08       6.9  1008    0.92        did not converge
MGH10         2       9.0   126    0.63       7.3   474    0.52       7.9   474    1.36
Eckerle4      1       7.2    18    0.41       7.2    63    0.33       7.2    63    0.28
Eckerle4      2       7.2     7    0.30       7.2    25    0.28       7.2    25    0.14
Rat43         1       5.1    26    0.42       5.1   102    0.41       5.1   102    0.46
Rat43         2       4.9     7    0.30       4.9    31    0.30       4.9    31    0.17
Bennett5      1        did not converge       5.0  3000   26.56        did not converge
Bennett5      2       4.1   203    8.74       5.4   790    7.29       5.7   787    4.15

The fit and its residuals

The curve through the data, and the residuals y - f(x) against x. NIST's residual standard deviation for Thurber is 13.71; the fit's is printed below it.

xs = np.linspace(p.x.min(), p.x.max(), 200)
fig = figure(
    "Thurber: electron mobility against log density", "log density", "mobility"
)
fig.scatter(p.x, p.y, "NIST data")
fig.plot(xs, N.rational3(xs, *r.popt), "vp.optimize fit")
fig.style("vp.optimize fit", "s2 nomarker")

res = p.y - N.rational3(p.x, *r.popt)
print(
    "residual standard deviation %.4f (NIST, certified: 1.3714600784E+01)"
    % np.sqrt(res @ res / (len(p.y) - len(p.certified)))
)
fig2 = figure("Residuals of the fit", "log density", "y - f(x)")
fig2.scatter(p.x, res, "residual")
fig2.hline("zero", 0.0)
residual standard deviation 13.7146 (NIST, certified: 1.3714600784E+01)
Thurber: electron mobility against log densityNIST datavp.optimize fit050010001500−3−2−1012log densitymobility
Residuals of the fitresidual−20020−3−2−1012log densityy - f(x)

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.

Source

gallery/math/plot_nist_thurber_fit.py, run as written when these pages were built. The output and the figures above are from that run.

The Viv program it calls is posts/2026-09-25-fitting-electron-mobility-with-minpack.viv: its atlas, the defs and the calls between them, or open it in the workbench.