Probability Distributions for ML

This deck walks you through the foundational probability distributions you'll encounter in machine learning, starting with the basics of what a probability distribution actually is and how discrete and continuous variables are handled differently. You'll get clear definitions of core tools like the probability mass function, probability density function, and cumulative distribution function, along with how they relate to one another. From there, it moves into summary statistics such as expected value, variance, and standard deviation, before introducing the specific distributions most relevant to ML: Bernoulli, Binomial, Categorical, and Multinomial.

It's a great fit if you're a machine learning learner, data science student, or practitioner who wants to strengthen the mathematical foundation underneath algorithms like logistic regression, naive Bayes, or generative models. If you've ever felt shaky about why a sigmoid outputs a probability, or what "sampling from a distribution" really means, working through these cards will make those ideas much more intuitive.

Because the topics build on each other, going through the deck in order the first time can help you see the connections between distributions rather than treating them as isolated formulas. After that, spaced repetition works really well here: review the cards over several days or weeks rather than cramming, and try to explain each concept in your own words before flipping the card. Linking each distribution back to a real ML use case, even a simple one, will also help the material stick far longer than memorizing formulas alone.

A probability distribution is a rule that assigns probabilities to outcomes in a sample space, with the total probability summing (for discrete cases) or integrating (for continuous cases) to one. For discrete random variables, this rule is called a probability mass function (PMF), written \(p(x)\), which must be non-negative and satisfy \(\sum_x p(x) = 1\). For continuous random variables, the analogous concept is the probability density function (PDF), written \(f(x)\), such that the probability of the variable falling in any interval is \(P(a \le X \le b) = \int_a^b f(x)\,dx\). Both PMFs and PDFs are non-negative and integrate or sum to one, but only the PDF can take values greater than one, because it represents density rather than probability directly.

Beyond the immediate description of probabilities, distributions can be characterized through cumulative and aggregated quantities. The cumulative distribution function (CDF) is defined as \(F(x) = P(X \le x)\); it is non-decreasing, right-continuous, starts at 0 as \(x \to -\infty\), and approaches 1 as \(x \to \infty\). For continuous random variables, the CDF is the integral of the PDF up to \(x\), and where \(f\) is continuous, the PDF is recovered as the derivative of the CDF: \(f(x) = dF/dx\). The expected value (or mean) summarizes central tendency: for discrete \(X\) it is \(E[X] = \sum_x x \cdot p(x)\), and for continuous \(X\) it is \(E[X] = \int_{-\infty}^{\infty} x \cdot f(x)\,dx\), with absolute convergence assumed.

Dispersion around the mean is captured by the variance, defined as \(\mathrm{Var}(X) = E[(X - E[X])^2]\), equivalently \(E[X^2] - (E[X])^2\). The standard deviation is the square root of the variance, restoring the original units of \(X\). These summary statistics — mean, variance, and the underlying CDF — provide the basic vocabulary for describing any distribution, whether discrete or continuous, single-variable or multivariate.

A family of closely related discrete distributions arises from independent Bernoulli trials. The Bernoulli distribution itself models a single binary outcome: \(P(X = 1) = p\), \(P(X = 0) = 1 - p\), with mean \(p\) and variance \(p(1 - p)\). The Binomial distribution counts the number of successes across \(n\) independent Bernoulli(\(p\)) trials, with PMF \(P(X = k) = \binom{n}{k} p^k (1-p)^{n-k}\) for \(k = 0, 1, \dots, n\), mean \(np\), and variance \(np(1-p)\). Generalizing from two outcomes to \(K\) categories gives the Categorical (or Multinoulli) distribution, which samples a single category with probabilities summing to one, and the Multinomial distribution, which counts outcomes across \(K\) categories in \(n\) independent draws with joint PMF \(\frac{n!}{\prod_i x_i!} \prod_i p_i^{x_i}\).

Several other discrete distributions handle waiting-time and event-counting problems. The Geometric distribution (in its "number-of-trials" form) gives the trial index of the first success, \(P(X = k) = (1-p)^{k-1} p\) for \(k = 1, 2, \dots\), with mean \(1/p\); the equivalent "failures-before-success" form is \(P(X = k) = (1-p)^k p\) starting at \(k = 0\), with mean \((1-p)/p\). Extending this idea, the Negative Binomial distribution counts failures before the \(r\)-th success, with PMF \(P(X = k) = \binom{r+k-1}{k} p^r (1-p)^k\), mean \(r(1-p)/p\), and variance \(r(1-p)/p^2\). The Poisson distribution models counts of independent events at constant rate \(\lambda\) over a fixed interval, \(P(X = k) = e^{-\lambda} \lambda^k / k!\), with both mean and variance equal to \(\lambda\). It is a good model for rare events such as decay events, web page hits per minute, or arrival counts.

