NIST's certified linear regressions, digit by digit

NIST's Statistical Reference Datasets include eleven linear least-squares problems with coefficients, standard errors, residual standard deviation and R^2 certified to 15 digits, from the trivial (Norris, a line) to the notorious (Filip, a degree-10 polynomial; Wampler4-5, degree 5 with large residuals). vp.stats.ols fits each with every solver of its slot, and the table gives the certified digits of the worst coefficient (NIST's log relative error, LRE) beside statsmodels 0.14 and R 4.3.3's lm.

One fit, checked now

Longley's macroeconomic data: six correlated predictors and an intercept.

import json
import sys

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

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

y, X, cert = N.read("Longley")
r = vp.stats.ols(y, X, names=["const"] + ["x%d" % k for k in range(1, 7)])
print(r.summary())
print(
    "certified digits: coefficients %.1f, standard errors %.1f, R^2 %.1f"
    % (
        N.lre(r.coef, cert["coef"]).min(),
        N.lre(r.bse, cert["se"]).min(),
        N.lre(r.rsquared, cert["r2"]),
    )
)

OLS regression results (vp.stats.ols, solver 'qr')
==============================================================================
No. Observations:       16     R-squared:          0.9955     F-statistic:       330.3
Df Residuals:            9     Adj. R-squared:     0.9925     Prob (F):       4.98e-10
Df Model:                6     Log-Likelihood:     -109.6     AIC:               233.2
Residual SE:         304.9     Rank:                    7     BIC:               238.6
==============================================================================
                  coef    std err         t     P>|t|      [0.025      0.975]
-----------------------------------------------------------------------------
const       -3.482e+06  8.904e+05    -3.911   0.00356  -5.497e+06  -1.468e+06 **
x1               15.06      84.91     0.177     0.863        -177       207.2 
x2            -0.03582    0.03349    -1.070     0.313     -0.1116     0.03994 
x3               -2.02     0.4884    -4.136   0.00254      -3.125     -0.9154 **
x4              -1.033     0.2143    -4.822  0.000944      -1.518     -0.5485 ***
x5             -0.0511     0.2261    -0.226     0.826     -0.5625      0.4603 
x6                1829      455.5     4.016   0.00304       798.8        2860 **
certified digits: coefficients 13.2, standard errors 13.0, R^2 15.0

All eleven datasets

Measured by python3 lib/stats/nist_lls.py (lib/stats/nist_lls.json). "qr" is the default: pivoted Householder QR on unit-norm columns and two steps of Bjorck's refinement of the augmented system, with both of its residuals summed exactly. "attainable" is the exact least-squares solution of the design matrix as stored in doubles (mpmath, 60 digits): for Filip, x^10 rounded to double has already lost the digits past 7.6, and the QR reaches that solution to 14.4 digits.

table = json.load(open("lib/stats/nist_lls.json"))["datasets"]
cols = [
    "vipyr qr",
    "vipyr qr, no refinement",
    "vipyr cholesky",
    "vipyr svd",
    "statsmodels",
    "R lm",
    "attainable",
]
short = ["qr", "qr-1", "chol", "svd", "statsm", "R lm", "attain"]
print("%-9s" % "" + "".join("%8s" % s for s in short))
for name, row in table.items():
    print("%-9s" % name + "".join("%8.1f" % row[c]["coef"] for c in cols))

               qr    qr-1    chol     svd  statsm    R lm  attain
Norris       14.4    12.0    12.2    12.7    13.0    12.7    14.1
Pontius      13.4    12.4    11.3    12.3     6.3    12.2    13.5
NoInt1       14.7    14.7    14.7    14.7    14.7    14.7    14.7
NoInt2       15.0    15.0    15.0    15.0    15.0    15.0    15.0
Filip         7.6     7.0     0.0     0.0     0.0     0.0     7.6
Longley      13.2    10.9     7.1    11.1    10.8    13.0    14.6
Wampler1     11.6     9.5     7.0    10.1     8.9     9.4    15.0
Wampler2     13.0    12.6     9.5    12.6    10.3    13.0    13.2
Wampler3     11.2     9.0     7.0     9.9     8.9     9.4    15.0
Wampler4     10.6     7.6     6.3     7.9     7.9     7.8    15.0
Wampler5     15.0     5.6     6.0     5.9     5.9     5.8    15.0

The picture

Certified digits of the worst coefficient, per dataset.

names = list(table)
fig = figure("Certified digits of the worst coefficient", "dataset", "LRE")
for k, (c, lab) in enumerate(
    [
        ("vipyr qr", "vp.stats qr"),
        ("statsmodels", "statsmodels"),
        ("R lm", "R lm"),
    ]
):
    fig.bars(
        list(range(len(names))),
        [table[n][c]["coef"] for n in names],
        label=lab,
    )
for i, n in enumerate(names):
    fig.xtick(n, i)
fig

Source

gallery/stats/plot_nist_linear_regression.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 lib/stats/gallery_longley.viv: its atlas, the defs and the calls between them, or open it in the workbench.