Simplest Bayesian regressions

Source: Simplest/simplest.Rmd

The Bayesian companion uses lapylace directly: a Gaussian likelihood, weakly informative priors, and posterior predictive draws.

Simulated regression

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(18901)
fake = pd.DataFrame({"x": rng.binomial(1, .5, size=120)})
fake["y"] = 1 + 2.5*fake.x + rng.normal(0, 1.2, size=len(fake))
fit = fit_glm("y ~ x", fake, seed=18902, prior_scale=5, intercept_scale=5, aux_scale=2)
fit.summary(["alpha", "beta", "sigma"])

	Mean	MCSE	StdDev	5%	50%	95%	N_Eff	N_Eff/s	R_hat
alpha	1.09432	0.005871	0.145242	0.845253	1.09222	1.32692	611.9790	11768.80000	1.00236
beta[1]	2.46720	0.008657	0.212531	2.126320	2.46324	2.80889	602.6660	11589.70000	1.00228
sigma	1.16629	0.003110	0.072350	1.057620	1.16385	1.29013	542.2801	10428.46343	1.00429

new = pd.DataFrame({"x": [0, 1]})
pred = fit.posterior_predict(new, rng=np.random.default_rng(18903))
pd.DataFrame({"x": new.x, "pred_mean": pred.mean(axis=0), "lo80": np.quantile(pred,.1,axis=0), "hi80": np.quantile(pred,.9,axis=0)})

	x	pred_mean	lo80	hi80
0	0	1.105348	-0.378289	2.571845
1	1	3.540268	1.999323	5.041428