Bayesian Student-t Linear Model — PyMC (continuous ν)¶
Companion to student_t_lm_python.ipynb (the from-scratch Geweke 1993 Gibbs sampler). There, ν is drawn from a discrete grid {2,…,40}; here PyMC treats ν as a continuous parameter with a weakly-informative Gamma prior and samples (β, σ, ν) jointly by NUTS via pm.StudentT. We validate against the discrete-ν Gibbs on the Section 1 synthetic data and on three Nelson-Plosser series.
| Discrete Gibbs (from-scratch) | PyMC / NUTS (here) | |
|---|---|---|
| ν | discrete {2,…,40} | continuous $(0,\infty)$, Gamma(2, 0.1) prior |
| sampler | 4-block Gibbs | NUTS (HMC) |
| latent λᵢ | explicit (outlier weights) | marginalised out |
In [1]:
import os, sys
conda_lib = os.path.join(os.environ.get('CONDA_PREFIX',''), 'Library', 'bin')
if os.path.isdir(conda_lib) and conda_lib not in os.environ.get('PATH',''):
os.environ['PATH'] = conda_lib + ';' + os.environ.get('PATH','') # Windows DLL fix (PyMC needs this on PATH)
import numpy as np, pandas as pd, matplotlib.pyplot as plt
import pymc as pm, arviz as az
from scipy.stats import t as t_dist, gaussian_kde
sys.path.insert(0, r'c:\Users\user\project')
from student_t_lm_gibbs import student_t_gibbs # the from-scratch Gibbs, for side-by-side overlay
print('Import OK')
g++ not available, if using conda: `conda install gxx`
Import OK
1. Synthetic data: continuous-ν NUTS vs discrete-ν Gibbs¶
In [2]:
# §1 synthetic data (identical to the from-scratch notebook) + run the discrete-ν Gibbs for comparison
rng1 = np.random.default_rng(0); n, k = 120, 3
beta_true = np.array([2.0, 0.8, -0.5]); sigma2_true, nu_true = 1.0, 4.0
X1 = np.column_stack([np.ones(n), rng1.standard_normal((n, k - 1))])
y1 = X1 @ beta_true + t_dist.rvs(df=nu_true, scale=np.sqrt(sigma2_true), size=n, random_state=rng1)
out_t = student_t_gibbs(y1, X1, R=15000, burn=3000, seed=1, nprint=0)
names1 = ['β[0]', 'β[1]', 'β[2]']
print('discrete-ν Gibbs: ν mean = %.1f' % out_t['nu'].mean())
Done in 0.0 min | kept draws: 12000 discrete-ν Gibbs: ν mean = 8.6
In [3]:
import pymc as pm
import arviz as az
# Use the same §1 data: y1, X1, beta_true, nu_true
with pm.Model() as model_t:
# Priors
beta_pm = pm.Flat('beta', shape=k) # flat prior on β
sigma = pm.HalfNormal('sigma', sigma=2.0) # half-normal on σ (>0)
nu_pm = pm.Gamma('nu', alpha=2.0, beta=0.1) # weakly informative; mean=20
mu = pm.math.dot(X1, beta_pm)
# Likelihood: Student-t with continuous ν
obs = pm.StudentT('obs', nu=nu_pm, mu=mu, sigma=sigma, observed=y1)
# Sample
idata = pm.sample(
draws=4000, tune=2000, chains=2,
target_accept=0.9,
random_seed=42, progressbar=False,
)
print(az.summary(idata, var_names=['beta', 'sigma', 'nu']))
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [beta, sigma, nu]
Sampling 2 chains for 2_000 tune and 4_000 draw iterations (4_000 + 8_000 draws total) took 4 seconds.
We recommend running at least 4 chains for robust computation of convergence diagnostics
mean sd eti89_lb eti89_ub ess_bulk ess_tail r_hat mcse_mean \
beta[0] 2.031 0.109 1.9 2.2 8991 6147 1.00 0.0012
beta[1] 0.791 0.118 0.6 0.98 8638 5822 1.00 0.0013
beta[2] -0.664 0.107 -0.83 -0.49 7637 5495 1.00 0.0012
sigma 1.062 0.126 0.87 1.3 5059 5468 1.00 0.0018
nu 8.9 7 3.1 22 5132 5356 1.00 0.099
mcse_sd
beta[0] 0.00081
beta[1] 0.00094
beta[2] 0.00089
sigma 0.0012
nu 0.19
In [4]:
# ── Side-by-side comparison: Gibbs (discrete) vs PyMC (continuous) ────────────
nu_gibbs = out_t['nu']
beta_gibbs = out_t['beta']
# PyMC posterior draws (flatten chains)
nu_pymc = idata.posterior['nu'].values.ravel()
beta_pymc = idata.posterior['beta'].values.reshape(-1, k)
sigma_pymc = idata.posterior['sigma'].values.ravel()
# β comparison table
names1 = ['β[0]', 'β[1]', 'β[2]']
print('β posteriors (true = [2.0, 0.8, -0.5])\n')
print(f'{"":6s} {"Gibbs mean":>12s} {"Gibbs sd":>10s} {"PyMC mean":>12s} {"PyMC sd":>10s}')
for j, nm in enumerate(names1):
gm, gs = beta_gibbs[:, j].mean(), beta_gibbs[:, j].std()
pm_, ps = beta_pymc[:, j].mean(), beta_pymc[:, j].std()
print(f'{nm:6s} {gm:12.3f} {gs:10.3f} {pm_:12.3f} {ps:10.3f}')
print(f'\nσ² posterior:')
print(f' Gibbs σ² mean = {out_t["sigma2"].mean():.3f} sd = {out_t["sigma2"].std():.3f}')
print(f' PyMC σ² mean = {(sigma_pymc**2).mean():.3f} sd = {(sigma_pymc**2).std():.3f} (true σ²=1.0)')
print(f'\nν posterior:')
print(f' Gibbs (discrete) mean={nu_gibbs.mean():.1f} median={np.median(nu_gibbs):.0f}'
f' mode={float(pd.Series(nu_gibbs).mode().iloc[0]):.0f}')
print(f' PyMC (continuous) mean={nu_pymc.mean():.1f} median={np.median(nu_pymc):.1f}'
f' q025={np.quantile(nu_pymc, 0.025):.1f} q975={np.quantile(nu_pymc, 0.975):.1f}')
print(f' True ν = {nu_true:.0f}')
β posteriors (true = [2.0, 0.8, -0.5])
Gibbs mean Gibbs sd PyMC mean PyMC sd
β[0] 2.028 0.111 2.031 0.109
β[1] 0.794 0.116 0.791 0.118
β[2] -0.663 0.108 -0.664 0.107
σ² posterior:
Gibbs σ² mean = 1.102 sd = 0.288
PyMC σ² mean = 1.143 sd = 0.270 (true σ²=1.0)
ν posterior:
Gibbs (discrete) mean=8.6 median=6 mode=4
PyMC (continuous) mean=8.9 median=6.6 q025=2.7 q975=28.4
True ν = 4
In [5]:
fig, axes = plt.subplots(1, 4, figsize=(16, 4))
# β[0], β[1], β[2] — overlaid posteriors
for j, ax in enumerate(axes[:3]):
lo = min(np.percentile(beta_gibbs[:, j], 0.5),
np.percentile(beta_pymc[:, j], 0.5))
hi = max(np.percentile(beta_gibbs[:, j], 99.5),
np.percentile(beta_pymc[:, j], 99.5))
bins = np.linspace(lo, hi, 60)
ax.hist(beta_gibbs[:, j], bins=bins, density=True,
alpha=0.55, color='orange', label='Gibbs (discrete ν)')
ax.hist(beta_pymc[:, j], bins=bins, density=True,
alpha=0.45, color='steelblue', label='PyMC (continuous ν)', histtype='step', lw=2)
ax.axvline(beta_true[j], color='k', ls='--', lw=1.5, label=f'True {beta_true[j]}')
ax.set_title(names1[j]); ax.legend(fontsize=7)
# ν posteriors — discrete bar vs continuous KDE
ax = axes[3]
nu_vals, nu_cnt = np.unique(nu_gibbs, return_counts=True)
ax.bar(nu_vals, nu_cnt / len(nu_gibbs), width=0.7, color='orange',
alpha=0.7, label='Gibbs (discrete)', zorder=2)
# KDE for PyMC continuous ν
from scipy.stats import gaussian_kde
kde = gaussian_kde(nu_pymc, bw_method=0.3)
nu_grid_plot = np.linspace(1, 40, 300)
ax.plot(nu_grid_plot, kde(nu_grid_plot), color='steelblue', lw=2,
label='PyMC (continuous KDE)')
ax.axvline(nu_true, color='k', ls='--', lw=1.5, label=f'True ν={nu_true:.0f}')
ax.set_xlabel('ν'); ax.set_title('ν posterior: Gibbs vs PyMC')
ax.legend(fontsize=7); ax.set_xlim(1, 30)
fig.suptitle('§3 Comparison — Geweke Gibbs vs PyMC NUTS on §1 synthetic data\n'
'β posteriors agree; ν — discrete Gibbs vs continuous KDE', fontsize=10)
plt.tight_layout(); plt.show()
The first three panels overlay the coefficient posteriors from the discrete-ν Gibbs (orange fill) and the continuous-ν PyMC/NUTS (blue step) — they coincide to well within a posterior SD, so the β inference doesn't care how ν is treated. The fourth panel contrasts the two ν posteriors: the Gibbs a discrete bar distribution (mode at the true 4), PyMC a continuous density (median ≈ 6.6) — different shapes, the same heavy-tailed verdict.
2. Nelson-Plosser: ν posterior, Gibbs vs PyMC (three representative series)¶
In [6]:
# Nelson-Plosser load (focus series) + discrete-ν Gibbs modes for the comparison titles
np_df = pd.read_csv(r'nelplo.csv')
log_series = ['cpi','ip','gnp.nom','vel','emp','nom.wages','gnp.def','money.stock','gnp.real','stock.prices','gnp.capita','real.wages']
series_labels = {'cpi':'CPI','ip':'Industrial production','gnp.nom':'Nominal GNP','vel':'Velocity','emp':'Employment',
'int.rate':'Interest rate','nom.wages':'Nominal wages','gnp.def':'GNP deflator','money.stock':'Money stock',
'gnp.real':'Real GNP','stock.prices':'Stock prices','gnp.capita':'GNP per capita','real.wages':'Real wages','unemp':'Unemployment rate'}
np_df = np_df[(np_df['year']>=1909)&(np_df['year']<=1970)].copy()
for s in log_series: np_df[s] = np.log(np_df[s])
focus_series = ['gnp.real', 'int.rate', 'cpi']
results_np = {}
for s in focus_series:
sub = np_df[['year', s]].dropna(); y_s = sub[s].values
t_s = (sub['year'].values - sub['year'].mean()) / sub['year'].std(); X_s = np.column_stack([np.ones(len(y_s)), t_s])
o = student_t_gibbs(y_s, X_s, R=12000, burn=3000, seed=hash(s) % (2**31), nprint=0)
results_np[s] = {'nu_mode': float(pd.Series(o['nu']).mode().iloc[0])}
print('focus series Gibbs ν modes:', {s: results_np[s]['nu_mode'] for s in focus_series})
Done in 0.0 min | kept draws: 9000
Done in 0.0 min | kept draws: 9000
Done in 0.0 min | kept draws: 9000
focus series Gibbs ν modes: {'gnp.real': 2.0, 'int.rate': 6.0, 'cpi': 19.0}
In [7]:
focus_series = ['gnp.real', 'int.rate', 'cpi']
pymc_np_results = {}
for s in focus_series:
sub = np_df[['year', s]].dropna()
y_s = sub[s].values
t_s = (sub['year'].values - sub['year'].mean()) / sub['year'].std()
X_s = np.column_stack([np.ones(len(y_s)), t_s])
n_s = len(y_s)
with pm.Model() as m:
beta_p = pm.Flat('beta', shape=2)
sigma_p = pm.HalfNormal('sigma', sigma=2.0)
nu_p = pm.Gamma('nu', alpha=2.0, beta=0.1)
mu_p = pm.math.dot(X_s, beta_p)
pm.StudentT('obs', nu=nu_p, mu=mu_p, sigma=sigma_p, observed=y_s)
idata_s = pm.sample(4000, tune=2000, chains=2, target_accept=0.9,
random_seed=hash(s) % (2**31),
progressbar=False)
nu_draws = idata_s.posterior['nu'].values.ravel()
pymc_np_results[s] = {
'nu_mean': nu_draws.mean(),
'nu_median': np.median(nu_draws),
'nu_q025': np.quantile(nu_draws, 0.025),
'nu_q975': np.quantile(nu_draws, 0.975),
'nu_draws': nu_draws,
}
r = pymc_np_results[s]
print(f'{series_labels[s]:20s} '
f'PyMC ν mean={r["nu_mean"]:.1f} median={r["nu_median"]:.1f}'
f' 95% CrI [{r["nu_q025"]:.1f}, {r["nu_q975"]:.1f}] | '
f'Gibbs mode={results_np[s]["nu_mode"]:.0f}')
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [beta, sigma, nu]
Sampling 2 chains for 2_000 tune and 4_000 draw iterations (4_000 + 8_000 draws total) took 4 seconds.
We recommend running at least 4 chains for robust computation of convergence diagnostics
Initializing NUTS using jitter+adapt_diag...
Real GNP PyMC ν mean=3.7 median=2.4 95% CrI [1.2, 16.4] | Gibbs mode=2
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [beta, sigma, nu]
Sampling 2 chains for 2_000 tune and 4_000 draw iterations (4_000 + 8_000 draws total) took 4 seconds.
We recommend running at least 4 chains for robust computation of convergence diagnostics
Initializing NUTS using jitter+adapt_diag...
Interest rate PyMC ν mean=16.1 median=13.0 95% CrI [3.7, 46.6] | Gibbs mode=6
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [beta, sigma, nu]
Sampling 2 chains for 2_000 tune and 4_000 draw iterations (4_000 + 8_000 draws total) took 4 seconds.
We recommend running at least 4 chains for robust computation of convergence diagnostics
CPI PyMC ν mean=24.0 median=20.9 95% CrI [5.8, 58.7] | Gibbs mode=19
In [8]:
fig, axes = plt.subplots(1, 3, figsize=(15, 4))
for ax, s in zip(axes, focus_series):
# Gibbs discrete bars
sub = np_df[['year', s]].dropna()
y_s = sub[s].values
t_s = (sub['year'].values - sub['year'].mean()) / sub['year'].std()
X_s = np.column_stack([np.ones(len(y_s)), t_s])
out_g = student_t_gibbs(y_s, X_s, R=12000, burn=3000,
seed=hash(s) % (2**31), nprint=0)
nu_g = out_g['nu']
nu_vals, nu_cnt = np.unique(nu_g, return_counts=True)
ax.bar(nu_vals, nu_cnt / len(nu_g), width=0.7,
color='orange', alpha=0.75, label='Gibbs (discrete)', zorder=2)
# PyMC continuous KDE
nu_p = pymc_np_results[s]['nu_draws']
kde = gaussian_kde(nu_p, bw_method=0.25)
xg = np.linspace(1, 42, 400)
ax.plot(xg, kde(xg), color='steelblue', lw=2,
label='PyMC (continuous KDE)')
ax.set_xlim(1, 42); ax.set_xticks([5, 15, 25, 35])
ax.set_xlabel('ν')
ax.set_title(f'{series_labels[s]}\n'
f'Gibbs mode={results_np[s]["nu_mode"]:.0f} '
f'PyMC median={pymc_np_results[s]["nu_median"]:.1f}')
ax.legend(fontsize=8)
fig.suptitle('§3b Nelson-Plosser — ν posterior: Gibbs (orange bars) vs PyMC NUTS (blue KDE)\n'
'Three representative series spanning heavy-tailed → near-Normal', fontsize=11)
plt.tight_layout(); plt.show()
Done in 0.0 min | kept draws: 9000
Done in 0.0 min | kept draws: 9000
Done in 0.0 min | kept draws: 9000
Three Nelson–Plosser series spanning the range — the Gibbs ν posterior (orange bars) against the PyMC continuous ν (blue KDE). Both pin Real GNP at low ν (strongly heavy-tailed: Gibbs mode 2, PyMC median 2.5). For the near-Normal CPI the two diverge in the expected way — the continuous ν is prior-sensitive (the Gamma(2, 0.1), mean-20 prior pulls it up) while the discrete grid caps at 40 — the same story, different behaviour at the Normal boundary.
Results¶
- β posteriors agree between the discrete-ν Gibbs and continuous-ν PyMC to well within a posterior SD — the inference on the regression coefficients is robust to how ν is treated.
- ν differs in shape, not story: the Gibbs gives a discrete distribution (mode near the truth), PyMC a continuous one; both identify heavy tails. For borderline/near-Normal series the continuous ν is prior-sensitive (the Gamma(2,0.1) mean-20 prior pulls toward Normal), while the discrete grid caps at 40.
- Trade-off: PyMC is more flexible (continuous ν, easy model extensions) but marginalises out the λᵢ, so it loses the per-observation outlier weights the Gibbs exposes directly. For robust diagnostics, the scale-mixture Gibbs is preferable; for flexible modelling, PyMC.