Master Bayesian Linear Regression with PyMC3: A Hands‑On Guide

Probability Programming Basics – Linear Regression

Probabilistic programming enables automatic Bayesian inference on user‑defined probability models. Modern MCMC methods such as Hamiltonian Monte Carlo (HMC) and its self‑tuning variant No‑U‑Turn Sampler (NUTS) require gradient information, which may not always be available. PyMC3 is an open‑source Python framework that leverages Theano for gradient computation and provides C‑accelerated operations, allowing models to be defined directly in Python.

PyMC3 includes advanced samplers like NUTS and HMC that perform well on high‑dimensional posterior distributions, offering faster convergence than traditional methods. These samplers automatically obtain gradients from the likelihood, and NUTS adapts step sizes without user‑specified tuning.

Example

Consider a simple Bayesian linear regression with normal priors on the intercept, coefficients, and observation noise. The model predicts observations as a linear combination of two predictors plus Gaussian noise.

Simulating Observation Data

We generate synthetic data using NumPy and then attempt to recover the parameters with PyMC3.

<code>import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

np.random.seed(123)
alpha = 1
sigma = 1
beta = [1, 2.5]
N = 100
X1 = np.random.randn(N)
X2 = np.random.randn(N)
Y = alpha + beta[0]*X1 + beta[1]*X2 + np.random.randn(N)*sigma

fig1, ax1 = plt.subplots(1, 2, figsize=(10,4))
ax1[0].scatter(X1, Y); ax1[0].set_xlabel('X1'); ax1[0].set_ylabel('Y')
ax1[1].scatter(X2, Y); ax1[1].set_xlabel('X2'); ax1[1].set_ylabel('Y')

fig2 = plt.figure(2)
ax2 = Axes3D(fig2)
ax2.scatter(X1, X2, Y)
ax2.set_xlabel('X1'); ax2.set_ylabel('X2'); ax2.set_zlabel('Y')
</code>

Defining the Model

PyMC3’s model definition syntax mirrors statistical formulas, with each statement representing a component of the model.

<code>import pymc3 as pm

basic_model = pm.Model()
with basic_model:
    alpha = pm.Normal('alpha', mu=0, sd=10)
    beta = pm.Normal('beta', mu=0, sd=10, shape=2)
    sigma = pm.HalfNormal('sigma', sd=1)
    mu = alpha + beta[0]*X1 + beta[1]*X2
    Y_obs = pm.Normal('Y_obs', mu=mu, sd=sigma, observed=Y)
</code>

The with block creates a context where all random variables are added to basic_model . The observed variable Y_obs links the data to the model.

Model Fitting

Two approaches can be used to estimate the posterior: (1) MAP estimation via optimization, and (2) full MCMC sampling.

Maximum A Posteriori (MAP)

PyMC3’s find_MAP function performs optimization (default BFGS) to locate the posterior mode.

<code>map_estimate = pm.find_MAP(model=basic_model)
</code>

MAP may fail for multimodal or high‑dimensional posteriors.

Sampling Algorithms

For full posterior inference, PyMC3 offers several step methods: NUTS, Metropolis, Slice, HamiltonianMC, and BinaryMetropolis. The library can auto‑select an appropriate sampler based on variable types, but users may manually specify one.

Binary variables → BinaryMetropolis

Discrete variables → Metropolis

Continuous variables → NUTS

Gradient‑Based Sampling

NUTS leverages gradient information from the log‑probability to achieve efficient exploration of continuous parameter spaces. PyMC3 uses Theano’s automatic differentiation for these gradients and can initialize NUTS with ADVI when no step method is provided.

The following code runs 2000 NUTS iterations on the model.

<code>with basic_model:
    trace = pm.sample(2000)
</code>

The resulting trace object stores sampled values for each variable.

<code>print("alpha", trace['alpha'].shape)
print(trace['alpha'][0:5], "\n")
print("beta", trace['beta'].shape)
print(trace['beta'], "\n")
print("sigma", trace['sigma'].shape)
print(trace['sigma'])
</code>

Specific samplers can be assigned manually; for example, using a Slice sampler for sigma while letting PyMC3 choose NUTS for the remaining variables.

<code>with basic_model:
    # Obtain MAP for initialization
    start = pm.find_MAP(fmin=optimize.fmin_powell)
    # Instantiate Slice sampler for sigma
    step = pm.Slice(vars=[sigma])
    # Sample 5000 iterations
    trace = pm.sample(5000, step=step, start=start)
</code>

Posterior Analysis

PyMC3 provides traceplot to visualize marginal histograms and trace lines for each sampled variable, as well as summary for numerical statistics.

Reference

https://blog.csdn.net/jackxu8/article/details/71080865