The Poisson distribution emerges naturally as a limit of the Binomial: when \(p \to 0\) and \(n \to \infty\) with \(np = \lambda\) held fixed, \(P(X = k) \approx e^{-\lambda} \lambda^k / k!\). This justifies using Poisson as the "rare events" approximation. As a useful implementation detail, \(k\) independent Bernoulli(0.5) samples can be combined into a uniformly distributed integer in \(\{0, \dots, 2^k - 1\}\) by interpreting the bit string as binary — an example of how simple distributions compose to produce more complex ones.

Many continuous distributions arise from natural generative processes. The Uniform distribution on \([a, b]\) has constant density \(f(x) = 1/(b - a)\) over the interval, with mean \((a + b)/2\) and variance \((b - a)^2/12\). The Exponential distribution with rate \(\lambda\) models the time between events in a Poisson process; its PDF is \(f(t) = \lambda e^{-\lambda t}\) for \(t \ge 0\), with mean \(1/\lambda\) and variance \(1/\lambda^2\). A defining property of the Exponential is its memorylessness: \(P(X > s + t \mid X > s) = P(X > t)\), meaning that the distribution of remaining waiting time is independent of how long one has already waited. The Gamma distribution generalizes the Exponential via a shape parameter \(k > 0\) and rate \(\theta > 0\) (or scale \(\beta = 1/\theta\)), with PDF \(f(x) = x^{k-1} e^{-x/\theta} / (\Gamma(k) \theta^k)\), mean \(k\theta\), and variance \(k\theta^2\). Setting \(k = 1\) recovers the Exponential, and Gamma(\(k, \theta\)) models the sum of \(k\) independent Exponential(\(1/\theta\)) waiting times.

A second cluster of continuous distributions lives naturally on bounded intervals or the positive reals. The Beta distribution on \((0, 1)\) has PDF \(f(p) = p^{\alpha - 1}(1 - p)^{\beta - 1} / B(\alpha, \beta)\), with mean \(\alpha/(\alpha + \beta)\); it is widely used in ML because it serves as the conjugate prior for Bernoulli, Binomial, and Geometric likelihoods, and provides a flexible family for modeling probabilities. The Dirichlet distribution extends Beta to the \(K\)-simplex, with density \(f(\pi) = (1/B(\alpha)) \prod_i \pi_i^{\alpha_i - 1}\) for \(\pi_i \ge 0\) with \(\sum \pi_i = 1\); it is the conjugate prior of the Multinomial and is central to Bayesian topic models. The Log-Normal distribution describes positive random variables whose logarithm is Normal, with mean \(e^{\mu + \sigma^2/2}\), and is appropriate when multiplicative rather than additive variation dominates.

Other continuous families are important in ML for modeling heavy-tailed behavior, sparse signals, or specialized physical processes. The Student's t-distribution has heavier tails than the Normal, with mean 0 (for \(\nu > 1\)) and variance \(\nu/(\nu - 2)\) for \(\nu > 2\), converging to Normal as \(\nu \to \infty\). The Cauchy (or Lorentz) distribution is the t with \(\nu = 1\); it has undefined mean and variance and is a classic heavy-tailed example. The Laplace (double-exponential) distribution \(f(x \mid \mu, b) = (1/(2b)) e^{-|x - \mu|/b}\) has heavier tails than Normal and is often used as a sparsity-inducing prior on regression coefficients. Heavy-tailed distributions in general decay more slowly than exponential — for instance as power laws \(P(X > x) \propto x^{-\alpha}\), exemplified by the Pareto distribution \(f(x) = \alpha x_m^{\alpha} / x^{\alpha + 1}\) for \(x \ge x_m\), which underlies the "80/20 rule" in wealth and city-size modeling, and the discrete Zipf distribution where the \(k\)-th most frequent item has probability proportional to \(1/k^s\). The chi-squared distribution arises as the sum of squared standard normals and is itself a Gamma(\(k/2, 2\)); in ratio form it gives the F-distribution used in ANOVA. Survival and hazard functions provide an alternative language for non-negative distributions: the survival function \(S(t) = P(T > t) = 1 - F(t)\), and the hazard \(h(t) = f(t)/S(t)\) is the instantaneous failure rate conditional on survival. For the Weibull distribution with shape \(k\), the hazard decreases for \(k < 1\) (infant mortality), is constant for \(k = 1\) (Exponential), and increases for \(k > 1\) (wear-out); the Rayleigh distribution, a Weibull with shape 2, models magnitudes of 2D Gaussian vectors, while the Maxwell–Boltzmann distribution, a chi distribution with 3 degrees of freedom, describes molecular speeds.

