When exact posterior computation is infeasible, Markov chain Monte Carlo (MCMC) provides a way to sample from a target distribution \(\pi\). A Markov chain satisfies the Markov property \(P(X_t \mid X_{1:t-1}) = P(X_t \mid X_{t-1})\) — the future depends on the past only through the present. A stationary distribution \(\pi\) satisfies \(\pi = \pi P\); under irreducibility and aperiodicity, \(\pi\) exists, is unique, and equals the long-run fraction of time the chain spends in each state. The detailed balance condition \(\pi_i P_{ij} = \pi_j P_{ji}\) is sufficient (with irreducibility) to guarantee that \(\pi\) is the stationary distribution. The Metropolis–Hastings algorithm exploits this by proposing \(y\) from \(q(y \mid x)\) and accepting with probability \(\min\left(1, \frac{\pi(y) q(x \mid y)}{\pi(x) q(y \mid x)}\right)\). Gibbs sampling is a special case that samples each variable from its full conditional given the others, with proposal equal to the conditional and acceptance probability one. Hamiltonian Monte Carlo (HMC) uses gradient information to propose distant moves that are accepted with high probability, making it effective for high-dimensional, smooth targets.
Gradients through expectations under distributions are central to training probabilistic models. The score function is \(\nabla_\theta \log p_\theta(x)\); the REINFORCE estimator uses this to obtain \(\nabla_\theta E_{p_\theta}[f(x)] = E_{p_\theta}[f(x) \nabla_\theta \log p_\theta(x)]\), which is unbiased but typically high-variance. Variance can be reduced with a baseline (control variate) \(b(x)\), since \(E[\nabla_\theta \log p_\theta(x) \cdot b(x)] = 0\) when \(b\) does not depend on \(\theta\) or on the sampling randomness. The reparameterized gradient estimator instead writes \(x = g_\theta(\varepsilon)\) with \(\varepsilon \sim q(\varepsilon)\) noise-free, giving \(\nabla_\theta E_q[f(x)] = E_\varepsilon[\nabla_\theta f(g_\theta(\varepsilon))]\), which often has substantially lower variance. Importance sampling offers another route: \(E_{p(x)}[f(x)] \approx \frac{1}{N}\sum_i f(x_i) p(x_i)/q(x_i)\) with samples drawn from a proposal \(q\).
These tools underpin a family of modern generative models. The variational autoencoder (VAE) is a latent-variable model \(p_\theta(x \mid z)\) with \(z \sim N(0, I)\) prior, trained by maximizing the ELBO via an inference network \(q_\phi(z \mid x)\) and gradients obtained through the reparameterization trick. Normalizing flows compose invertible, differentiable transformations \(f = f_K \circ \cdots \circ f_1\) that map a simple base density like \(N(0, I)\) to a complex one, with log-densities computed via the change-of-variables formula \(\log p(x) = \log p_z(f^{-1}(x)) + \sum_i \log |\det(\partial f_i^{-1}/\partial x)|\). Diffusion models learn to reverse a fixed noising Markov chain that gradually adds Gaussian noise to data; training minimizes a denoising score-matching loss on the noise-conditional distributions \(p(x_t \mid x_0) \sim N(\sqrt{\bar{\alpha}_t} x_0, (1 - \bar{\alpha}_t) I)\). Score-matching objectives match \(\nabla_x \log p_\theta(x)\) to the data score \(\nabla_x \log p_{\text{data}}(x)\); denoising score matching provides a tractable surrogate by corrupting data with Gaussian noise and matching the score of the noise-conditional model. Energy-based models (EBMs) take the form \(p_\theta(x) = e^{-E_\theta(x)} / Z(\theta)\), generalizing the Boltzmann (Gibbs) distribution \(p(x) = e^{-E(x)/T}/Z\) from physics to machine learning; restricted Boltzmann machines (RBMs) are a two-layer undirected version with visible and binary hidden units, trained by contrastive divergence. Latent Dirichlet Allocation (LDA) is a generative topic model in which each document has a topic-mixing vector \(\theta_d \sim\) Dirichlet(\(\alpha\)), each topic has a word distribution \(\phi_k \sim\) Dirichlet(\(\beta\)), and each word is drawn by choosing \(z \sim\) Categorical(\(\theta_d\)) and \(w \sim\) Categorical(\(\phi_z\)). Finally, density estimation can be non-parametric via the kernel density estimator (KDE), which places a smooth kernel \(K\) (with bandwidth \(h\)) on each data point: \(\hat{f}(x) = \frac{1}{nh}\sum_i K((x - x_i)/h)\). The bandwidth controls the classic bias–variance tradeoff: smoother models yield higher bias but lower variance, while complex models reduce bias at the cost of higher variance and overfitting risk. AIC and BIC provide penalized log-likelihood criteria for model selection: \(\text{AIC} = 2k - 2\log \hat{L}\) and \(\text{BIC} = k \log n - 2\log \hat{L}\), with BIC penalizing complexity more strongly as sample size \(n\) grows.