What generative modeling is

A classifier answers a question about an object: is this image a cat? A generative model produces the object itself: a new image of a cat that was not in the training set. This note sets up what that means precisely, why it is hard, and how the rest of the course attacks it.

The problem: examples of an unknown distribution#

We treat a data set as independent draws from some distribution:

$$x^{(1)}, \dots, x^{(N)} \sim p_{\text{data}}(x).$$

For images, $x$ is a vector of pixel values; for text, a sequence of tokens; for audio, a waveform. The distribution $p_{\text{data}}$ is unknown: nobody can write down the probability density of “photographs of cats”. All we have is the sample.

The goal is a model $p_\theta(x)$ with parameters $\theta$ such that $p_\theta \approx p_{\text{data}}$, and from which we can draw new samples $x \sim p_\theta$. If the approximation is good, the new samples look like the data without copying it.

Figure

Generative modeling: from examples to a model and back to examples

So yes: generative modeling is learning to sample from a distribution we do not know, by approximating it from examples. Two questions follow. What does “approximate” mean, so that we can train? And how does a computer sample from a distribution at all? The rest of the note answers both.

What a generative model is for#

Sampling new objects is the most visible use, but not the only one:

• Generation. New images, text, audio, molecules that resemble the training data.
• Conditional generation. A model of $p(x \mid c)$ produces data that matches a condition $c$: an image for a text prompt, the missing part of a photograph, a high-resolution version of a small image, a translation of a sentence.
• Density estimation. If the model can evaluate $p_\theta(x)$, the density is one possible anomaly score: low density can flag unusual objects. The score depends on the representation (a density changes under a change of units or variables) and on how well the model fits, and a high density at a point is neither a high probability of that point nor a guarantee that it looks typical; the anomaly detection note discusses these limits. Good density models are also good compressors, since an optimal code spends about $-\log_2 p(x)$ bits on $x$.
• Representation learning. Models with a latent variable (the VAE) map data to a low-dimensional code whose directions often carry meaning, learned without labels.

What “approximate” means: a divergence to minimize#

To train $p_\theta$ we need a number that says how far it is from $p_{\text{data}}$. The standard choice is the Kullback–Leibler divergence (Entropy, cross-entropy and KL divergence):

$$D_{\mathrm{KL}}(p_{\text{data}} \,\|\, p_\theta) = \mathbb{E}_{x \sim p_{\text{data}}}\big[\log p_{\text{data}}(x)\big] - \mathbb{E}_{x \sim p_{\text{data}}}\big[\log p_\theta(x)\big].$$

The first term does not depend on $\theta$. Minimizing the divergence therefore means maximizing the expected log-likelihood of the data under the model, and with the sample in place of the expectation this is maximum likelihood:

$$\theta^\star = \arg\max_\theta \frac{1}{N}\sum_{i=1}^{N} \log p_\theta\big(x^{(i)}\big).$$

The catch is that this needs $p_\theta(x)$ as a number. For flexible models that is often impossible to compute, and the families in this course differ mainly in how they get around it: an exact formula, a lower bound, or a different objective altogether.

The course starts with this objective on the simplest case: a 2-D Gaussian fitted to Gaussian data by gradient descent on the KL divergence (Fitting distributions with KL).

How a computer samples at all#

A model is only generative if we can draw from it. Computers do not draw from arbitrary distributions directly. Everything is built from one primitive, uniform random numbers, by a chain of transformations. The whole section is one idea:

Start with a distribution that is easy to sample, and transform its samples into samples of a harder one.

The chain has five links, and every later model rests on them, so we go through them slowly. Each step gives the idea, the key formula in a frame, why it works, and the same step in a few lines of code. Readers who already know inverse-CDF and Box–Muller sampling can skip to why this is hard and the six families.

Step 1: uniform numbers#

A computer cannot flip a coin. A processor is deterministic: the same program on the same input does the same thing every time. Random numbers in programs come from two sources.

