At the heart of every generative model lies a seemingly simple goal: estimate the probability density function of data. Given samples $x_1, x_2, \ldots, x_n$ drawn from some unknown distribution $P^*(x)$, we want to learn an approximation $P_\theta(x)$ that captures the true data-generating process.

This is density estimation—arguably the most fundamental problem in statistics. It sounds straightforward: we have data, we want to know where that data comes from. But the apparent simplicity conceals profound challenges, especially for high-dimensional, structured data like images, audio, or text.

In this page, we develop a rigorous understanding of density estimation: what it means mathematically, how classical approaches work, where they fail, and how modern neural methods attempt to overcome these limitations.

The Density Estimation Problem

Let $X$ be a random variable taking values in some space $\mathcal{X}$ (e.g., $\mathbb{R}^d$ for continuous data). We assume data points $\mathcal{D} = {x_1, \ldots, x_n}$ are drawn i.i.d. from an unknown true distribution with density $p^*(x)$:

$$x_i \stackrel{\text{i.i.d.}}{\sim} p^*(x)$$

Our goal is to construct an estimator $\hat{p}(x)$ such that $\hat{p}$ is 'close' to $p^*$ in some meaningful sense.

What does 'close' mean?

Different notions of closeness lead to different estimation procedures:

1. KL Divergence (Maximum Likelihood): $$D_{KL}(p^* | \hat{p}) = \mathbb{E}_{x \sim p^}\left[\log \frac{p^(x)}{\hat{p}(x)}\right]$$ Minimizing KL divergence is equivalent to maximizing log-likelihood.

2. Total Variation Distance: $$TV(p^, \hat{p}) = \frac{1}{2} \int |p^(x) - \hat{p}(x)| dx$$ Upper bounds classification error.

3. Wasserstein Distance: $$W_1(p^, \hat{p}) = \inf_{\gamma \in \Pi(p^, \hat{p})} \mathbb{E}_{(x,y) \sim \gamma}[|x - y|]$$ Captures geometric similarity; used in Wasserstein GANs.

Each metric has different properties. KL divergence is asymmetric and can be infinite if densities have different supports. Wasserstein is more robust but computationally harder.

Why is density estimation hard?

The challenge scales dramatically with dimensionality. Consider:

• 1D: Estimate the distribution of human heights. We might need ~100 samples to get a reasonable histogram.

• 10D: Estimate the joint distribution of 10 medical measurements. The 'space' we need to cover grows exponentially—with 10 bins per dimension, we have $10^{10}$ cells to estimate.

• 1,000,000D: Estimate the distribution of 1000×1000 RGB images. Each image is a point in million-dimensional space. We will never have enough data to 'fill' this space.

This is the curse of dimensionality—as dimensions increase, the volume of space grows exponentially, and data becomes sparse. Every point is far from every other point. Naive density estimation becomes impossible.

Modern generative models succeed by exploiting structure: natural data lies on low-dimensional manifolds, exhibits local regularities, and follows hierarchical patterns. Deep learning finds and leverages this structure.

Parametric density estimation assumes the data comes from a known family of distributions with unknown parameters. We estimate the parameters from data.

The Parametric Model:

$$p(x) = p_\theta(x)$$

where $\theta \in \Theta$ are parameters to be learned (e.g., mean $\mu$ and variance $\sigma^2$ for a Gaussian).

Advantages:

• Sample-efficient (fewer parameters to estimate)
• Often computationally tractable
• Provides interpretable parameters

Disadvantages:

• Assumes a specific functional form
• Model misspecification leads to systematic bias
• Limited flexibility

Common Parametric Models:

Maximum Likelihood Estimation (MLE)

The workhorse of parametric estimation is MLE. Given data $\mathcal{D} = {x_1, \ldots, x_n}$, we choose parameters that maximize the probability of observing this data:

$$\hat{\theta}{MLE} = \arg\max\theta \prod_{i=1}^n p_\theta(x_i) = \arg\max_\theta \sum_{i=1}^n \log p_\theta(x_i)$$

The log transformation converts products to sums, improving numerical stability and optimization.

MLE for Gaussian:

For data from $\mathcal{N}(\mu, \sigma^2)$:

$$\hat{\mu}{MLE} = \frac{1}{n} \sum{i=1}^n x_i = \bar{x}$$

$$\hat{\sigma}^2_{MLE} = \frac{1}{n} \sum_{i=1}^n (x_i - \bar{x})^2$$

The MLE for the mean is the sample mean—intuitive and optimal. The variance estimator is biased (divide by $n-1$ for unbiased version), but consistent.

