`vp.stats.ols`

Least squares with statsmodels' inference table, by pivoted QR, checked against NIST's certified results.

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

vp.stats.ols(y, X) fits y = X b + e by least squares and returns what statsmodels' OLS(y, X).fit() reports: coefficients, standard errors, t and p-values, 95% intervals, R^2 and adjusted R^2, the F test, the residual standard deviation, the log-likelihood, AIC and BIC. result.summary() prints the table. weights= gives weighted least squares (statsmodels' WLS).

The solver is a named slot (solver=): "qr", the default, is Householder QR with column pivoting on X with its columns scaled to unit norm, followed by one step of iterative refinement in which the residual y - X b is summed exactly (each product split by fma, the terms added by fsum); a column whose pivot falls below tol times the first is dropped, its coefficient NaN, as R's lm reports NA. "cholesky" solves the normal equations and "svd" takes the pseudo-inverse (statsmodels' default). They are there to compare: on NIST's Filip and Wampler problems the normal equations lose the digits the QR keeps (lib/stats/nist_lls.json).

A constant column is the intercept; add one yourself (no formula interface). R^2 is centered when there is one, as in statsmodels and R. AIC and BIC follow statsmodels (R's AIC counts the residual variance as a parameter and is larger by 2).

Signature

vp.stats.ols(y, X, *, weights=None, solver='qr', refine=True, tol=1e-10, names=None)

Parameters

y : ndarray[float64], required: The response, n values.
X : ndarray[float64], required: The n by k design matrix; a constant column is the intercept (add one yourself).
weights : ndarray[float64], default None: Weights w_i: weighted least squares, minimizing sum w_i r_i^2.
solver : str, default 'qr': The least-squares slot: 'qr' (pivoted Householder QR), 'cholesky' (normal equations) or 'svd' (pseudo-inverse).
refine : bool, default True: With 'qr', one step of iterative refinement with an exactly summed residual.
tol : float, default 1e-10: Rank: a pivot below tol times the first (QR), or a singular value below tol times the largest (SVD), drops its column.
names : list, default None: Names of the columns (a DataFrame's columns by default).

Returns

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

coef : ndarray[float64]: The coefficients (NaN for a dropped column).
bse : ndarray[float64]: Their standard errors.
tvalues : ndarray[float64]: Coef / bse.
pvalues : ndarray[float64]: Two-sided p-values from Student's t with df_resid.
conf_int : ndarray[float64]: 95% confidence intervals, k by 2.
rsquared : float: R^2 (centered when there is a constant column, as statsmodels).
rsquared_adj : float: Adjusted R^2.
fvalue : float: The F statistic of the model against the constant.
f_pvalue : float: Its p-value.
ssr : float: The residual sum of squares (weighted).
centered_tss : float: The total sum of squares about the (weighted) mean.
scale : float: The residual variance ssr / df_resid.
resid_sd : float: The residual standard deviation, sqrt(scale).
df_model : int: Rank minus the constant.
df_resid : int: N minus rank.
nobs : int: N.
rank : int: The rank found by the solver.
llf : float: The Gaussian log-likelihood.
aic : float: Akaike's criterion, statsmodels' convention (R's AIC counts sigma: + 2).
bic : float: The Bayesian criterion, statsmodels' convention.
resid : ndarray[float64]: Y - X coef.
fittedvalues : ndarray[float64]: X coef.
cov_params : ndarray[float64]: Scale (X^T W X)^-1.
names : list: The column names.
solver : str: The solver used.
weighted : bool: Whether weights were given.
k_constant : int: 1 when X has a constant column.

Raises

ValueError: A solver that is not registered, or y and X of different lengths.
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.

Examples

import numpy as np
import vipyr as vp
rng = np.random.default_rng(1)
x = rng.uniform(0, 10, 50)
X = np.column_stack([np.ones(50), x])
y = 2.0 + 0.5 * x + rng.normal(0, 1, 50)
r = vp.stats.ols(y, X, names=["const", "x"])
print(r.summary())
assert r.check()

OLS regression results (vp.stats.ols, solver 'qr')
==============================================================================
No. Observations:       50     R-squared:          0.7435     F-statistic:       139.2
Df Residuals:           48     Adj. R-squared:     0.7382     Prob (F):       8.64e-16
Df Model:                1     Log-Likelihood:     -61.54     AIC:               127.1
Residual SE:        0.8456     Rank:                    2     BIC:               130.9
==============================================================================
                  coef    std err         t     P>|t|      [0.025      0.975]
-----------------------------------------------------------------------------
const            1.854     0.2541     7.296  2.59e-09       1.343       2.365 ***
x               0.5089    0.04314    11.797  8.64e-16      0.4222      0.5956 ***

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

Checked against

Reference: NIST StRD's 11 linear regression datasets (Norris, Pontius, NoInt1, NoInt2, Filip, Longley, Wampler1-5): certified digits of every coefficient, standard error, the residual SD and R^2, beside statsmodels 0.14 and R 4.3.3's lm (lib/stats/nist_lls.py, lib/stats/nist_lls.json).

vp.stats.ols