• Physical noise. The operating system collects events nobody can predict, such as the exact timing of keyboard, network and disk interrupts and electrical noise inside the processor, and mixes them into a pool that programs read (/dev/urandom on Linux and macOS, os.urandom in Python). These bits are truly unpredictable, but they are slow to produce. They are used for cryptographic keys and for picking a starting point, not for the millions of numbers a simulation needs.
• Pseudo-random generators. A short deterministic rule turns one number into the next. The generator keeps an internal state $s$, and each call updates it, $s_{n+1} = T(s_n)$, and returns a number computed from it. The whole sequence is determined by the starting state, which is set from a seed. A good generator passes large published batteries of statistical tests; no generator passes every conceivable test, since a test that knows the rule can always predict the next number.

A toy generator. The simplest rule is the linear congruential generator: multiply, add, keep the remainder (Knuth, 1997, §3.2.1). With 16 possible states,

$$s_{n+1} = (5\,s_n + 3) \bmod 16, \qquad u_n = s_n / 16 .$$

Start from the seed $s_0 = 7$. Then $5 \cdot 7 + 3 = 38 = 2 \cdot 16 + 6$, so $s_1 = 6$. Next, $5 \cdot 6 + 3 = 33 = 2 \cdot 16 + 1$, so $s_2 = 1$. Continuing the same way gives the states

$$7,\ 6,\ 1,\ 8,\ 11,\ 10,\ 5,\ 12,\ 15,\ 14,\ 9,\ 0,\ 3,\ 2,\ 13,\ 4,\ 7,\ \dots$$

and dividing by 16 gives the numbers $0.4375, 0.375, 0.0625, 0.5, 0.6875, \dots$ After 16 steps the generator is back at 7 and repeats. On the way it has visited every state exactly once, so over one cycle each of the 16 values $0, \tfrac1{16}, \dots, \tfrac{15}{16}$ appears exactly once: no value is preferred. That is what “uniform” means for a generator.

Figure

Where the "random" numbers come from

What real generators change. The toy has two visible defects. It has only 16 values and repeats after 16 steps, and each number is tied to the previous one by a simple rule, so consecutive pairs $(u_n, u_{n+1})$ fall on a few straight lines. Real generators keep the idea of a state updated by a fixed rule and remove both defects.

• An enormous state. PCG64, NumPy’s default, has a 128-bit state and repeats only after $2^{128}$ numbers (O’Neill, 2014). The Mersenne Twister, PyTorch’s generator on the CPU, has a 19,937-bit state and a period of $2^{19937} - 1$ (Matsumoto and Nishimura, 1998).
• A scrambling output step. The state is mixed by bit shifts, rotations and multiplications before it is returned, so no pattern shows between consecutive numbers. PCG64 passes BigCrush, the largest battery of the TestU01 library, and the Mersenne Twister fails only its two tests of linear complexity (L’Ecuyer and Simard, 2007).
• Counters on GPUs. Thousands of GPU threads need numbers at once, so PyTorch uses Philox there. It computes the $n$-th number directly as a scrambled function of $n$ and the seed, so threads do not have to wait for each other (Salmon et al., 2011).

From an integer to a number in $[0, 1)$. A generator returns 64-bit integers $k \in \{0, 1, \dots, 2^{64} - 1\}$. A double-precision number has a 53-bit mantissa, so libraries keep the top 53 bits of $k$ and scale them:

$$u = \big\lfloor k / 2^{11} \big\rfloor \cdot 2^{-53} .$$

The result is one of $2^{53}$ equally spaced values in $[0, 1)$, each equally likely. With seed 0, NumPy’s first raw output is $k = 11\,749\,869\,230\,777\,074\,271$, and the formula gives $u = 0.63696\ldots$, exactly the first number that np.random.default_rng(0).random() returns.