The Normal (Gaussian) distribution is the centerpiece of continuous probability for ML, with PDF \(f(x) = \frac{1}{\sigma\sqrt{2\pi}} \exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right)\), parameterized by mean \(\mu\) and variance \(\sigma^2\). The standard Normal with \(\mu = 0\) and \(\sigma^2 = 1\) has CDF denoted \(\Phi\) and PDF \(\phi(z) = \frac{1}{\sqrt{2\pi}} e^{-z^2/2}\). A useful heuristic is the 68-95-99.7 rule: roughly 68% of the probability mass lies within 1\(\sigma\) of the mean, 95% within 2\(\sigma\), and 99.7% within 3\(\sigma\). The Central Limit Theorem (CLT) explains the Normals ubiquity: the sum (properly normalized) of \(n\) i.i.d. random variables with finite mean \(\mu\) and variance \(\sigma^2\) converges in distribution to Normal(\(\mu, \sigma^2/n\)) as \(n\) grows. This justifies treating averages, losses, and estimator errors as approximately Normal in machine learning, enabling confidence intervals, z-tests, and Gaussian approximations of aggregate behavior.

Sampling from a standard Normal is often done via the Box–Muller transform, which maps two independent U(0, 1) uniforms to two independent standard Normals: \[Z_0 = \sqrt{-2 \ln U_1} \cos(2\pi U_2), \quad Z_1 = \sqrt{-2 \ln U_1} \sin(2\pi U_2).\] The multivariate extension, the multivariate Gaussian, places a \(d\)-dimensional Normal with density \(f(x) = (2\pi)^{-d/2} |\Sigma|^{-1/2} \exp\left(-\tfrac{1}{2}(x - \mu)^T \Sigma^{-1} (x - \mu)\right)\), where \(\Sigma\) must be positive definite. Distances in this space are measured by the Mahalanobis distance \(d(x, \mu) = \sqrt{(x - \mu)^T \Sigma^{-1} (x - \mu)}\), which scales deviations by the covariance structure. An isotropic Gaussian has \(\Sigma = \sigma^2 I_d\), making its density depend only on Euclidean distance and giving it rotational symmetry. A singular (or degenerate) Gaussian has a singular \(\Sigma\), so its mass lies on a lower-dimensional affine subspace and the density is undefined in the ambient space. Whitening applies a linear transformation \(z = L^{-1}(x - \mu)\), where \(L\) is typically the Cholesky factor of \(\Sigma\) with \(LL^T = \Sigma\), mapping a centered multivariate Normal to a standard Normal.

Several matrix-valued distributions extend the Gaussian framework to covariance estimation and Bayesian modeling. The Wishart distribution generalizes chi-squared to matrices, governing the distribution of sums of outer products of multivariate Gaussian vectors; it serves as a conjugate prior for the inverse covariance of a multivariate Normal. The inverse-Wishart plays the conjugate role for the covariance \(\Sigma\) itself: with \(\Sigma \sim\) Inv-Wishart(\(\Psi, \nu\)) and i.i.d. multivariate Normal data, the posterior is also inverse-Wishart. The LKJ prior is a distribution over correlation matrices indexed by a shape parameter \(\eta \ge 0\), with \(\eta = 1\) producing a uniform distribution over correlation matrices and larger \(\eta\) concentrating mass near the identity. Finally, the reparameterization trick — writing \(x = \mu + \sigma \odot \varepsilon\) with \(\varepsilon \sim N(0, I)\) — allows gradients of expectations \(E_{x \sim N(\mu, \sigma^2)}[f(x)]\) to flow through \(\mu\) and \(\sigma\), and is foundational to variational autoencoders and other gradient-based generative models.

