From Models to Simulators

What is the correspondence between a probabilistic model and a simulator?

Every model in this course is a probabilistic recipe: a description of how data could arise.

Today: for each of several examples, we see

1. What do the random variables in the model actually represent?
2. How does that recipe become code, line by line?

Ask the same questions of every model

1. What constitutes one complete simulated dataset?
2. What inputs and parameters are supplied?
3. Which quantities are drawn randomly?
4. Which quantities are calculated from previously generated values?
5. What is shared across observations or times?
6. What is returned, and what remains unobserved?

We’ll go through five examples, roughly in order of increasing complexity:

• independent draws
• draws that depend on each other
• a model with a variable we generate but never observe (a latent variable model)
• a hierarchical mixture (also with latent variables)
• a system that evolves over time.

Example: Bernoulli trials and their sum

A rather silly example, just to illustrate the principles.

Experiment to match: someone flips a coin \(n\) times and tells us only how many came up heads.

\[ X_i\overset{\text{iid}}\sim\operatorname{Bernoulli}(p),\quad i=1,\ldots,n, \qquad Y=\sum_{i=1}^n X_i. \]

x <- rbinom(n, size = 1, prob = p)
y <- sum(x)

x = rng.binomial(1, p, size=n)
y = x.sum()

We simulate every trial, then return only the total, matching what the experiment actually reports. \(Y\) is random even though the sum is a deterministic function because that function takes random variables as arguments.

Model-to-code correspondence: Bernoulli example

By the way: Summing \(n\) iid \(\operatorname{Bernoulli}(p)\) trials is common enough that its exact distribution has been worked out: \(Y\) is exactly \(\operatorname{Binomial}(n,p)\), and both R (rbinom(1, n, p)) and Python (rng.binomial(n, p)) let you draw it directly. Much of what we build this semester won’t be common enough to have that kind of built-in answer, in which case we’ll have to build the models more explicitly.

Example: a Markov (bigram) language model

Markov language models are examples of a generative model for which a fairly parsimonious description can still generate long, varied output.

\(W_t\) is the actual word generated at position \(t\). \(A\) is a transition table, estimated from a passage of text; its row \(A_{w,\cdot}\) gives the conditional distribution of the next word given that the current word is \(w\).

\[ W_1 = w_{\text{start}}, \qquad W_t \mid W_{t-1} \sim \operatorname{Categorical}\big(A_{W_{t-1},\,\cdot}\big), \quad t=2,\ldots,T. \]

\(w_{\text{start}}\) is a supplied starting word, and so \(W_1\) is not random here.

We would normally estimate \(A\) from real text (shown in code associated with the lecture), then treat it as a fixed, supplied parameter. The estimation is not normally done with the sorts of methods described in this course.

generated[1] <- start_word
for (t in 2:n_tokens) {
  current <- generated[t - 1]
  probs <- transition_prob[current, ]
  generated[t] <- sample(vocab, 1, prob = probs)
}

generated = [start_word]
for _ in range(n_tokens - 1):
    current = generated[-1]
    words, probs = zip(
        *transition_prob[current].items())
    generated.append(rng.choice(words, p=probs))

Model-to-code correspondence: Markov language model

The current word → sample next word → append → repeat loop is how we simulate from the model.

Latent variable models

In the Markov model, everything we generated was returned. That’s not always true.

A generative model can include latent variables (quantities that are defined by the model/simulator but never returned/observed) either because they make the rule simpler or more realistic to state.

The hidden Markov model is our first example: an unobserved state sequence drives what we actually observe.

Example: a hidden Markov model

Concrete version: \(Z_t\) is which regime the economy is in this quarter (recession or expansion) and \(Y_t\) is a noisy indicator we actually observe, like GDP growth, whose typical level and volatility depend on the regime. We see \(Y_t\); the regime \(Z_t\) itself is never directly observed, only inferred. This is related to a real macroeconomic model.

\[ Z_1\sim\operatorname{Categorical}(\pi), \qquad Z_t\mid Z_{t-1}\sim\operatorname{Categorical}(A_{Z_{t-1},\cdot}), \]

\[ Y_t\mid Z_t=k\sim\operatorname{Normal}(\mu_k,\sigma_k^2). \]

z[1] <- sample.int(n_states, 1, prob = pi)
y[1] <- rnorm(1, mu[z[1]], sigma[z[1]])
for (t in 2:n_steps) {
  z[t] <- sample.int(n_states, 1,
                      prob = A[z[t - 1], ])
  y[t] <- rnorm(1, mu[z[t]], sigma[z[t]])
}

z[0] = rng.choice(n_states, p=pi)
y[0] = rng.normal(mu[z[0]], sigma[z[0]])
for t in range(1, n_steps):
    z[t] = rng.choice(n_states,
                       p=A[z[t - 1], :])
    y[t] = rng.normal(mu[z[t]], sigma[z[t]])

Model-to-code correspondence: HMM

Example: a topic model

Topic models are generative models that don’t give a great description of all of the structure in the data they are attempting to model (documents), but they nonetheless can give good summaries.