The seed. The seed sets the starting state; the library expands a small integer into the full state. The same seed, with the same generator, the same calls with the same arguments in the same order and the same build of the library, gives the same numbers, which is what makes an experiment repeatable. That is a narrow guarantee. NumPy keeps the raw streams of its generators stable across versions, but the samplers built on them, such as normal or multivariate_normal, may change between feature releases or builds (NumPy documentation), and a whole GPU training run adds sources of variation of its own (a training run in practice). Different seeds do not automatically give independent streams either: in the toy generator above every seed lands on the same cycle of 16 states, so another seed only shifts the same sequence. Large generators make overlaps very unlikely, and for parallel streams NumPy derives child seeds with SeedSequence.spawn (NumPy documentation).

In code. The library does all of the above in one call:

import numpy as np
rng = np.random.default_rng(0)       # seed 0 fixes the whole sequence
u = rng.random(100_000)              # 100,000 numbers from U(0, 1)
print(u[:3].round(4))                # [0.637  0.2698 0.041 ]
print((u < 0.25).mean())             # 0.25015: a quarter of them fall below 0.25

Step 2: any one-dimensional distribution, by inverting the CDF#

A discrete example first. Suppose we need an outcome $a$, $b$ or $c$ with probabilities $0.2$, $0.5$ and $0.3$. Cut the interval $[0, 1)$ into three pieces of exactly these lengths, $[0, 0.2)$, $[0.2, 0.7)$ and $[0.7, 1)$. Then draw $u$ and return the outcome whose piece contains it. A uniform number falls into a piece with probability equal to the piece’s length, so each outcome comes out with the probability we asked for. For $u = 0.637$, the first number above, the result is $b$.

Figure

Sampling a discrete distribution: cut the interval [0, 1)

The cut points $0.2$, $0.7$ and $1$ are the running totals of the probabilities, $\Pr(\text{outcome} \le a)$, $\Pr(\text{outcome} \le b)$ and $\Pr(\text{outcome} \le c)$. That is, they are the values of the cumulative distribution function $F$. Sampling is a lookup: find the first outcome whose running total exceeds $u$.

The continuous case is the same idea with infinitely many pieces. For a continuous random variable the CDF $F(x) = \Pr(X \le x)$ rises smoothly from 0 to 1: it says how much probability lies to the left of $x$. Assume it is continuous and strictly increasing, so that every level $u$ is reached at exactly one point. Here is the trick: draw a uniform $u$ and ask which value of $x$ has accumulated exactly that much probability. That value is the sample:

$$\boxed{\,X = F^{-1}(U) \sim F\,} \qquad \text{where } F^{-1}(u) \text{ is the } x \text{ with } F(x) = u .$$

This is inverse transform sampling (Devroye, 1986, §II.2; the functions of random variables note of the probability course treats it formally).

Why it works. Fix any number $x$ and compute the CDF of the returned value in three lines:

$$\begin{aligned} \Pr(X \le x) &= \Pr\big(F^{-1}(U) \le x\big) \\ &= \Pr\big(U \le F(x)\big) \\ &= F(x) . \end{aligned}$$

The first line is the definition of $X$. The second holds because $F$ is increasing: the point where the CDF reaches $u$ lies at or left of $x$ exactly when $u$ is at most $F(x)$. The third holds because a uniform number is at most $F(x)$ with probability equal to the length of $[0, F(x)]$. So $X$ has CDF $F$: exactly the distribution we wanted.

Worked example: the exponential distribution. Waiting times with rate $\lambda$ have the CDF $F(x) = 1 - e^{-\lambda x}$ for $x \ge 0$. Solving $F(x) = u$ for $x$ takes three steps:

$$\begin{aligned} 1 - e^{-\lambda x} &= u, \\ e^{-\lambda x} &= 1 - u, \\ -\lambda x &= \ln(1 - u), \end{aligned} \qquad \boxed{\,x = -\frac{\ln(1 - u)}{\lambda}\,}$$

With $\lambda = 1$:

$u$	0.1	0.5	0.9	0.99
$x = -\ln(1 - u)$	0.105	0.693	2.303	4.605