Parameter estimation can be framed either in frequentist or Bayesian terms. Maximum likelihood estimation (MLE) chooses parameters \(\theta\) to maximize the log-likelihood \(L(\theta) = \sum_n \log p(x_n \mid \theta)\), equivalently minimizing the KL divergence between the empirical distribution and \(p(\cdot \mid \theta)\) up to a constant. Maximum a posteriori (MAP) estimation adds a prior term: it maximizes \(\log p(\theta \mid D) \propto \log p(D \mid \theta) + \log p(\theta)\). For instance, when estimating a Gaussian mean with known variance \(\sigma^2\) under a Gaussian prior \(N(\mu_0, \tau^2)\), MLE returns the sample mean \(\bar{x}\), while MAP returns the precision-weighted shrinkage estimator \((n \tau^2 \bar{x} + \sigma^2 \mu_0) / (n \tau^2 + \sigma^2)\), interpolating between data and prior. Frequentists treat the true parameter as fixed and reason about sampling distributions, whereas Bayesians treat parameters as random variables, place priors over them, and update to posteriors via Bayes' rule, yielding full predictive distributions.

Conjugate priors are a cornerstone of Bayesian computation: a prior \(p(\theta)\) is conjugate to a likelihood if the posterior \(p(\theta \mid D)\) lies in the same family as the prior. Standard pairings include Beta–Bernoulli, Dirichlet–Multinomial, Gamma–Poisson, and Normal–Normal, as well as inverse-Wishart for multivariate Normal covariances. For example, with a Bernoulli likelihood and Beta(\(\alpha, \beta\)) prior, the posterior is Beta(\(\alpha +\) successes, \(\beta +\) failures); analogously, a Multinomial likelihood combined with a Dirichlet(\(\alpha\)) prior yields a Dirichlet(\(\alpha +\) counts) posterior, a relationship that underpins Bayesian text models like Latent Dirichlet Allocation. When conjugacy fails, posterior computation requires approximate methods: MCMC, variational inference, or Laplace / expectation-propagation approximations.

The choice of prior encodes prior beliefs and affects posterior behavior. The Jeffreys prior is non-informative in a precise sense: it is proportional to \(\sqrt{|I(\theta)|}\), where \(I(\theta)\) is the Fisher information, and is invariant under reparameterization. The posterior predictive distribution \(p(x_{\text{new}} \mid D) = \int p(x_{\text{new}} \mid \theta) p(\theta \mid D)\, d\theta\) averages likelihoods over the posterior to make predictions under uncertainty, while the prior predictive \(p(x) = \int p(x \mid \theta) p(\theta)\, d\theta\) gives the marginal distribution induced by the prior alone. Bayesian linear regression illustrates these ideas concretely: under \(y = w^T x + \varepsilon\) with \(\varepsilon \sim N(0, \sigma^2)\) and prior \(w \sim N(0, \tau^2 I)\), the posterior over weights is itself Gaussian with closed-form mean \((X^T X + (\sigma^2/\tau^2) I)^{-1} X^T y\), and the precision term \(\sigma^2/\tau^2\) acts as an L2 regularizer.

Mixture distributions combine several component distributions to model data with heterogeneous structure: an index \(z\) is drawn from Categorical(\(\pi\)), then \(x\) is drawn from \(p(x \mid z)\), giving marginal density \(p(x) = \sum_k \pi_k p_k(x)\). A Gaussian Mixture Model (GMM) is the most common instance, mixing \(K\) multivariate Normals with weights \(\pi_k\), means \(\mu_k\), and covariances \(\Sigma_k\). The component assignment \(z\) is a latent variable — an unobserved quantity that influences observed data — and is learned rather than observed directly. The standard algorithm for fitting GMMs is Expectation-Maximization (EM): the E-step computes posterior responsibilities \(\gamma_{nk} = \pi_k N(x_n; \mu_k, \Sigma_k) / \sum_j \pi_j N(x_n; \mu_j, \Sigma_j)\), and the M-step re-estimates \(\pi_k\), \(\mu_k\), \(\Sigma_k\) from \(\gamma\)-weighted sufficient statistics. EM is a general recipe for maximum likelihood in latent-variable models, and it underlies countless applications from clustering to speech recognition.

Information theory provides the tools to compare distributions and bound quantities that arise in latent-variable modeling. KL divergence \(D_{\text{KL}}(p \| q) = E_{x \sim p}[\log p(x)/q(x)]\) is non-negative, asymmetric, and not a true metric, but it captures how one distribution \(p\) differs from a reference \(q\). Jensen's inequality states that for a convex function \(\varphi\), \(\varphi(E[X]) \le E[\varphi(X)]\); this simple fact drives many variational bounds. Most importantly, the Evidence Lower Bound (ELBO) is defined as \(\text{ELBO} = E_q[\log p(x, z)] + H(q) = \log p(x) - D_{\text{KL}}(q(z) \| p(z \mid x))\); it lower-bounds the log evidence \(\log p(x)\) and is tractable whenever the joint \(p(x, z)\) and a tractable \(q\) can be evaluated, making it the workhorse of variational inference.

