Linear regression across Python tools

Source: DifferentSoftware/linear.Rmd

The comparison page now centers on lapylace for the Bayesian formula interface and NumPy for the closed-form least-squares check.

Simulated data

from pathlib import Path
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import lapylace as lp

root = Path("../../ROS-Examples")

def coef_medians(fit):
    return pd.Series(
        [np.median(fit.alpha_draws()), *np.median(fit.beta_draws(), axis=0)],
        index=["Intercept", *fit.columns],
    )

def fit_glm(formula, data, family=None, seed=1, prior_scale=2.5, intercept_scale=5, aux_scale=10, **kwargs):
    return lp.stan_glm(
        formula,
        data=data,
        family=family or lp.gaussian(),
        prior=lp.normal(0, prior_scale),
        prior_intercept=lp.normal(0, intercept_scale),
        prior_aux=lp.exponential(aux_scale),
        chains=2,
        parallel_chains=2,
        iter_warmup=300,
        iter_sampling=500,
        seed=seed,
        refresh=100,
        **kwargs,
    )

rng = np.random.default_rng(19001)
dat = pd.DataFrame({"x1": rng.normal(size=80), "x2": rng.normal(size=80)})
dat["y"] = 1 + 2*dat.x1 - dat.x2 + rng.normal(0, .8, size=len(dat))
X = np.column_stack([np.ones(len(dat)), dat[["x1", "x2"]]])
np_coef = np.linalg.lstsq(X, dat.y, rcond=None)[0]
fit = fit_glm("y ~ x1 + x2", dat, seed=19002, prior_scale=5, intercept_scale=5, aux_scale=2)
pd.DataFrame({"numpy_lstsq": np_coef, "lapylace_median": coef_medians(fit).to_numpy()}, index=["Intercept", "x1", "x2"]).round(3)

	numpy_lstsq	lapylace_median
Intercept	0.840	0.841
x1	1.994	1.989
x2	-0.989	-0.994