Introduction

Bayesian modeling provides a flexible way to incorporate prior beliefs and quantify uncertainty in your data analysis. If you’re new to Bayesian thinking, a simple linear regression model is often the best place to start.

In this article, we’ll walk through your first Bayesian model, covering prior specification, Markov Chain Monte Carlo (MCMC) sampling, and essential diagnostic plots using ArviZ. We will also utilize PyMC, a probabilistic programming library allowing us to build Bayesian models and fit them using Markov Chain Monte Carlo (MCMC) techniques.

By the end, you’ll have a concise overview of how to build, fit, and check a Bayesian linear regression model in Python.

Building a Simple Bayesian Linear Regression

The foundation of Bayesian modeling is combining prior beliefs with observed data to obtain a posterior distribution. Here, our linear regression setup is:

y = α + βX + ϵ

where α (intercept) and β (slope) are unknown parameters, and ϵ is normally distributed noise. We assume ϵ ∼ N (0, σ²).

Prior Specification

In a Bayesian framework, we place priors on α, β, and σ. If you are unsure about your parameters’ likely values, you can start with weakly informative or uninformative priors. In this example:

• α: We place a Normal(0, 10) prior, meaning we expect the intercept to be around 0 but allow wide variation
• β: We also place a Normal(0, 10) prior, reflecting minimal knowledge about the slope’s magnitude
• σ: We use a HalfNormal(1) prior, which restricts σ to be positive and typically near 0–2 but allows for uncertainty

Here’s the basic Python code to specify these priors and the model in PyMC.

import pymc as pm
import arviz as az

# Suppose X and y are your data arrays
X = ...
y = ...

with pm.Model() as linear_model:
    # Priors
    alpha = pm.Normal('alpha', mu=0, sigma=10)
    beta = pm.Normal('beta', mu=0, sigma=10)
    sigma = pm.HalfNormal('sigma', sigma=1)
    
    # Linear model
    mu = alpha + beta * X
    
    # Likelihood
    y_obs = pm.Normal('y_obs', mu=mu, sigma=sigma, observed=y)

MCMC Sampling Basics

Markov Chain Monte Carlo (MCMC) is a family of algorithms for sampling from complex probability distributions. In practice, we often use Hamiltonian Monte Carlo (HMC) or the No-U-Turn Sampler (NUTS) to efficiently explore the posterior. The idea is to create a chain of samples that converge to the target distribution, yielding estimates for the parameters.

After defining the model, we sample from the posterior distribution. Here’s how you might do it in code:

with linear_model:
    trace = pm.sample(draws=2000, tune=1000, target_accept=0.95, random_seed=42)

Key points in the sampling step:

• draws: How many posterior samples to keep after tuning.
• tune: Extra iterations to allow the sampler to “warm up,” which you usually discard.
• target_accept: Adjusts the acceptance rate in HMC/NUTS; a value around 0.80–0.95 often works well.
• random_seed: Ensures reproducibility.

You’ll end up with a trace object containing posterior samples of α, β, and σ. These samples represent plausible values of the parameters given the data and priors.

Assessing Model Performance with ArviZ

ArviZ is a companion library designed for analyzing MCMC outputs. It provides plots like the trace plot, rank plot, and posterior distribution plots to help ensure that your chains converge and mix well. Here’s a typical workflow:

az.plot_trace(trace)
az.plot_posterior(trace)
az.summary(trace)

Explanation of the above workflow code:

• Trace Plots: Show the parameter values over sampling iterations. Look for “fuzzy caterpillars” with no major trends.
• Posterior Plots: Visualize the posterior distributions of your parameters.
• ArviZ Summary: Provides numerical summaries (mean, standard deviation, highest posterior density intervals) of each parameter.

A Full Example

import numpy as np
import matplotlib.pyplot as plt
import pymc as pm
import arviz as az

# 1. Generate synthetic data for a simple linear regression
np.random.seed(42)
alpha_true = 2.5
beta_true = 1.2
sigma_true = 0.5

# Generate X and corresponding y
X = np.random.randn(100)
y = alpha_true + beta_true * X + np.random.normal(0, sigma_true, size=len(X))

# 2. Specify and build the Bayesian model using PyMC
with pm.Model() as linear_model:
    # Priors
    alpha = pm.Normal('alpha', mu=0, sigma=10)
    beta = pm.Normal('beta', mu=0, sigma=10)
    sigma = pm.HalfNormal('sigma', sigma=1)
    
    # Linear model
    mu = alpha + beta * X
    
    # Likelihood
    y_obs = pm.Normal('y_obs', mu=mu, sigma=sigma, observed=y)
    
    # 3. Run MCMC sampling
    trace = pm.sample(
        draws=2000,      # total posterior samples to collect
        tune=1000,       # "warm-up" iterations discarded after tuning
        target_accept=0.95,
        random_seed=42
    )

# 4. Use ArviZ to visualize and summarize the posterior, save visualizations
az.plot_trace(trace, var_names=['alpha', 'beta', 'sigma'])
plt.savefig("trace_plot.png")
plt.clf()  # Clear the current figure

az.plot_posterior(trace, var_names=['alpha', 'beta', 'sigma'])
plt.savefig("posterior_plot.png")
plt.clf()  # Clear the current figure

summary_df = az.summary(trace, var_names=['alpha', 'beta', 'sigma'])
print(summary_df)

Output will be something similar to the following:

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [alpha, beta, sigma]
Sampling 4 chains, 0 divergences ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 100% 0:00:00 / 0:00:00
Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 1 seconds.
        mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_bulk  ess_tail  r_hat
alpha  2.503  0.049   2.414    2.597      0.001      0.0    7295.0    5310.0    1.0
beta   1.129  0.053   1.037    1.237      0.001      0.0    7626.0    5377.0    1.0
sigma  0.480  0.034   0.416    0.542      0.000      0.0    7285.0    5839.0    1.0

And the trace and posterior plots, respectively:

Figure 1: Trace plot

Figure 2: Posterior plot

How it works:

1. Synthetic Data: We generate random 𝑋 values from a standard normal distribution. The true model is 𝑦 = 𝛼_true + 𝛽_true × 𝑋 plus normally distributed noise.
2. Model Definition: We define priors on α, β, and σ, and specify a normal likelihood for y.
3. MCMC Sampling: PyMC’s sample() function uses modern algorithms (like NUTS) to draw samples from the posterior distribution.
4. Diagnostics & Summary: ArviZ provides plots and numerical summaries to help you check convergence and interpret the posterior distributions.

Run this script from start to finish, and you’ll have a complete demonstration of the Bayesian workflow — from data generation to final model insights.

Conclusion

By combining prior information with observed data, Bayesian linear regression allows you to quantify uncertainty in your estimates and model parameters more explicitly than traditional approaches. You start by specifying priors, run MCMC to explore the posterior, and diagnose results using ArviZ. This simple example is a starting point, and once you’re comfortable with these basics, you can extend Bayesian techniques to more complex scenarios such as hierarchical models, time series, or generalized linear models.