In this course, you will need to be able to read and write probability models in a particular mathematical language. We will interpret that language in terms of simulation. The fundamental issue is that there are many different ways you could simulate random values, and we need a consistent way of describing the sorts of values a simulator should give. We are going to use the language of probability to do this.
We think of a simulator as a procedure that produces a random output. In general, we will be simulating datasets in pieces that are assembled in a way that either corresponds to a scientific model or which we otherwise think will be useful. We therefore need a way to refer to these pieces in the abstract, without referring to the result of any particular run.
A piece of a simulator that produces a real number as output will correspond to a random variable, which we will denote by a capital letter, for example \(Y\). We can think of a random variable \(Y\) as a placeholder for the value produced by that step before we have actually run the simulator. When we do run the simulator and get a particular value, we denote that realized value by the corresponding lowercase letter, e.g. \(y\). So we often see equalities that look like \(Y = y\), which just means that the piece of the simulator corresponding to the random variable \(Y\) actually produced the value \(y\).
The next issue is describing what sorts of values the simulator actually generates and how frequently it generates them. If \(Y\) is random variable, then the way we define what fraction of the time the outcome of \(Y\) should fall in various sets is given by a probability law.
For example, the law tells us the probability that \(Y\) falls in an interval:
\[ \Pr(a<Y\le b). \]
It does this for every pair of real numbers \(a<b\). Thus, if we know the law of \(Y\), it tells us things like:
\[ \Pr(0<Y\le 1), \qquad \Pr(-2<Y\le 3), \qquad \Pr(Y>10). \]
The description above leaves out a lot of technical details. For instance, for a probability law to be useful, the probabilities assigned to various sets must fit together in a coherent way. If you take probability theory, you will learn that it is not actually possible to create a probability law that will assign coherent probabilities to all possible sets. The construction of coherent probability laws is interesting mathematically, but in this course, we will usually be assuming that someone else has already defined a probability law for us that we can use, and so we don’t need to worry too much about it.
It does turn out that a law that gives coherent probabilities to all intervals is enough to determine the probability of any reasonable set of real numbers that we will encounter. We will therefore use the following as our working description:
A probability law for a real-valued quantity specifies the probability that the quantity falls in each interval.
The notation we use to indicate that \(Y\) follows the probability law \(P\) is \[ Y \sim P, \] which we read as “\(Y\) is distributed according to \(P\).”
Thus, if we say that we have a simulator or a piece of a simulator corresponding to a random variable \(Y\) that follows a probability law \(P\), what we mean is that the probability that the simulator/piece of the simulator generates a value that falls in any interval \((a,b]\) is given by \(P((a,b])\).
Probability laws are much more intuitive for random variables that take either a finite or countably infinite number of values, so we briefly describe them here. Suppose that \(Y\) can take values only in a finite or countable set, such as
\[ \{0,1\}, \qquad \{0,1,2,\ldots\}, \qquad \text{or} \qquad \{1,2,\ldots,K\}. \]
In this case, specifying the law is equivalent to specifying the probability of every individual possible outcome. If the possible values are \(y_1,y_2,\ldots\), write
\[ p_j=\Pr(Y=y_j). \]
These probabilities must satisfy
\[ p_j\ge 0 \qquad\text{and}\qquad \sum_j p_j=1. \]
Once the individual probabilities are known, the probability of any set of possible outcomes is found by adding:
\[ \Pr(Y\in A)=\sum_{j:y_j\in A}p_j. \]
For example, suppose \(Y\) takes values \(0\), \(1\), and \(2\), with
\[ \Pr(Y=0)=0.2, \qquad \Pr(Y=1)=0.5, \qquad \Pr(Y=2)=0.3. \]
Then
\[ \Pr(0<Y\le 2) =\Pr(Y=1)+\Pr(Y=2) =0.8. \]
Thus, the point probabilities determine all the interval probabilities. For a discrete distribution, it is therefore usually simpler to describe the law using those individual probabilities rather than all of its interval probabilities. This will come up a bit later when we describe probability mass functions, but for now we are still just assuming that someone else has written down a valid law that we can use.
So far, we have used a random variable such as \(Y\) to represent one real-valued piece of a simulator. A dataset will usually contain many values, however, and those values may be related to one another. To describe a simulator for an entire dataset, we need notation that specifies both the probability behavior of the individual pieces and how those pieces are related.
Suppose, for example, that a simulator produces two values, represented by the random variables \(Y_1\) and \(Y_2\). We can collect them into a single random vector (which we usually denote by a bold-face uppercase letter),
\[ \mathbf{Y}=(Y_1,Y_2). \]
A run of the simulator produces a particular pair of values,
\[ \mathbf{y}=(y_1,y_2). \]
Just as a single random variable has a probability law, the random vector \(\mathbf{Y}\) has a joint probability law.
The joint law of \(Y_1\) and \(Y_2\) describes how the two values vary together. For example, it tells us probabilities such as
\[ \Pr(a_1<Y_1\le b_1,\;a_2<Y_2\le b_2). \]
This is the probability that the first output falls in \((a_1,b_1]\) and, in the same run of the simulator, the second output falls in \((a_2,b_2]\).
Knowing these probabilities for all pairs of intervals is enough to determine the joint law of two real-valued random variables. We can write
\[ (Y_1,Y_2)\sim P \]
to indicate that the pair has joint law (P). More generally, a dataset represented by
\[ \mathbf{Y}=(Y_1,\ldots,Y_n) \]
has a joint law that specifies probabilities for simultaneous statements about all \(n\) of its components.
The law of one component considered by itself is called its marginal law. For example, the joint law of \((Y_1,Y_2)\) determines the marginal laws of \(Y_1\) and \(Y_2\). However, the reverse is not true: knowing the two marginal laws does not generally determine the joint law.
For example, suppose that \(Y_1\) and \(Y_2\) both take values in \(\{0,1\}\), with
\[ \Pr(Y_1=1)=\Pr(Y_2=1)=\frac12. \]
One possible joint law assigns equal probability to all four possible pairs:
\[ \begin{aligned} \Pr(Y_1=0,Y_2=0)&=\frac14, & \Pr(Y_1=0,Y_2=1)&=\frac14,\\ \Pr(Y_1=1,Y_2=0)&=\frac14, & \Pr(Y_1=1,Y_2=1)&=\frac14. \end{aligned} \]
Another possible joint law is
\[ \Pr(Y_1=0,Y_2=0)=\frac12, \qquad \Pr(Y_1=1,Y_2=1)=\frac12, \]
with the other two pairs having probability zero. Under this law, \(Y_1\) and \(Y_2\) are always equal.
Both joint laws give
\[ \Pr(Y_1=1)=\Pr(Y_2=1)=\frac12, \]
but they describe very different simulators. Specifying the law of every component separately therefore does not specify the law of the complete dataset. We must also say how the components are related.
In principle, we could do this by writing down the full joint law directly. For a large dataset, however, this quickly becomes unwieldy. A collection of \(n\) binary random variables alone has \(2^n\) possible joint outcomes. Probability models therefore usually specify joint laws using simpler building blocks. Two particularly important building blocks are independence and conditional distributions.
The simplest way to assemble several random variables is to generate them independently. Informally, this means that each part of the simulator is run without using the values produced by the other parts.
Random variables \(Y_1,\ldots,Y_n\) are independent if probabilities of simultaneous events can be obtained by multiplying their separate probabilities. That is,
\[ \Pr(Y_1\in A_1,\ldots,Y_n\in A_n) = \prod_{i=1}^n \Pr(Y_i\in A_i) \]
for every collection of appropriate sets \(A_1,\ldots,A_n\).
For example, if \(Y_1\) and \(Y_2\) are independent and
\[ \Pr(Y_1=1)=\Pr(Y_2=1)=\frac12, \]
then
\[ \Pr(Y_1=1,Y_2=1) = \Pr(Y_1=1)\Pr(Y_2=1) = \frac14. \]
Thus, independence combines the separate laws into a particular joint law.
If \(Y_1,\ldots,Y_n\) are independent and all follow the same probability law \(P\), we write
\[ Y_1,\ldots,Y_n\overset{\mathrm{iid}}{\sim}P. \]
The abbreviation “iid” means independent and identically distributed. As a simulation instruction, this says:
Run the same simulator \(n\) times independently, and call the resulting outputs \(Y_1,\ldots,Y_n\).
The two parts of the statement do different jobs. “Identically distributed” says that every run follows the same law \(P\). “Independent” says how the different runs are related. It is possible for random variables to have the same distribution without being independent, as in the example where \(Y_1=Y_2\) on every run.
Independence is useful when you want the pieces of the simulator to not affect one another. In many models, however, a value generated at one step changes what should happen at a later step. We describe this using a conditional probability law.
Suppose that we first generate \(Y_1\), and then generate \(Y_2\) using a law that depends on the value produced for \(Y_1\). We can write the simulation instructions as
\[ Y_1\sim P_1, \]
followed by
\[ Y_2\mid Y_1=y_1 \sim P_2(\mathord{\cdot}\mid y_1). \]
The second line is read as “conditional on \(Y_1=y_1\), \(Y_2\) follows the law \(P_2(\mathord{\cdot}\mid y_1)\).” The dot marks the place where the possible values of \(Y_2\) would go.
These two lines describe the following simulation procedure:
For a discrete example, suppose that \(Y_1,Y_2\in\{0,1\}\) and
\[ \Pr(Y_1=1)=0.5. \]
After generating \(Y_1\), suppose we generate \(Y_2\) according to
\[ \Pr(Y_2=1\mid Y_1=0)=0.2 \]
and
\[ \Pr(Y_2=1\mid Y_1=1)=0.8. \]
The probability of generating \(Y_1=0\) and then \(Y_2=1\) is
\[ \begin{aligned} \Pr(Y_1=0,Y_2=1) &= \Pr(Y_1=0) \Pr(Y_2=1\mid Y_1=0)\\ &= 0.5(0.2)\\ &= 0.1. \end{aligned} \]
Similarly,
\[ \Pr(Y_1=1,Y_2=1) = 0.5(0.8) = 0.4. \]
The marginal law for the second output is obtained by considering both ways in which \(Y_2=1\) can occur:
\[ \begin{aligned} \Pr(Y_2=1) &= \Pr(Y_1=0,Y_2=1) + \Pr(Y_1=1,Y_2=1)\\ &= 0.1+0.4\\ &= 0.5. \end{aligned} \]
Thus, both \(Y_1\) and \(Y_2\) have marginal probability \(0.5\) of equaling 1, but they are not independent. Once we know the value of \(Y_1\), our probability for \(Y_2\) changes.
Conditional laws let us construct much larger simulators in stages. For example, we might specify
\[ Y_1\sim P_1, \]
\[ Y_2\mid Y_1=y_1 \sim P_2(\mathord{\cdot}\mid y_1), \]
and then
\[ Y_3\mid Y_1=y_1,Y_2=y_2 \sim P_3(\mathord{\cdot}\mid y_1,y_2). \]
Continuing in this way gives a simulation procedure for the entire dataset. These sequential instructions determine a joint law for
\[ (Y_1,\ldots,Y_n), \]
even if we never write that joint law down directly.
Independence can be viewed as a special case of conditional simulation. If the conditional law of \(Y_2\) is the same for every possible value of \(Y_1\), so that
\[ Y_2\mid Y_1=y_1\sim P_2 \]
regardless of \(y_1\), then learning the result of the first step does not change the law used at the second step. In that case, \(Y_1\) and \(Y_2\) are independent.
This gives us two basic ways to assemble a simulator:
In both cases, the complete collection of simulation instructions specifies a joint probability law for the resulting dataset.
We now have the language needed to specify a probability law for an entire dataset. We can describe the joint law directly, or we can construct it from independent and conditionally distributed pieces.
A single joint law specifies one particular probability distribution for the dataset. In a statistical problem, however, we usually do not know exactly which law generated the observed data. Instead, we identify a collection of laws that we are willing to consider as reasonable possibilities. This collection is called a probability model or statistical model.
We will often write a probability model as
\[ \mathcal{P} = \{P_\theta:\theta\in\Theta\}. \]
Here:
If \(\mathbf{Y}\) is a random variable representing the complete dataset, we write
\[ \mathbf{Y}\sim P_\theta \]
to say that, when the parameter is \(\theta\), the we think the dataset has come from the joint law \(P_\theta\).
From the perspective of simulation, each value of \(\theta\) gives us a different simulator for the dataset. We can read
\[ \mathbf{Y}\sim P_\theta \]
as the instruction:
Set the simulator’s parameter to \(\theta\), run the simulator, and call the resulting dataset \(\mathbf{Y}\).
As \(\theta\) varies over \(\Theta\), the behavior of the simulator changes. The probability model is the collection of all the probability laws that can be produced in this way. As we discussed in the first class, the task of estimating \(\theta\) boils down to the task of finding \(\theta\) that makes simulated datasets look like the real data.
Suppose that a dataset consists of \(n\) binary observations,
\[ \mathbf{Y}=(Y_1,\ldots,Y_n), \]
where \(Y_i=1\) represents a success and \(Y_i=0\) represents a failure. One possible model is
\[ Y_1,\ldots,Y_n \overset{\mathrm{iid}}{\sim} \operatorname{Bernoulli}(p), \qquad p\in[0,1]. \]
The Bernoulli law with parameter \(p\) is the probability law on \(\{0,1\}\) defined by
\[ \Pr(Y_i=1)=p \qquad\text{and}\qquad \Pr(Y_i=0)=1-p. \]
The model makes several distinct claims:
For any fixed value of \(p\), these claims determine a joint law for the complete dataset. For example, if a particular binary dataset \(\mathbf y\) contains \(s\) successes and \(n-s\) failures, independence gives
\[ \Pr_p(\mathbf{Y}=\mathbf y) = p^s(1-p)^{n-s}. \]
Changing \(p\) changes the joint law. A simulator with \(p=0.1\) will usually produce datasets containing relatively few successes, whereas a simulator with \(p=0.9\) will usually produce datasets containing many successes.
Thus, this probability model is the collection
\[ \mathcal P = \left\{ P_p:p\in[0,1] \right\}, \]
where \(P_p\) is the joint law of \(n\) independent Bernoulli observations with success probability \(p\).
Frequently used probability laws are given names. For example,
\[ Y\sim\operatorname{Binomial}(n,p) \]
says that \(Y\) has a binomial distribution. One generative construction of this law is:
In R or python, one draw can be generated with
rbinom(1, size = n, prob = p)or
import numpy as np
np.random.binomial(n,p)The R/python functions need not internally simulate all \(n\) Bernoulli trials. They only need to return values with the correct binomial law. Whether the individual trials are meaningful parts of the scientific story depends on the application.
Some distributions that we will use are:
| Distribution | Possible outputs | Meaning of the parameters |
|---|---|---|
| \(\operatorname{Bernoulli}(p)\) | \(0\) or \(1\) | success probability \(p\) |
| \(\operatorname{Binomial}(n,p)\) | \(0,1,\ldots,n\) | number of trials \(n\), success probability \(p\) |
| \(\operatorname{Poisson}(\lambda)\) | \(0,1,2,\ldots\) | expected event count \(\lambda\) |
| \(\operatorname{Categorical}(\pi)\) | \(1,\ldots,K\) | category probabilities \(\pi_1,\ldots,\pi_K\) |
| \(\operatorname{Normal}(\mu,\sigma^2)\) | any real number | center \(\mu\), variance \(\sigma^2\) |
To use a distribution as part of a simulator, you initially need to know:
Its probability mass function or density will become useful later, but it is not required merely to simulate from the distribution.
\[ E(Y_i\mid X_i=x_i)=\beta_0+\beta_1x_i \]
specifies the conditional mean but not the variation around it or dependence among observations.
\[ Y_i\sim P_\theta \]
for every \(i\) does not specify whether the observations are independent.
A simulator for infections is incomplete if the dataset contains reported cases and reporting is imperfect.
Being able to predict each variable from all the others does not necessarily provide a coherent way to generate all variables together.
| Notation | Simulation interpretation |
|---|---|
| \(Y\sim P\) | generate \(Y\) with probability law \(P\) |
| \(Y\sim P_\theta\) | supply \(\theta\), then generate \(Y\) with law \(P_\theta\) |
| \((X,Y)\sim P_{X,Y}\) | jointly generate and return the pair \((X,Y)\) |
| \(Y_i\overset{\mathrm{ind}}{\sim}P_i\) | run the specified simulators independently and return all the outputs |
| \(Y_i\overset{\mathrm{iid}}{\sim}P_\theta\) | generate independent observations with the same distribution |
| \(Y\mid X=x;\theta\sim P_\theta(\cdot\mid x)\) | supply \(x\) and \(\theta\), then generate \(Y\) from a distribution depending on \(X\) and \(\theta\) |
| \(X\sim P_X;\ Y\mid X=x\sim P(\cdot\mid x)\) | generate \(X\), then generate \(Y\) using the realized \(X\) |
| \(Z\sim P_\theta^Z;\ Y\mid Z=z\sim P_\theta^Y(\cdot\mid z)\) | generate a latent quantity, then generate the observation |
Before treating a model as complete, verify that you can answer all of the following:
The last question is the most important. Probability notation is useful because it is a compact language for describing probability laws and their relationships. In this course, the notation has not done its job until you can translate it into a coherent simulator.
For each model, we will ask the same questions:
Most statistics students will have encountered linear regression in the form
\[ Y_i=\beta_0+\beta_1x_i+\varepsilon_i, \qquad \varepsilon_i\overset{\mathrm{iid}}{\sim}N(0,\sigma^2). \]
This is already a complete probability model for the responses conditional on the covariate values. It can therefore be read directly as a simulation procedure: hold the \(x_i\)’s fixed, generate independent normal errors, and add those errors to the regression line. This familiar example shows how standard statistical notation specifies a simulator, what roles the different parts of a model play, and why we may model the responses conditional on supplied covariates without specifying a probability law for the covariates themselves.
Suppose that \(x_1,\ldots,x_n\) are supplied covariate values. A simple linear regression model is
\[ Y_i\mid x_i;\beta_0,\beta_1,\sigma \mathrel{\overset{\mathrm{ind}}{\sim}} \operatorname{Normal}(\beta_0+\beta_1x_i,\sigma^2), \qquad i=1,\ldots,n. \]
An equivalent noise representation is
\[ Y_i=\beta_0+\beta_1x_i+\epsilon_i, \qquad \epsilon_i\mathrel{\overset{\mathrm{iid}}{\sim}} \operatorname{Normal}(0,\sigma^2). \]
The first representation emphasizes the probability law of each response. The second emphasizes a deterministic signal plus random deviation. They define the same conditional joint law for \(Y=(Y_1,\ldots,Y_n)\).
| Piece | Role in the simulator | Substantive or statistical significance |
|---|---|---|
| \(x_i\) | Supplied input for observation \(i\) | The model is conditional on these values; it does not explain how they arose. |
| \(\beta_0\) | Shifts the center of every response | Conditional mean at \(x=0\), if that value is meaningful. |
| \(\beta_1\) | Changes the center with \(x_i\) | Change in the conditional mean associated with a one-unit change in \(x\). |
| \(\sigma\) | Controls the spread around the line | Typical scale of deviations not explained by the linear mean. |
| Normal law | Specifies the shape of the deviations | Adds much more than the claim that the deviations have mean zero. |
| Independence | Specifies the joint relationship among responses | Rules out residual dependence after conditioning on the supplied \(x_i\). |
| \(Y_i\) | Generated and observed quantity | One response in the simulated dataset. |
This model is generative for \(Y\) conditional on the supplied design \(x\). If the scientific question also concerns how the \(X_i\) arise, we need to add a model for them:
\[ X_i\sim P_X, \qquad Y_i\mid X_i=x_i\sim \operatorname{Normal}(\beta_0+\beta_1x_i,\sigma^2). \]
Whether that addition is necessary depends on what one complete replicate of the study is supposed to include.
This is a nice example from the Bayesian Workflow book (Chapter 25). We will see here that two models can generate the same kind of observed data while making very different commitments about the process that produced them. An empirical model may provide a flexible description of how success varies with distance. A mechanistic model uses a proposed physical explanation to determine the shape of that relationship. Mechanistic interpretation can make a model more informative, but it also creates additional substantive ways for the model to be wrong.
Suppose putts are attempted at distances \(x_1,\ldots,x_J\). At distance \(x_j\), there are \(n_j\) attempts and the observed count of successful putts is \(Y_j\).
Both models begin with the observation stage
\[ Y_j\mid n_j,p_j \mathrel{\overset{\mathrm{ind}}{\sim}} \operatorname{Binomial}(n_j,p_j), \qquad j=1,\ldots,J. \]
This says that the distances and numbers of attempts are supplied, success has probability \(p_j\) at distance \(x_j\), and the attempts within each distance are treated as independent and exchangeable. The remaining question is how to specify \(p_j\).
One option is
\[ \operatorname{logit}(p_j)=a+bx_j, \]
or equivalently
\[ p_j=\operatorname{logit}^{-1}(a+bx_j) =\frac{1}{1+\exp\{-(a+bx_j)\}}, \] where \(\operatorname{logit}(p) = \log\left(\frac{p}{1-p}\right), 0<p<1\) and \(\operatorname{logit}^{-1}(\eta) = \frac{1}{1+\exp(-\eta)}, \eta\in\mathbb{R}.\)
Together with the binomial observation model, this is a complete conditional generator for the success counts.
| Piece | Role | Significance |
|---|---|---|
| \(x_j,n_j\) | Supplied design quantities | Distance and number of attempted putts. |
| \(a\) | Controls the height and location of the curve | Log-odds of success at distance zero; usually an extrapolative intercept rather than a direct physical quantity. |
| \(b\) | Controls how rapidly success changes with distance | An empirical slope on the log-odds scale. |
| \(p_j\) | Deterministic intermediate quantity | Success probability implied by the curve at distance \(x_j\). |
| \(Y_j\) | Generated observation | Number of successes among \(n_j\) attempts. |
The model describes a smooth decreasing curve when \(b<0\). It does not say why the curve should be logistic, nor does it describe the geometry or motor errors involved in putting.
Let \(R\) be the radius of the hole and \(r\) the radius of the ball. In a simplified geometric model, a putt from distance \(x_j\) succeeds if its angular error is small enough that the center of the ball passes within \(R-r\) of the center line. The angular tolerance is
\[ \delta(x_j)=\sin^{-1}\left(\frac{R-r}{x_j}\right). \]
For putt \(\ell\) at distance \(x_j\), suppose
\[ A_{j\ell}\mid\sigma \mathrel{\overset{\mathrm{iid}}{\sim}} \operatorname{Normal}(0,\sigma^2), \]
and define the success indicator
\[ S_{j\ell}=\mathbf{1}\{|A_{j\ell}|<\delta(x_j)\}. \]
The recorded count is
\[ Y_j=\sum_{\ell=1}^{n_j}S_{j\ell}. \]
If we record only the number of successes and do not keep the individual angular errors, the same simulator implies
\[ Y_j\mid n_j,x_j,\sigma \sim \operatorname{Binomial}(n_j,p_j), \]
where
\[ p_j =\Pr\{|A_{j\ell}|<\delta(x_j)\}, \qquad A_{j\ell}\sim\operatorname{Normal}(0,\sigma^2). \]
We do not yet need a formula for evaluating this probability. The simulator can generate each \(A_{j\ell}\) and check whether it falls within the tolerance.
| Piece | Role | Significance |
|---|---|---|
| \(R-r\) | Fixed geometric clearance | Determines how far the ball’s center may miss the center line. |
| \(\delta(x_j)\) | Distance-dependent angular tolerance | Converts the geometry into a success criterion. |
| \(\sigma\) | Spread of angular errors | A potentially interpretable measure of directional accuracy. |
| \(A_{j\ell}\) | Latent variable for one putt | Unobserved aiming error generated anew on every attempt. |
| \(S_{j\ell}\) | Deterministic function of geometry and error | Indicates whether a particular putt succeeds. |
| \(Y_j\) | Aggregated observation | What remains after individual attempts are counted. |
The logistic model says that \((a,b)\) describe the success curve. The angular-error model says that the curve arises from geometry plus a distance-invariant normal distribution of aiming errors. Its parameter \(\sigma\) has a clearer physical interpretation, and the geometry determines the nonlinear shape of the distance effect.
That interpretation depends on assumptions that may be false. For example:
The mechanistic model is not automatically better because it is more mechanistic. It is useful when its additional claims are scientifically defensible and when it actually fits the data.
In this example, we see how a set of local conditional rules can define a complete joint generator when the rules are arranged in a noncircular order. A simple Markov language model also introduces autoregressive generation: a newly generated output becomes an input to the next simulation step.
Let the vocabulary contain \(K\) tokens, including a special end-of-sequence token. Let \(W_t\) be the token at position \(t\). A first-order Markov model is
\[ W_1\sim\operatorname{Categorical}(\pi), \]
\[ W_t\mid W_{t-1}=j \sim \operatorname{Categorical}(A_{j,\cdot}), \qquad t=2,3,\ldots. \]
Here \(\pi\) is a vector of initial-token probabilities and \(A\) is a \(K\times K\) transition matrix. Its entry
\[ A_{jk}=\Pr(W_t=k\mid W_{t-1}=j) \]
is the probability of token \(k\) following token \(j\). Every row of \(A\) is nonnegative and sums to one.
For a fixed sequence length \(T\), the ordered conditional model defines the joint law
\[ \Pr(W_{1:T}=w_{1:T}) =\pi_{w_1}\prod_{t=2}^T A_{w_{t-1},w_t}. \]
If the end token is reached with probability one, the same model defines a distribution over finite variable-length sequences. Without that condition, the law may also assign positive probability to an infinite sequence.
| Piece | Role | Significance |
|---|---|---|
| \(\pi\) | Parameter supplied to the simulator | Controls how sequences begin. |
| \(A\) | Parameter shared over all positions | Records local token-to-token transition probabilities. |
| \(W_t\) | Generated token | Both an output at step \(t\) and an input at step \(t+1\). |
| Markov assumption | Conditional-independence claim | Once \(W_{t-1}\) is known, earlier tokens do not affect \(W_t\). |
| End token | Generated stopping event | Turns the model into a distribution over sequence lengths as well as contents. |
This model is plainly unrealistic as a model of language: after conditioning on the immediately preceding token, it forgets the entire earlier context. Nevertheless, the ordering makes it a coherent generator. It also makes clear that joint does not mean “generated all at once.” A joint sample can be built sequentially from ordered conditional laws.
In this example, we see that a conditional prediction for every variable is not necessarily a joint probability model. When every conditional rule uses quantities that have not yet been generated, there may be no coherent starting point. More fundamentally, an arbitrary collection of proposed full conditional distributions may be mathematically incompatible: there may be no joint law having all of them as its conditionals.
Suppose \(Y_{jk}\) indicates whether student \(j\) answered exam question \(k\) correctly. Imagine fitting a separate model for each question using the student’s answers to all the other questions:
\[ Y_{jk}\mid Y_{j,-k} \sim \operatorname{Bernoulli} \left[ \operatorname{logit}^{-1} \left\{a_k+b_k s(Y_{j,-k})\right\} \right]. \]
Here:
This can be useful for predicting a missing response when all the other responses have been observed. It does not by itself tell us how to generate a new student’s complete response vector \((Y_{j1},\ldots,Y_{jK})\). To generate the first answer, the rule asks for all the other answers; each of those rules asks for the first answer in return.
There are two separate issues:
One might initialize all the answers arbitrarily and repeatedly update one answer at a time. That creates a Markov-chain algorithm, but it does not automatically solve the problem. The update order, initialization, and number of sweeps are then additional parts of the generator. If the chain has a stationary distribution, that distribution need not have the proposed regressions as its full conditionals unless the compatibility conditions hold.
Some specially constrained collections of binary conditional models do correspond to a joint model, as in an Ising or autologistic model. The lesson is not that circular conditional specifications can never be repaired. It is that they do not define a joint generator merely because a predictive model has been written for every component.
Here, we see a generative model that is not supposed to correspond to a mechanistic model of how the data were generated, but which can nonetheless describe useful structure in the data.
Suppose a corpus contains documents \(d=1,\ldots,D\). Document \(d\) contains \(N_d\) word tokens. Let there be \(K\) topics and a vocabulary of size \(V\).
The Dirichlet distribution used below generates a vector of nonnegative values that sum to one. Its parameter \(\alpha\) controls which kinds of topic-proportion vectors are common.
For each document, generate topic proportions
\[ \Theta_d\mid\alpha \sim \operatorname{Dirichlet}(\alpha). \]
For each token position \(n=1,\ldots,N_d\), generate a latent topic and then a word:
\[ Z_{dn}\mid\Theta_d \sim \operatorname{Categorical}(\Theta_d), \]
\[ W_{dn}\mid Z_{dn}=k;\beta_{1:K} \sim \operatorname{Categorical}(\beta_k). \]
The vector \(\beta_k=(\beta_{k1},\ldots,\beta_{kV})\) gives the word probabilities for topic \(k\). Each \(\beta_k\) is nonnegative and sums to one. In this version of the model, \(\alpha\) and \(\beta_{1:K}\) are parameters supplied to the simulator. A Bayesian version could add another level that generates the \(\beta_k\) from a prior distribution.
| Piece | Role | Significance |
|---|---|---|
| \(N_d\) | Supplied document length | The model shown is conditional on document lengths. |
| \(\alpha\) | Shared parameter | Controls how concentrated or diffuse document topic mixtures tend to be. |
| \(\Theta_d\) | Document-level latent variable | Gives the topic proportions for one document and is generated anew for each document. |
| \(\beta_k\) | Shared topic parameter | Gives the distribution of words associated with topic \(k\). |
| \(Z_{dn}\) | Token-level latent variable | Selects which topic generates a particular token. |
| \(W_{dn}\) | Observed categorical variable | The identity of the word at position \(n\). |
Conditional on \(\Theta_d\), the topic indicators within a document are independent. Conditional on their topic indicators, the words are also independent. The model therefore treats a document as a bag of words: permuting the word order does not change its probability.
The topic parameters \(\beta_k\) may reveal recognizable groups of co-occurring words, and \(\Theta_d\) may summarize the mixture of those groups within a document. Those quantities are closer to potential objects of substantive interpretation than the individual weights of a neural language model.
However, the model itself guarantees only that they play particular probabilistic roles. It does not guarantee that a topic corresponds to one real social, scientific, or semantic concept. Topics may split a single concept, combine several concepts, reflect writing style rather than subject matter, or change across model specifications. Topic labels are interchangeable without further constraints, and the fitted decomposition can be unstable or weakly identified.
The topic model illustrates a useful middle ground: a generator can be knowingly unrealistic at the level of complete text while still being valuable because its latent structure targets a useful simplification.
In this example, we introduce “latent states,” which are quantities generated within each simulated dataset but not directly observed. A parameter is ordinarily held fixed while a dataset is generated. Both can be unknown when we analyze observed data; the distinction is their role in the generative procedure, not simply whether we know them.
As an example, suppose an animal moves through a landscape. At time \(t\), let \(Z_t\) denote an unobserved behavioral state such as resting, foraging, or traveling. Let \(Y_t\) contain measured features of movement, such as step length and turning angle.
A hidden Markov model can be written as
\[ Z_1\mid\pi \sim \operatorname{Categorical}(\pi), \]
\[ Z_t\mid Z_{t-1}=j;A \sim \operatorname{Categorical}(A_{j,\cdot}), \qquad t=2,\ldots,T, \]
\[ Y_t\mid Z_t=k;\phi_{1:K} \sim F_k(\phi_k), \qquad t=1,\ldots,T. \]
The notation \(F_k(\phi_k)\) stands for the observation distribution associated with state \(k\). For animal movement, one might use a Gamma distribution for step length and a circular distribution for turning angle, with state-specific parameters collected in \(\phi_k\).
| Piece | Role | Significance |
|---|---|---|
| \(\pi\) | Initial-state parameter | Gives the probabilities of the possible state at the beginning of a trajectory. |
| \(A\) | Transition parameter | Row \(j\) describes how the next state is generated when the current state is \(j\). Large diagonal entries produce persistent states. |
| \(\phi_k\) | Observation parameter for state \(k\) | Describes the movement measurements expected in that state. |
| \(Z_t\) | Latent state | Generated anew at each time but linked over time through the Markov model. |
| \(Y_t\) | Observation | Noisy information about the current state. |
The model makes two central conditional-independence claims:
These claims are what allow the joint law to factor as
\[ P(z_{1:T},y_{1:T}) =P(z_1) \prod_{t=2}^T P(z_t\mid z_{t-1}) \prod_{t=1}^T P(y_t\mid z_t). \]
While simulating one trajectory, \(A\) and \(\phi_{1:K}\) are held fixed, whereas \(Z_{1:T}\) is generated. After observing \(Y_{1:T}\), we may be uncertain about all of them. We may want to infer the shared transition behavior, the observation distributions, and the particular state sequence that occurred.
The state labels themselves are created by the model. Calling a state “foraging” is a scientific interpretation, not a consequence of the equations. A state may instead capture unmodeled habitat differences, measurement artifacts, or several behaviors combined. As in a topic model, the labels can be permuted without changing the probability law.
A neural language model generates text one token at a time. Given a prompt \(w_{1:m}\), it specifies
\[ W_t\mid W_{<t}=w_{<t} \sim \operatorname{Categorical}\!\left(\boldsymbol{\pi}_\theta(w_{<t})\right), \qquad t=m+1,\ldots,T, \]
and therefore
\[ P_\theta(w_{m+1:T}\mid w_{1:m}) = \prod_{t=m+1}^T P_\theta(w_t\mid w_{<t}). \]
The next-token probabilities are constructed by a neural network. Schematically,
\[ z_s=E_{w_s}+p_s, \]
\[ h_t=\operatorname{Transformer}_\theta(z_1,\ldots,z_{t-1}), \]
and
\[ \boldsymbol{\pi}_\theta(w_{<t}) = \operatorname{softmax}(W_{\mathrm{out}}h_t+b_{\mathrm{out}}). \]
Here \(E_{w_s}\) is the learned embedding of token \(w_s\), while \(p_s\) records its position. The transformer combines the preceding token representations into a context-dependent representation \(h_t\), using a causal restriction that prevents it from accessing future tokens. The output map and softmax turn this representation into a categorical probability law over the next token. Once the context and parameters \(\theta\) are fixed, these probabilities are calculated deterministically; randomness enters when the next token is drawn and appended to the context. The parameters \(\theta\) comprise a large collection of learned embeddings and network weights. They are useful because together they produce effective conditional distributions, not because individual parameters ordinarily have scientific interpretations.
| Example | Main generated object | Important latent quantities | Role of parameters | Central lesson |
|---|---|---|---|---|
| Linear regression | Response vector conditional on covariates | Usually none in the elementary form | Coefficients and noise scale summarize a conditional law | A mean function is only one part of a generator. |
| Golf: logistic | Success counts by distance | None after success probabilities are computed | Empirical curve parameters | A complete empirical generator need not explain the physical mechanism. |
| Golf: angular error | Individual aiming errors and success counts | Angular errors for individual putts | Directional-error scale has mechanistic meaning | Mechanism supplies structure but creates additional assumptions to check. |
| Markov language model | Token sequence | Generated tokens become later inputs | Transition probabilities are shared local rules | Ordered conditionals define a joint law. |
| Exam-response regressions | One response conditional on all others | None specified | Predictive coefficients | Circular componentwise predictions need not define a joint generator. |
| Topic model | Documents as bags of words | Document mixtures and token topics | Topics may be substantively interpreted | An unrealistic generator can define useful latent summaries. |
| Hidden Markov model | Time series of observations | Persistent state sequence | Transition and emission parameters are shared across a trajectory | Parameter versus latent state is a distinction of generative role. |
| Neural language model | Token continuation conditional on a prompt | Internal activations are computed, not ordinarily treated as sampled latent states | Vast collections of instrumental weights | Generative models need not have individually interpretable parameters. |
Across all of these cases, calling a model generative tells us that it defines a coherent probability law and a way to simulate the relevant data. It does not by itself tell us that the model is realistic, causal, scientifically interpretable, identifiable, computationally convenient, or appropriate for the question at hand. Those are separate claims that require separate arguments and checks.