Half of the uniform numbers ($u < 0.5$) land in the short stretch $[0, 0.693)$, and the last 1% are spread over everything beyond 4.6. The shape of the exponential, many short waits and a long tail of long ones, appears by itself.

In code. One line turns uniform numbers into exponential ones (log1p(-u) computes $\ln(1 - u)$ accurately for small $u$):

import numpy as np
u = np.random.default_rng(0).random(100_000)
x = -np.log1p(-u)                                 # x = -ln(1 - u): exponential with rate 1
print(x.mean().round(3), np.median(x).round(3))   # 0.998 0.691 (theory: 1 and ln 2 = 0.693)

The normal distribution. For the standard normal, $F$ is $\Phi(x) = \tfrac12\big(1 + \operatorname{erf}(x/\sqrt2)\big)$. The recipe is the same, but $\Phi^{-1}$ has no formula in elementary functions. Libraries compute it with rational approximations accurate to about 16 digits, for example Wichura’s algorithm (Wichura, 1988). The figure shows why the recipe produces a bell: equal stretches of $u$ land on short stretches of $x$ where $\Phi$ is steep and on long ones where it is flat.

Figure

From uniform to normal through the inverse CDF

Step 3: the standard normal, by Box–Muller#

Inverting $\Phi$ numerically works, but it is slow. Box and Muller found an exact recipe that turns two independent uniforms $U_1, U_2$ into two independent standard normals (Box and Muller, 1958). First build a distance and an angle,

$$R = \sqrt{-2\ln U_1}, \qquad \varphi = 2\pi U_2,$$

then convert them to coordinates:

$$\boxed{\begin{aligned} Z_1 &= R\cos\varphi, \\ Z_2 &= R\sin\varphi . \end{aligned}}$$

The recipe looks arbitrary, but it is step 2 applied in polar coordinates. Think of $(Z_1, Z_2)$ as a random point in the plane, and describe it by its angle $\Theta$ and its distance $R$ from the origin.

1. The angle is uniform. Two independent standard normals have the joint density $\frac{1}{2\pi}e^{-z_1^2/2}e^{-z_2^2/2} = \frac{1}{2\pi}e^{-r^2/2}$, where $r^2 = z_1^2 + z_2^2$. It depends only on the distance from the origin, so the cloud of points is circularly symmetric. Every direction is equally likely, so $\Theta$ is uniform on $[0, 2\pi)$, and it does not depend on $R$. One uniform number gives the angle: $\Theta = 2\pi U_2$, the $\varphi$ above.

2. The squared distance is exponential. The probability that the point lands within distance $\sqrt t$ of the origin is the integral of the density over that disc. In polar coordinates the area element is $r\,dr\,d\theta$, and the integral over the angle gives $2\pi$, which cancels the $\frac{1}{2\pi}$:

   $$\Pr(R^2 \le t) = \int_0^{2\pi}\!\!\int_0^{\sqrt t} \frac{1}{2\pi}\, e^{-r^2/2}\, r\,dr\,d\theta = \Big[-e^{-r^2/2}\Big]_0^{\sqrt t} = 1 - e^{-t/2} .$$

   This is the exponential CDF from the worked example, with rate $\lambda = \tfrac12$.

3. Sample the squared distance by step 2. The worked example gives $R^2 = -\ln(1 - U_1)/\tfrac12 = -2\ln(1 - U_1)$. The number $1 - U_1$ is uniform too, so we may use $U_1$ itself: $R^2 = -2\ln U_1$ and $R = \sqrt{-2\ln U_1}$. (The distance $R$ itself has the Rayleigh distribution.) In exact arithmetic $U_1 = 0$ has probability zero, but a floating-point generator on $[0, 1)$ can return exactly 0, where the logarithm is $-\infty$; the code below therefore feeds the logarithm $1 - U_1$, which lies in $(0, 1]$.

4. Go back to coordinates. $Z_1 = R\cos\Theta$ and $Z_2 = R\sin\Theta$, which is the boxed formula.

