Bayesian Machine Learning

Bayesian vs Frequentist Thinking

Two fundamentally different interpretations of probability:

	Frequentist	Bayesian
Probability means	Long-run frequency	Degree of belief
Parameters	Fixed unknowns	Random variables with distributions
Data	Repeated samples	Fixed, observed
Inference output	Point estimate + confidence interval	Full posterior distribution
Prior knowledge	Not used	Explicitly encoded as prior
Uncertainty	Frequentist CI (misinterpreted as credible interval)	Credible interval with correct probability interpretation

The Bayesian Update

$\underbrace{P(\theta \mid \mathcal{D})}_{\text{posterior}} = \frac{\overbrace{P(\mathcal{D} \mid \theta)}^{\text{likelihood}} \cdot \overbrace{P(\theta)}^{\text{prior}}}{\underbrace{P(\mathcal{D})}_{\text{marginal likelihood (evidence)}}}$

• Prior $P(\theta)$: what we believe about parameters before seeing data
• Likelihood $P(\mathcal{D} \mid \theta)$: how probable the observed data is given parameters
• Posterior $P(\theta \mid \mathcal{D})$: updated belief after observing data
• Evidence $P(\mathcal{D}) = \int P(\mathcal{D}|\theta)P(\theta)d\theta$: normalizing constant (often intractable)

Conjugate Priors and Closed-Form Posteriors

When the prior and likelihood are conjugate, the posterior has the same distributional family as the prior — yielding an analytic solution.

Likelihood	Conjugate Prior	Posterior
Bernoulli (coin flip)	Beta($\alpha, \beta$)	Beta($\alpha + \text{heads}, \beta + \text{tails}$)
Poisson (counts)	Gamma($\alpha, \beta$)	Gamma($\alpha + \sum y_i, \beta + n$)
Gaussian (known $\sigma$)	Gaussian($\mu_0, \sigma_0^2$)	Gaussian (weighted combination)
Multinomial	Dirichlet	Dirichlet (add counts)

Example: Coin Flip (Beta-Binomial)

Observe 7 heads in 10 flips. Prior: $\text{Beta}(2, 2)$ (weakly biased toward 0.5):

$\theta \mid \mathcal{D} \sim \text{Beta}(2 + 7,\ 2 + 3) = \text{Beta}(9, 5)$

Posterior mean: $\frac{9}{9+5} \approx 0.64$ (pulled toward prior — regularization effect).

Bayesian regularization: a Gaussian prior on weights $\equiv$ L2 regularization. A Laplace prior $\equiv$ L1 (LASSO).

Gaussian Processes (GPs)

A Gaussian Process is a distribution over functions:

$f(\cdot) \sim \mathcal{GP}(m(\cdot),\ k(\cdot, \cdot))$

• $m(x) = \mathbb{E}[f(x)]$: mean function (often set to 0)
• $k(x, x') = \text{Cov}(f(x), f(x'))$: kernel / covariance function — encodes assumptions about smoothness, periodicity, etc.

Key insight: any finite collection of function values is jointly Gaussian: $\mathbf{f} = [f(x_1), \ldots, f(x_n)]^\top \sim \mathcal{N}(\mathbf{m}, \mathbf{K})$

where $K_{ij} = k(x_i, x_j)$.

GP Regression (Prediction)

Given observations $\mathcal{D} = \{(x_i, y_i)\}$ with $y_i = f(x_i) + \varepsilon_i$, $\varepsilon_i \sim \mathcal{N}(0, \sigma_n^2)$:

$f_* \mid X_*, X, \mathbf{y} \sim \mathcal{N}(\boldsymbol{\mu}_*, \boldsymbol{\Sigma}_*)$

$\boldsymbol{\mu}_* = K_{*X}(K_{XX} + \sigma_n^2 I)^{-1}\mathbf{y}$

$\boldsymbol{\Sigma}_* = K_{**} - K_{*X}(K_{XX} + \sigma_n^2 I)^{-1}K_{X*}$

The posterior gives a predictive mean and uncertainty bands at every test point.

Common Kernels

Kernel	Assumption	Use Case
RBF / Squared Exponential	Infinitely smooth	Default; smooth functions
Matérn 5/2	Twice-differentiable	More realistic; robust
Periodic	Repeating patterns	Seasonal data
Linear	Linear functions	Bayesian linear regression
Rational Quadratic	Multi-scale smoothness	Mixture of RBF scales

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, Matern, WhiteKernel, ConstantKernel

