Electron mobility fitted to NIST's certified digits

Adam Zettel · September 25, 2026 · written with an agent

b min i \sum (\frac{b _{1} + b _{2} x _{i} + b _{3} x _{i}^{2} + b _{4} x _{i}^{3}}{1 + b _{5} x _{i} + b _{6} x _{i}^{2} + b _{7} x _{i}^{3}} - y_{i})^{2}

A seven-parameter rational model fitted to NIST's semiconductor electron-mobility data by lmder, the Levenberg-Marquardt routine of MINPACK, translated to Viv step by step. Every trial step matches a gfortran build of the original Fortran to the last bit, the answer matches scipy.optimize.leastsq, and the parameters land within 0.001 standard deviations of NIST's certified values.

Open in the workbench Download .vivOutputs below are from the site build’s run; Run all runs it again in your browser.

Look inside this program

Atlas
The program’s map: every function and which calls which.
Explore
Run it and step through what each part did, with the code, values and plots beside it.

The data are 37 measurements from a NIST study of electron mobility in a semiconductor, as a function of the natural log of the carrier density (the Thurber data set of NIST's Statistical Reference Datasets). NIST rates the fit "higher difficulty" and publishes certified parameter values to 11 digits. The model is a ratio of two cubics in x with seven parameters,

y (x) = \frac{b _{1} + b _{2} x + b _{3} x ^{2} + b _{4} x ^{3}}{1 + b _{5} x + b _{6} x ^{2} + b _{7} x ^{3}},

and the fit minimizes the sum of squares of the 37 residuals $f_{i} = y (x_{i}) - y_{i}$ .

The fit is done by lmder from lib/minpack: a port of MINPACK's Levenberg-Marquardt code (Moré, Garbow and Hillstrom, Argonne, 1980), the same code scipy.optimize.leastsq calls. Each pass computes the Jacobian $J_{ik} = \partial f_{i} / \partial b_{k}$ , factors it, and tries steps $p$ inside a trust region $∥ D p ∥ \leq Δ$ , a radius around the current parameters where the linear model $f (b + p) \approx f (b) + J p$ of the residuals is trusted. A step that reduces the sum of squares about as much as the model predicted grows the radius; a step that does worse shrinks it and is retried. The library's intent files (lib/minpack/*.intent) explain each step, and viv story lib/minpack/story_osborne1.viv prints one run iteration by iteration.

Three references check the result:

the original Fortran, built with gfortran and instrumented to record every trial step (lib/minpack/oracle); the port does the same operations in the same order, and on this problem, which uses only arithmetic, every recorded number agrees to the last bit;
scipy.optimize.leastsq, with the same settings;
NIST's certified values.

Printed

info: 1

LM_TOL: 1.49012e-08

ntrial: 32

counts: [33, 27]

b: [1288.14, 1491.07, 583.233, 75.4155, 0.966289, 0.39797, 0.0497256]

5642.71

The fit

info 1 means lmder stopped because both the actual and the predicted relative reduction of the sum of squares fell below the tolerance $ε \approx 1.5 \times 1 0^{- 8}$ , where $ε$ is the machine precision. It took 33 evaluations of the residuals and 27 of the Jacobian; each of the 32 trial steps after the first evaluation is one row of hist.

NIST's measurementsthe fitted rational model

How it converged

The first two trial steps, starting with the Gauss-Newton step, make the fit worse and are rejected, and each rejection shrinks the trust radius. Three more are rejected in the next few passes, 5 of the 32 in all. After that every step is accepted, and the norm of the residuals approaches its final value roughly geometrically, with the trust radius no longer the limit: the plot shows, on a log scale, how far the norm still is from its final value.

Convergence, one point per trial step

How far the norm of the residuals is from its final value

norm of residuals - final norm

The trust radius

Tests 8 of 8 pass

the Fortran stops the same way: info 1, 33 evaluations of f and 27 of J (passes)
NIST's data has 37 points, and the tolerance sqrt(eps) is 1.5e-8 (passes)
the norm after each of the 32 trial steps equals the Fortran's to the last bit (passes)
the parameters equal the Fortran's to the last bit (passes)
scipy.optimize.leastsq stops with the same info and counts, at parameters within 1e-13 (passes)
every parameter is within 0.001 standard deviations of NIST's certified value (passes)
the residual sum of squares is NIST's 5642.7082397 to 1e-8 (passes)
5 of the 32 trial steps are rejected, and the norm never rises (passes)