Worked example. Take $u_1 = 0.3$ and $u_2 = 0.8$.

• Distance: $R^2 = -2\ln 0.3 = 2.408$, so $R = 1.552$. As a check, $\Pr(R^2 \le 2.408) = 1 - e^{-1.204} = 0.7 = 1 - u_1$.
• Angle: $\Theta = 2\pi \cdot 0.8$, that is $288^\circ$, with $\cos\Theta = 0.309$ and $\sin\Theta = -0.951$.
• Point: $z = (1.552 \cdot 0.309,\ 1.552 \cdot (-0.951)) = (0.480, -1.476)$.

Figure

Box–Muller step by step: two uniform numbers → a point of the normal cloud

In code. The recipe gives standard deviation 1 in both coordinates and no correlation between them. In practice we call the library’s own normal sampler, which uses a faster variant of the same principle (NumPy’s uses the ziggurat method):

import numpy as np
rng = np.random.default_rng(1)
u1, u2 = 1 - rng.random(100_000), rng.random(100_000)   # 1 - U is in (0, 1]: log never sees 0
r, phi = np.sqrt(-2 * np.log(u1)), 2 * np.pi * u2
z1, z2 = r * np.cos(phi), r * np.sin(phi)
print(z1.std().round(3), z2.std().round(3), np.corrcoef(z1, z2)[0, 1].round(3))   # 0.998 1.002 -0.002
print((np.abs(z1) < 1).mean().round(3))   # 0.685: the familiar 68% within one standard deviation

Step 4: any Gaussian, by an affine map#

Steps 1–3 give standard normal noise $Z \sim \mathcal N(\mathbf 0, \mathbf I)$: a round cloud centred at zero. A Gaussian with any mean $\boldsymbol\mu$ and covariance $\boldsymbol\Sigma$ comes from it by a shift and a linear map. Take any matrix $\mathbf A$ with $\mathbf A\mathbf A^{\top} = \boldsymbol\Sigma$; then

$$\boxed{\,X = \boldsymbol\mu + \mathbf A Z \sim \mathcal N(\boldsymbol\mu, \boldsymbol\Sigma)\,}$$

The mean shifts the cloud, and the matrix stretches, rotates and correlates its coordinates. The next note, Multivariate Gaussian distributions, shows how to choose $\mathbf A$ and why it works.

In code. The Cholesky factor is one valid choice of $\mathbf A$:

import numpy as np
mu = np.array([1.0, 2.0])
Sigma = np.array([[2.0, 0.8], [0.8, 1.0]])
A = np.linalg.cholesky(Sigma)                 # lower triangular, A @ A.T == Sigma
z = np.random.default_rng(2).standard_normal((100_000, 2))
x = mu + z @ A.T
print(A.round(3).ravel())                     # [1.414 0.    0.566 0.825]
print(x.mean(axis=0).round(2), np.cov(x.T).round(2).ravel())   # [1. 2.] [2.  0.8 0.8 1. ]

Step 5: anything else, by a learned nonlinear map#

Gaussians are only ellipsoids, while images, speech and text have distributions far more complicated than an ellipse. To reach the shape of real data we push the noise through a nonlinear map $g$:

$$\boxed{\,Z \sim \mathcal N(\mathbf 0, \mathbf I) \;\longrightarrow\; X = g_\theta(Z)\,}$$

The distribution of $X$ is called the pushforward of the noise distribution through $g$. In a generative model, $g = g_\theta$ is a neural network with parameters $\theta$, and training chooses $\theta$ so that the pushforward matches the data.

Figure

Sampling as a chain of transformations of noise

This is the central picture of the course: a generative model is simple noise plus a learned transformation. For a direct generator, sampling is easy by construction (draw noise, apply the network once); some families below pay more for each sample, for example a chain of many network evaluations. The hard part is the other direction: choosing $\theta$, which needs some way to compare the pushforward with the data.

When $g$ is invertible and smooth, the density of $X$ even has a formula, the change of variables