# ---- True latent function (unknown to the model) ----
def f_true(x):
    return np.sin(3 * x) + 0.5 * x

np.random.seed(42)
X_train = np.sort(np.random.uniform(0, 5, 10)).reshape(-1, 1)
y_train = f_true(X_train.ravel()) + np.random.normal(0, 0.3, 10)

X_test = np.linspace(-0.5, 5.5, 200).reshape(-1, 1)

# ---- GP with RBF kernel ----
kernel = ConstantKernel(1.0) * RBF(length_scale=1.0) + WhiteKernel(noise_level=0.1)
gp = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=10, normalize_y=True)
gp.fit(X_train, y_train)

mu, sigma = gp.predict(X_test, return_std=True)
print(f"Optimized kernel: {gp.kernel_}")
print(f"Log-marginal-likelihood: {gp.log_marginal_likelihood(gp.kernel_.theta):.3f}")

# ---- GP with Matérn kernel (more robust) ----
kernel_m = ConstantKernel(1.0) * Matern(length_scale=1.0, nu=2.5) + WhiteKernel(0.1)
gp_m = GaussianProcessRegressor(kernel=kernel_m, n_restarts_optimizer=10, normalize_y=True)
gp_m.fit(X_train, y_train)
mu_m, sigma_m = gp_m.predict(X_test, return_std=True)

# ---- Draw posterior samples ----
samples = gp.sample_y(X_test, n_samples=5, random_state=1)