Properties of MLE:

• Consistency: $\hat{\theta}_{MLE} \to \theta^*$ as $n \to \infty$ (if model is correctly specified)
• Asymptotic efficiency: Achieves the Cramér-Rao lower bound
• Asymptotic normality: $\sqrt{n}(\hat{\theta}_{MLE} - \theta^) \to \mathcal{N}(0, I^{-1}(\theta^))$

These properties make MLE the gold standard for parametric estimation.

Non-parametric density estimation makes no fixed assumptions about the functional form of the distribution. The model complexity grows with data.

Philosophy:

Instead of assuming $p^*(x)$ belongs to a parametric family, we let the data 'speak for itself.' The density estimate is directly derived from the observed samples.

Histogram Estimator

The simplest non-parametric approach. Divide the space into bins and estimate density as the proportion of samples in each bin:

$$\hat{p}(x) = \frac{\text{count in bin containing } x}{n \cdot \text{bin width}}$$

Properties:

• Extremely simple to implement
• Density is normalized (integrates to 1)
• Discontinuous (density jumps at bin boundaries)
• Highly sensitive to bin placement and width

Kernel Density Estimation (KDE)

KDE smooths the histogram by placing a 'kernel' (smooth bump) at each data point:

$$\hat{p}(x) = \frac{1}{n h} \sum_{i=1}^n K\left(\frac{x - x_i}{h}\right)$$

where:

• $K$ is a kernel function (e.g., Gaussian: $K(u) = \frac{1}{\sqrt{2\pi}} e^{-u^2/2}$)
• $h$ is the bandwidth—a crucial hyperparameter

Bandwidth Selection:

• Too small $h$: Overfitting—density is spiky, each point creates its own peak
• Too large $h$: Underfitting—density is overly smooth, loses detail
• Optimal $h$: Balances bias (smoothness) and variance (sensitivity to individual points)

Common approaches:

• Silverman's rule: $h = 0.9 \cdot \min(\hat{\sigma}, IQR/1.34) \cdot n^{-1/5}$ (for 1D Gaussian kernel)
• Cross-validation: Leave-one-out likelihood optimization

k-Nearest Neighbors Density Estimation

An alternative non-parametric approach: instead of fixing bandwidth, fix the number of neighbors:

$$\hat{p}(x) = \frac{k}{n \cdot V_d(r_k(x))}$$

where $r_k(x)$ is the distance to the $k$-th nearest neighbor, and $V_d(r)$ is the volume of a $d$-dimensional ball of radius $r$.

Intuition: In dense regions, the $k$-th neighbor is close, so the ball has small volume, giving high density. In sparse regions, the ball is large, giving low density.

Properties:

• Adaptive bandwidth (effectively)
• Works better in regions with varying density
• Still suffers from curse of dimensionality
• Not guaranteed to integrate to 1

Comparison of Approaches:

Mixture models combine parametric flexibility with non-parametric-like expressiveness. They model the distribution as a weighted sum of simpler component distributions:

$$p(x) = \sum_{k=1}^K \pi_k \cdot p_k(x|\theta_k)$$

where:

• $K$ is the number of components (clusters)
• $\pi_k$ are mixing weights ($\sum_k \pi_k = 1$, $\pi_k \geq 0$)
• $p_k(x|\theta_k)$ are component densities

Gaussian Mixture Models (GMMs):

The most common mixture model uses Gaussian components:

$$p(x) = \sum_{k=1}^K \pi_k \cdot \mathcal{N}(x | \mu_k, \Sigma_k)$$

Why mixtures are powerful:

1. Universal approximation: With enough components, GMMs can approximate any continuous density arbitrarily well.
2. Interpretable structure: Each component can represent a distinct 'mode' or cluster in the data.
3. Latent variable interpretation: The mixture can be viewed as a generative process:
   • First, sample component $z \sim \text{Categorical}(\pi_1, \ldots, \pi_K)$
   • Then, sample $x \sim p_z(x|\theta_z)$

This latent variable view foreshadows more complex generative models.

Expectation-Maximization (EM) Algorithm

Directly maximizing likelihood for mixtures is hard—we don't know which component generated each point. EM solves this by alternating:

E-step (Expectation): Compute soft assignments—probability that each point came from each component:

$$\gamma_{ik} = \frac{\pi_k \cdot p_k(x_i|\theta_k)}{\sum_j \pi_j \cdot p_j(x_i|\theta_j)}$$

M-step (Maximization): Update parameters using weighted statistics:

$$\pi_k^{new} = \frac{1}{n} \sum_i \gamma_{ik}$$

$$\mu_k^{new} = \frac{\sum_i \gamma_{ik} x_i}{\sum_i \gamma_{ik}}$$

$$\Sigma_k^{new} = \frac{\sum_i \gamma_{ik} (x_i - \mu_k^{new})(x_i - \mu_k^{new})^T}{\sum_i \gamma_{ik}}$$

Properties of EM:

• Guarantees likelihood increases (or stays same) each iteration
• Converges to local maximum (not necessarily global)
• Simple to implement, but may be slow to converge
• Sensitive to initialization

Classical density estimation methods—histograms, KDE, GMMs—work well in low dimensions but fail catastrophically for high-dimensional data like images. Neural density estimation uses deep learning to learn expressive density models that can scale to complex data.

The Core Challenge:

To define a valid probability density, we need:

1. $p(x) \geq 0$ for all $x$ (non-negativity)
2. $\int p(x) dx = 1$ (normalization)

For neural networks outputting arbitrary values, ensuring these constraints is non-trivial.

Approach 1: Autoregressive Models

Factor the joint density using the chain rule of probability:

$$p(x) = p(x_1) \cdot p(x_2|x_1) \cdot p(x_3|x_1, x_2) \cdots = \prod_{i=1}^d p(x_i | x_{<i})$$

Each conditional $p(x_i | x_{<i})$ is modeled by a neural network. If each conditional is a valid density (e.g., categorical for discrete, Gaussian for continuous), the product is automatically valid.

Examples:

• PixelCNN/PixelRNN: Model image pixels left-to-right, top-to-bottom
• WaveNet: Model audio samples autoregressively
• GPT: Model text tokens autoregressively

Advantages: Exact likelihood computation, powerful expressiveness Disadvantages: Sequential generation (slow), arbitrary ordering choice

Approach 2: Normalizing Flows

Transform a simple base distribution (e.g., standard Gaussian) through invertible transformations to create a complex distribution:

$$z \sim \mathcal{N}(0, I) \quad \rightarrow \quad x = f_\theta(z)$$

If $f$ is invertible and differentiable, the density of $x$ follows from the change of variables formula:

$$p_X(x) = p_Z(f^{-1}(x)) \cdot \left|\det \frac{\partial f^{-1}}{\partial x}\right|$$

The Jacobian determinant accounts for how $f$ stretches or compresses space.

Key insight: By stacking many simple invertible transformations, we can build arbitrarily complex distributions while maintaining tractable density evaluation.

Examples: RealNVP, Glow, Neural Spline Flows

Advantages: Exact likelihood, efficient sampling, latent space Disadvantages: Architectural constraints (must be invertible), Jacobian computation can be costly

Approach 3: Energy-Based Models

Define an unnormalized density (energy function) and train without computing the normalizing constant:

$$p_\theta(x) = \frac{\exp(-E_\theta(x))}{Z_\theta} \quad \text{where} \quad Z_\theta = \int \exp(-E_\theta(x)) dx$$

The energy $E_\theta(x)$ can be any neural network—no constraints! Low energy = high probability.

Problem: $Z_\theta$ is intractable (integral over all possible $x$).

Solutions:

• Contrastive Divergence: Approximate gradient using MCMC samples
• Score Matching: Match gradients of log-density instead
• Noise Contrastive Estimation: Discriminate real from noise samples

Deep connection: Score-based diffusion models (state-of-the-art image generation) are intimately related to energy-based models through the score function $\nabla_x \log p(x)$.

Approach 4: Variational Autoencoders (VAEs)

Learn a latent variable model by optimizing a tractable lower bound on log-likelihood:

$$\log p(x) \geq \mathbb{E}{q(z|x)}[\log p(x|z)] - D{KL}(q(z|x) | p(z))$$

VAEs will be covered in depth in Module 2.

Maximum likelihood estimation is so fundamental to density estimation that it deserves deeper examination. Understanding its properties, limitations, and connections to information theory is essential for modern generative modeling.

The Likelihood Function:

Given data $\mathcal{D} = {x_1, \ldots, x_n}$ and model $p_\theta$, the likelihood is:

$$L(\theta) = \prod_{i=1}^n p_\theta(x_i)$$

The log-likelihood is:

$$\ell(\theta) = \sum_{i=1}^n \log p_\theta(x_i)$$

MLE as KL Minimization:

A beautiful connection: minimizing KL divergence from true distribution to model is equivalent to maximizing likelihood:

$$D_{KL}(p^* | p_\theta) = \mathbb{E}{x \sim p^}[\log p^(x)] - \mathbb{E}{x \sim p^*}[\log p_\theta(x)]$$

The first term is the (negative) entropy of the true distribution—constant with respect to $\theta$. Minimizing KL is thus equivalent to maximizing $\mathbb{E}{x \sim p^*}[\log p\theta(x)]$, which is estimated by the sample log-likelihood.

Cross-Entropy Loss:

In practice, we minimize the cross-entropy:

$$H(p^, p_\theta) = -\mathbb{E}_{x \sim p^}[\log p_\theta(x)] \approx -\frac{1}{n}\sum_{i=1}^n \log p_\theta(x_i)$$

This is exactly the negative log-likelihood (per sample). Whether you call it 'maximizing likelihood' or 'minimizing cross-entropy,' it's the same optimization.

MLE Pathologies:

Despite its elegance, MLE has important limitations:

1. Mode-Seeking Behavior

Minimizing $D_{KL}(p^* | p_\theta)$ (forward KL) is mode-seeking: the model strongly penalizes assigning low probability to high-probability regions of $p^*$, but can ignore modes (assign them probability near zero).

This contrasts with reverse KL, $D_{KL}(p_\theta | p^*)$, which is mode-covering but rarely used in density estimation.

2. Overfitting

With powerful models (many parameters), MLE can memorize training data. For exact histogram density estimation, the estimate assigns probability $1/n$ to each unique training point—useless for generalization.

3. Infinite Likelihood Trap

For continuous data, likelihood can be made arbitrarily high by placing delta functions at training points. Regularization (priors, early stopping) is essential.

4. Blind to Perceptual Quality

MLE optimizes for log-probability, which doesn't directly correspond to human perception. A model can have high likelihood but generate blurry images, or low likelihood but great samples. This motivated alternative objectives like adversarial training (GANs).

A radical departure from classical density estimation: implicit density models define distributions through a sampling procedure rather than an explicit density function.

The Key Insight:

Instead of learning $p_\theta(x)$ explicitly, learn a function $G_\theta: z \mapsto x$ that transforms noise into samples:

$$z \sim \mathcal{N}(0, I), \quad x = G_\theta(z)$$

The distribution of $x$ is implicitly defined by $G_\theta$ and the noise distribution. We can sample from it easily (forward pass through $G$), but we cannot evaluate $p_\theta(x)$ for arbitrary $x$.

Why abandon explicit density?

1. Tractability: Computing normalizing constants and Jacobians is hard. Implicit models sidestep this entirely.

2. Flexibility: Without density constraints, $G_\theta$ can be any neural network.

3. Sample quality: Empirically, implicit models (especially GANs) produce sharper samples than likelihood-based models.

The Trade-off:

• Cannot compute $p_\theta(x)$—no likelihood evaluation
• Cannot easily detect out-of-distribution inputs
• Training requires alternative objectives (adversarial, moment matching)

Generative Adversarial Networks (GANs):

The most famous implicit model. A generator $G$ tries to produce realistic samples, while a discriminator $D$ tries to distinguish real from fake:

$$\min_G \max_D \mathbb{E}{x \sim p{data}}[\log D(x)] + \mathbb{E}_{z \sim p_z}[\log(1 - D(G(z)))]$$

At equilibrium (in theory), $G$ produces samples indistinguishable from real data, implying $p_G = p_{data}$.

Moment Matching:

Another implicit approach: match statistics (moments) of generated and real distributions:

$$\min_\theta \left| \mathbb{E}{x \sim p{data}}[\phi(x)] - \mathbb{E}{z \sim p_z}[\phi(G\theta(z))] \right|^2$$

where $\phi$ extracts features. This is the basis of Maximum Mean Discrepancy (MMD) training.

Hybrid Approaches:

Some models combine implicit and explicit elements:

• BiGAN/ALI: GAN + encoder, enabling approximate inference
• Flow-GAN: Normalizing flow with adversarial loss
• VAE-GAN: VAE with discriminator for sharper samples

Density estimation is the mathematical backbone of generative modeling. We've traversed from classical approaches to modern neural methods, illuminating how different techniques trade off expressiveness, tractability, and computational cost.

What's next:

Density estimation tells us how to evaluate the probability of data. But generative models must also sample—create new data instances from the learned distribution. In the next page, we'll explore sampling techniques: how to convert learned distributions into tangible new examples, from simple ancestral sampling to sophisticated MCMC methods.