$$p_X(x) = p_Z\big(g^{-1}(x)\big)\,\left|\det \frac{\partial g^{-1}(x)}{\partial x}\right|,$$

which is the starting point of normalizing flows. For a general network $g_\theta$ there is no such formula, and that is where the families diverge.

Why this is hard#

Inverting a CDF solves sampling in one dimension, and a Gaussian has a formula in any dimension. Real data has neither.

• Dimension. A $256 \times 256$ colour image has $256 \cdot 256 \cdot 3 = 196{,}608$ numbers. A density on that space cannot be tabulated, and there is no tractable general inverse-sampling procedure. A joint CDF exists, but it has no single inverse; sampling coordinate by coordinate would need every conditional CDF, which is exactly what we do not know.
• Concentration. Natural images fill a tiny part of that space, close to a much lower-dimensional surface (a manifold): almost every random pixel array is noise. The model has to put nearly all its probability on that thin set.
• Normalization. A flexible function $f_\theta(x) \ge 0$ becomes a density only after dividing by $\int f_\theta(x)\,dx$, an integral over the whole space that is intractable for neural networks. Without it, the likelihood in the maximum-likelihood objective cannot be evaluated.

Six ways around it#

Each family in the course keeps the “noise plus transformation” picture and differs in what it gives up to make training possible.

Figure

Six families: noise becomes data; what differs is what can be computed

Family	How a sample is drawn	What is learned	What can be computed	Course notes
Variational autoencoder (VAE)	$z \sim \mathcal N(\mathbf 0, \mathbf I)$, then decode	decoder $p_\theta(x \mid z)$ and encoder $q_\phi(z \mid x)$	a lower bound on $\log p_\theta(x)$, the ELBO	VAE
Generative adversarial network (GAN)	$z \sim p(z)$, then $x = G_\theta(z)$ in one pass	generator $G_\theta(z)$, trained against a discriminator	no density; only samples	GAN, WGAN, scaling, StyleGAN, α-GAN, CycleGAN, StarGAN, beyond CycleGAN
Normalizing flow	$z \sim \mathcal N(\mathbf 0, \mathbf I)$, then a chain of invertible maps	invertible $f_\theta$ with a tractable Jacobian	the exact $\log p_\theta(x)$	Flows, architectures
Energy-based model (EBM)	start from noise, run a noisy Markov chain (Langevin) that explores $e^{-E_\theta}$; many steps, mixing can be slow	an energy $E_\theta(x)$, with $p_\theta(x) = e^{-E_\theta(x)}/Z_\theta$	$\log p_\theta(x)$ up to the unknown constant $\log Z_\theta$; samples by MCMC	MCMC samplers, EBMs
Diffusion model (DDPM)	start from noise $x_T$, denoise step by step to $x_0$	a network that undoes one small noising step	a lower bound on $\log p_\theta(x)$	Forward process, training, score matching, latent diffusion, drag editing
Autoregressive model	one value or token at a time, each given the previous ones	conditionals $p_\theta(x_i \mid x_{<i})$, each a softmax over values or tokens	the exact $\log p_\theta(x)$	PixelRNN and PixelCNN, transformers

The last note carries these tools from images to sound: a flow that turns text into speech, and an autoregressive transformer over the tokens of a neural audio codec that generates music.

Because not every model gives a likelihood, and a likelihood does not always reflect how good samples look, the course also covers metrics computed from samples alone: Inception Score, FID, LPIPS, precision and recall and human evaluation.

The plan of the course#