Variational inference (VI) approximates an intractable posterior \(p(z \mid x)\) with a tractable \(q_\phi(z)\) by minimizing \(D_{\text{KL}}(q_\phi \| p(\cdot \mid x))\), equivalently maximizing the ELBO. Mean-field VI further restricts \(q\) to a factorized form \(q(z) = \prod_i q_i(z_i)\); this leads to coordinate-ascent updates where each \(q_i\) is set to the log-normalized conditional expectation of \(\log p\) under the current values of the other \(q_j\)'s. When closed-form optimization is impossible, importance sampling offers an alternative: \(E_{p(x)}[f(x)]\) is approximated by \(\frac{1}{N} \sum_i f(x_i) p(x_i)/q(x_i)\) with samples drawn from a proposal \(q\). Variance depends critically on the choice of \(q\), which should ideally have heavier tails than \(|f| p\).

Classification algorithms can be organized by how they treat the joint distribution of features and labels. Generative classifiers model the joint \(p(x, y)\) by specifying \(p(y)\) and \(p(x \mid y)\); discriminative classifiers model \(p(y \mid x)\) directly. Naive Bayes is the canonical generative model: assuming features are conditionally independent given the class label \(y \sim\) Categorical, the joint becomes a product of per-feature likelihoods. A Bernoulli Naive Bayes uses Bernoulli(\(p_{c,i}\)) for each feature given class \(c\), while a Gaussian Naive Bayes uses \(N(x_i; \mu_{c,i}, \sigma^2_{c,i})\). Predictions are made by the maximum posterior class. Gaussian Discriminant Analysis (GDA) is a richer generative model with full multivariate Gaussian class-conditionals \(p(x \mid y = k) = N(x; \mu_k, \Sigma_k)\) and class priors \(p(y = k) = \pi_k\); the Bayes-optimal decision boundary is quadratic. Adding the assumption that all classes share a common covariance \(\Sigma\) turns GDA into Linear Discriminant Analysis (LDA), with a linear decision boundary, while removing that assumption gives Quadratic Discriminant Analysis (QDA) with a quadratic boundary. LDA is closely related to Fisher's linear discriminant, since both seek a projection that maximizes class separability — Fisher's criterion yields a closed-form projection, while LDA arises from the equal-covariance generative assumption.

Discriminative classification typically uses a softmax or sigmoid parameterization. The softmax function maps real-valued scores to a Categorical distribution via \(\text{softmax}(z_i) = e^{z_i} / \sum_j e^{z_j}\); it arises naturally from the multinomial logistic regression likelihood. The sigmoid \(\sigma(z) = 1/(1 + e^{-z})\) maps a real number to \((0, 1)\) and parameterizes Bernoulli probabilities in logistic regression. The logistic (sigmoid) distribution is the corresponding continuous distribution whose CDF equals \(\sigma\), with PDF \(f(x) = e^{-(x - \mu)/s} / (s(1 + e^{-(x - \mu)/s})^2)\), mean \(\mu\), and variance \(s^2 \pi^2 / 3\). The softmax itself can be viewed as a Gibbs distribution over discrete states with energies \(-z_i\) and temperature 1, and softmax with temperature \(T\) corresponds to \(\exp(z_i / T)\).

The corresponding loss functions are negative log-likelihoods: categorical cross-entropy \(-\sum_n \sum_k y_{n,k} \log p(y_{n,k} \mid x_n)\) for \(K\)-class softmax classifiers, binary cross-entropy \(-\frac{1}{N}\sum_n [y_n \log p(x_n) + (1 - y_n) \log(1 - p(x_n))]\) for sigmoid-parameterized Bernoullis, and MSE for a Gaussian with fixed variance up to a constant. The choice between discriminative and generative approaches reflects modeling priorities: generative classifiers can synthesize new examples and incorporate prior knowledge about how \(x\) is generated, while discriminative classifiers often achieve higher accuracy when the goal is purely prediction and the conditional \(p(y \mid x)\) is all that matters. In both cases, the parameterizations (softmax, sigmoid, multivariate Gaussian) and the associated losses connect directly back to the underlying probability distributions.

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.