print(f"
At x=2.5:")
idx = np.argmin(np.abs(X_test.ravel() - 2.5))
print(f"  Predictive mean:   {mu[idx]:.3f}")
print(f"  Predictive std:    {sigma[idx]:.3f}")
print(f"  95% CI: [{mu[idx]-1.96*sigma[idx]:.3f}, {mu[idx]+1.96*sigma[idx]:.3f}]")
print(f"  True value:        {f_true(2.5):.3f}")

# Uncertainty increases in regions far from training data (x < 0, x > 5)
far_idx = np.argmin(np.abs(X_test.ravel() - 5.3))
print(f"
At x=5.3 (extrapolation):")
print(f"  Predictive std:    {sigma[far_idx]:.3f}  (should be higher)")

Bayesian Optimization (BO)

Bayesian Optimization efficiently optimizes expensive black-box functions — those with no gradient, slow to evaluate, and potentially noisy.

$x^* = \arg\max_{x \in \mathcal{X}} f(x)$

Classic use: hyperparameter tuning where $f$ = validation accuracy after full training run.

The BO Loop

1. Initialize: evaluate f at a few random points
2. Fit a surrogate model (usually a GP) to observed {x_i, f(x_i)}
3. Maximize acquisition function α(x) to select next x_candidate
4. Evaluate f(x_candidate) (the expensive step)
5. Update surrogate with new observation
6. Repeat until budget exhausted

Acquisition Functions

Function	Idea	Exploration vs Exploitation
Expected Improvement (EI)	$\mathbb{E}[\max(f(x) - f^+, 0)]$	Balanced (most popular)
Upper Confidence Bound (UCB)	$\mu(x) + \kappa \sigma(x)$	Tunable via $\kappa$
Probability of Improvement (PI)	$P(f(x) > f^+)$	More exploitative
Thompson Sampling	Sample from posterior, pick max	Randomized exploration

EI is preferred in practice — it naturally balances exploration (high $\sigma$) and exploitation (high $\mu$) without a tuning parameter.

import numpy as np
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import Matern
from scipy.stats import norm
from scipy.optimize import minimize

# ---- Black-box function (simulate expensive ML training) ----
def black_box(x):
    """Noisy function with global max near x=2.2"""
    return -(x - 2.2)**2 + np.sin(5*x) + np.random.normal(0, 0.05)

# ---- Expected Improvement acquisition ----
def expected_improvement(x, gp, y_best, xi=0.01):
    x = np.atleast_2d(x)
    mu, sigma = gp.predict(x, return_std=True)
    z = (mu - y_best - xi) / (sigma + 1e-9)
    ei = (mu - y_best - xi) * norm.cdf(z) + sigma * norm.pdf(z)
    return -ei.ravel()  # minimize negative EI

# ---- Bayesian Optimization loop ----
np.random.seed(7)
X_obs = np.random.uniform(0, 5, 3).reshape(-1, 1)
y_obs = np.array([black_box(x[0]) for x in X_obs])

kernel = Matern(nu=2.5)
gp = GaussianProcessRegressor(kernel=kernel, alpha=1e-6, n_restarts_optimizer=5)

n_iter = 15
for i in range(n_iter):
    gp.fit(X_obs, y_obs)
    y_best = y_obs.max()

    # Maximize EI by multi-start minimization
    best_x, best_ei = None, np.inf
    for x0 in np.random.uniform(0, 5, 20):
        result = minimize(expected_improvement, x0, args=(gp, y_best),
                          bounds=[(0, 5)], method='L-BFGS-B')
        if result.fun < best_ei:
            best_ei = result.fun
            best_x = result.x[0]

    y_new = black_box(best_x)
    X_obs = np.vstack([X_obs, [[best_x]]])
    y_obs = np.append(y_obs, y_new)
    print(f"Iter {i+1:2d}: x={best_x:.3f}, f(x)={y_new:.3f}, best so far={y_obs.max():.3f}")

print(f"
BO found maximum at x={X_obs[y_obs.argmax(), 0]:.3f} (true: ~2.2)")

# ---- Same thing with Optuna (production BO) ----
import optuna
optuna.logging.set_verbosity(optuna.logging.WARNING)

def objective(trial):
    x = trial.suggest_float('x', 0, 5)
    return -(x - 2.2)**2 + np.sin(5*x)  # no noise for Optuna demo

study = optuna.create_study(direction='maximize',
                            sampler=optuna.samplers.TPESampler(seed=42))
study.optimize(objective, n_trials=20)
print(f"
Optuna best: x={study.best_params['x']:.3f}, f={study.best_value:.3f}")

Approximate Bayesian Inference

For complex models, the posterior $P(\theta | \mathcal{D})$ is intractable. Two major approximation families:

Markov Chain Monte Carlo (MCMC)

Samples from the posterior by constructing a Markov chain whose stationary distribution is $P(\theta | \mathcal{D})$.

• Metropolis-Hastings: propose → accept/reject based on likelihood ratio
• Hamiltonian Monte Carlo (HMC): uses gradient information for efficient exploration
• NUTS (No-U-Turn Sampler): adaptive HMC — used in PyMC and Stan
• Pros: asymptotically exact, full posterior
• Cons: slow for large datasets/models, convergence diagnostics needed

Variational Inference (VI)

Approximate the posterior with a simpler distribution $q_\phi(\theta)$ from a tractable family by minimizing KL divergence:

$q^*(\theta) = \arg\min_{q \in \mathcal{Q}} \text{KL}(q(\theta) \| P(\theta | \mathcal{D}))$

Equivalently, maximize the ELBO (Evidence Lower BOund):

$\text{ELBO} = \mathbb{E}_{q}[\log P(\mathcal{D} | \theta)] - \text{KL}(q(\theta) \| P(\theta))$

• Mean-field VI: $q(\theta) = \prod_i q_i(\theta_i)$ — factored, fast
• Stochastic VI: mini-batch gradients for large data
• Pros: fast, scalable, GPU-friendly (used in VAEs, Bayesian neural nets)
• Cons: underestimates uncertainty, mode-seeking behavior

Practical Summary

Method	Accuracy	Speed	Scalability
Conjugate Bayes	Exact	Instant	Small models only
MCMC (NUTS)	Near-exact	Slow	Medium data
Variational Inference	Approximate	Fast	Large data/models
Laplace Approximation	Approximate	Fast	Point estimate + curvature

Summary

• Bayesian learning: treat parameters as distributions; update from prior → posterior using Bayes' theorem
• Conjugate priors give analytic posteriors; Gaussian prior on weights = L2 regularization
• Gaussian Processes: nonparametric Bayesian models — distributions over functions with principled uncertainty quantification
• Bayesian Optimization: use a GP surrogate + acquisition function to efficiently optimize expensive black-box functions (hyperparameter tuning, NAS, AutoML)
• MCMC: asymptotically exact posterior sampling — use when precision matters and data is moderate
• Variational Inference: fast approximate inference — use for large-scale models and deep learning

Bayesian methods underpin many modern ML systems: uncertainty-aware predictions, active learning, AutoML, and safe reinforcement learning.