1. Foundations. Gaussians and sampling, entropy and KL divergence, fitting a distribution by minimizing KL, the training background (batch normalization, gradients, convolutions) and autoencoders.
2. Latent-variable models. The VAE: a latent code, a lower bound on the likelihood and the reparameterization trick; β-VAE and the discrete codes of VQ-VAE.
3. Adversarial models. GANs, the Jensen–Shannon divergence they implicitly minimize, and Wasserstein GANs; scaling GANs up (ProGAN, SAGAN, BigGAN) and controlling them (StyleGAN); autoencoders trained with adversarial losses (α-GAN); image-to-image translation without paired examples (CycleGAN), across many domains (StarGAN) and beyond (MUNIT, CUT, SPADE, super-resolution).
4. Evaluation. IS, FID and LPIPS, precision and recall, and human evaluation.
5. Exact likelihood. Normalizing flows and the architectures that keep the Jacobian cheap.
6. Energy-based models. Samplers for densities known up to a constant (MALA, Hamiltonian Monte Carlo, NUTS), and energy-based models trained with contrastive divergence.
7. Diffusion. The forward noising process and the training objective of DDPM, score matching and Langevin dynamics, latent diffusion with Stable Diffusion, and point-based editing with DragGAN and DragDiffusion.
8. Autoregressive models. The chain rule as a model, PixelRNN and PixelCNN, and transformers over pixels and image tokens.
9. Audio. Text-to-speech with a flow and music with an autoregressive transformer over neural-codec tokens.

What you should remember#

• Data is a sample from an unknown $p_{\text{data}}$; a generative model is a $p_\theta$ we can sample from, fitted so that $p_\theta \approx p_{\text{data}}$.
• Minimizing $D_{\mathrm{KL}}(p_{\text{data}} \,\|\, p_\theta)$ is maximum likelihood, which needs $p_\theta(x)$ as a number.
• All sampling starts from uniform numbers: inverse CDF in one dimension, Box–Muller for normals, an affine map for any Gaussian, a learned nonlinear map for everything else.
• Every model in the course is noise plus a learned transformation; the families differ in what they can compute: a bound (VAE, diffusion), nothing (GAN), the likelihood up to a constant (energy-based models) or the exact likelihood (flows, autoregressive models).

Next#

The first step of the ladder that the course builds on is the Gaussian: how its covariance shapes the data and how an affine map turns standard noise into any Gaussian. That is Multivariate Gaussian distributions.

References#

• I. Goodfellow, Y. Bengio, A. Courville. Deep Learning. MIT Press, 2016. Ch. 20: deep generative models.
• K. P. Murphy. Probabilistic Machine Learning: Advanced Topics. MIT Press, 2023. Part IV: generation, with an overview of the model families.
• C. M. Bishop. Pattern Recognition and Machine Learning. Springer, 2006. §11.1: basic sampling algorithms, including the inverse transform and Box–Muller.
• L. Devroye. Non-Uniform Random Variate Generation. Springer, 1986.
• G. E. P. Box, M. E. Muller. A note on the generation of random normal deviates. Annals of Mathematical Statistics, 29(2), 1958.
• D. E. Knuth. The Art of Computer Programming, vol. 2: Seminumerical Algorithms, 3rd ed. Addison-Wesley, 1997. §3.2.1: linear congruential generators and when they visit every state.
• P. L’Ecuyer, R. Simard. TestU01: a C library for empirical testing of random number generators. ACM Transactions on Mathematical Software, 33(4), 2007.
• M. J. Wichura. Algorithm AS 241: the percentage points of the normal distribution. Applied Statistics, 37(3), 1988.
• M. Matsumoto, T. Nishimura. Mersenne Twister: a 623-dimensionally equidistributed uniform pseudo-random number generator. ACM Transactions on Modeling and Computer Simulation, 8(1), 1998.
• M. E. O’Neill. PCG: a family of simple fast space-efficient statistically good algorithms for random number generation. Technical report HMC-CS-2014-0905, Harvey Mudd College, 2014.
• J. K. Salmon, M. A. Moraes, R. O. Dror, D. E. Shaw. Parallel random numbers: as easy as 1, 2, 3. SC ‘11, 2011. The Philox generator.
• NumPy developers. Random sampling: compatibility policy and parallel random number generation. NumPy documentation. What a seed does and does not reproduce; SeedSequence.spawn for parallel streams.