\(\theta_d\) represents how much document \(d\) is about each topic (e.g., 70% sports, 30% cooking). \(z_{d,n}\) represents which single topic word \(n\) actually came from. \(w_{d,n}\) is the word itself.

\[ \underbrace{\theta_d\sim\operatorname{Dirichlet}(\alpha)}_{\text{topic mixture}} \quad\longrightarrow\quad \underbrace{z_{dn}\sim\operatorname{Categorical}(\theta_d)}_{\text{topic for word }n} \quad\longrightarrow\quad \underbrace{w_{dn}\sim\operatorname{Categorical}(\beta_{z_{dn}})}_{\text{observed word}}. \]

theta_d <- rdirichlet1(alpha)
z <- sample.int(n_topics, n_words_per_doc,
                replace = TRUE, prob = theta_d)
for (n in seq_len(n_words_per_doc))
  w[n] <- sample(colnames(beta), 1,
                 prob = beta[z[n], ])

theta_d = rng.dirichlet(alpha)
z = rng.choice(n_topics, size=n_words_per_doc,
               p=theta_d)
w = [rng.choice(vocab, p=beta[z[n], :])
     for n in range(n_words_per_doc)]

Model-to-code correspondence: topic model

Base R has no rdirichlet(); we draw one by drawing independent Gammas and normalizing — a standard construction.

Example: a discrete-time SIR epidemic model

This model correspondsn to an actual assumed scientific mechanism of how infection spreads through a population.

\(S_t, I_t, R_t\) are counts of people in each health state at time \(t\). \(\Delta I_t\) is how many people newly became infected this step; \(\Delta R_t\) is how many newly recovered.

\[ \underbrace{S}_{\text{Susceptible}} \xrightarrow{\ \beta SI/N\ } \underbrace{I}_{\text{Infectious}} \xrightarrow{\ \gamma I\ } \underbrace{R}_{\text{Recovered}} \]

This continuous-time picture is the usual way SIR is drawn, but it hides exactly how event times are generated. The discrete-time version here keeps the same story and is more transparent to simulate from.

A discrete-time SIR model: two binomial draws per step

At each time step \(t = 0, 1, \ldots\), with population size \(N = S_t+I_t+R_t\) fixed:

\[ \Delta I_t \mid S_t, I_t \sim \operatorname{Binomial}\!\Big(S_t, 1-(1-\beta/N)^{I_t}\Big), \qquad \Delta R_t \mid I_t \sim \operatorname{Binomial}(I_t, \gamma), \]

Intuition for the success probability: each of the \(I_t\) infectious people independently infects a given susceptible person with probability \(\beta/N\) this step, so \((1-\beta/N)^{I_t}\) is the chance that person escapes all of them, and \(1-(1-\beta/N)^{I_t}\) is the chance they get infected by at least one.

\[ S_{t+1} = S_t - \Delta I_t, \qquad I_{t+1} = I_t + \Delta I_t - \Delta R_t, \qquad R_{t+1} = R_t + \Delta R_t. \]

for (t in seq_len(n_steps)) {
  p_infect <- 1 - (1 - beta / N)^I[t]
  delta_I <- rbinom(1, S[t], p_infect)
  delta_R <- rbinom(1, I[t], gamma)
  S[t + 1] <- S[t] - delta_I
  I[t + 1] <- I[t] + delta_I - delta_R
  R[t + 1] <- R[t] + delta_R
}

for t in range(n_steps):
    p_infect = 1 - (1 - beta / N) ** I[t]
    delta_I = rng.binomial(S[t], p_infect)
    delta_R = rng.binomial(I[t], gamma)
    S[t + 1] = S[t] - delta_I
    I[t + 1] = I[t] + delta_I - delta_R
    R[t + 1] = R[t] + delta_R

Model-to-code correspondence: discrete-time SIR

Today’s lab: a single-server queue

• One server, first-come-first-served.
• System starts empty.
• No abandonment, no capacity limit.
• Customer \(i\)’s interarrival duration \(A_i\sim\operatorname{Exponential}(\lambda)\) and service duration \(S_i\sim\operatorname{Exponential}(\mu)\), all independent.
• A fixed number of customers arrive; every one is followed through departure.

Probabilistic description of the model

\[ \begin{aligned} T_i^{\text{arr}} &= \textstyle\sum_{j\le i} A_j, \\ T_i^{\text{start}} &= \max\!\left(T_i^{\text{arr}},\, T_{i-1}^{\text{dep}}\right), \\ T_i^{\text{dep}} &= T_i^{\text{start}} + S_i, \\ W_i &= T_i^{\text{start}} - T_i^{\text{arr}}. \end{aligned} \]

The \(\max(\cdot)\) is the key line: customer \(i\)’s service can’t start before either their own arrival or the previous customer’s departure — whichever comes later.

What the lab asks you to do

1. Read the simulator — match each line of code to a piece of the recurrence.
2. Explore behavior — vary the arrival rate, watch what happens to waiting times.
3. Repeat a setting — look at the distribution of mean waiting time across queues, not just one run.

Optional: hit a waiting-time target, push \(\lambda\ge\mu\), try constant service times.