Simulated Method of Moments

Today

How can we use the knowledge we gained by simulating under different parameter values to estimate parameter values from real data?

Section 1: Wrapping up the previous lab

What did we learn from the queue?

• Simulation can reveal behavior that isn’t obvious from the code alone.
• A stochastic model produces a range of possible outcomes, not one.
• System behavior can change nonlinearly as parameters change — waiting times blow up sharply as utilization approaches 1.

Forward simulation is useful even without inference

• Queueing: waiting-time behavior near saturation.
• Epidemics: whether outbreaks take off, and how large they get.
• Ecological models: extinction probability.
• Demographic models: how heritability of a trait like fertility shapes population growth or decline across generations.
• Generative models for proteins or language: what the model “imagines” before it’s ever fit to real data.
• Mechanistic scientific simulators: qualitative consequences of assumptions.

Sometimes the question is simply: “If the world worked this way, what would happen?”

Forward simulation as a component of model fitting

In this course, we want to fit models. To do so, we need to know about the relationship between the parameters and observable aspects of the data.

For instance, we need to know things like:

• Which summaries move when the parameters change?
• Are any summaries more informative than others?
• Are there any parameters of the model that don’t lead to any observable differences in the simulations?

For simple models, it is possible to do out the math to find the answers to these questions.

In more complicated models, the only way to find out is to look at simulations across a range of parameter values.

Section 2: Investigating the behavior of a generative model

The SIR model

Overall idea: have a population of individuals, they start out susceptible, potentially become infected, and then recover.

\[ S \longrightarrow I \longrightarrow R \]

• \(S_t\): number susceptible at time \(t\).
• \(I_t\): number infected at time \(t\).
• \(R_t\): number recovered at time \(t\).

\[ \beta = \text{transmission parameter}, \qquad \gamma = \text{recovery parameter}, \qquad R_0 = \frac{\beta}{\gamma} \]

\(R_0\): roughly, how many new infections one infectious person generates while population immunity is still negligible (not immediately obvious from the framing of the model).

How the epidemic evolves

Discrete time steps, two binomial draws per step:

\[ \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). \]

\(\Delta I_t\) people move from S to I; \(\Delta R_t\) people move from I to R:

\[ 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. \]

One simulated epidemic

What actually happens inside simulate_sir, one time step at a time:

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

What does that produce?

Different parameters, different behavior

Two different \((\beta, \gamma)\) pairs.

The shape of the trajectory clearly changes as the parameters change.

However: variability from run to run with the same parameters makes it difficult to investigate systematically over different parameter values.

Therefore: generally take some summary statistics that initial investigation suggests might be informative and look at how they vary with the parameters in a more systematic way.

Some candidate summaries

Based on initial plots, we might use:

• final epidemic size
• early growth rate
• time to peak
• peak number infected

Let’s just investigate one parameter

To make things simpler, we will fix \(\gamma\) (a story would be that it is known from clinical data on infectious duration).

We will try to investigate the behavior of the model as a function of \(R_0 = \beta / \gamma\).

We will:

• Look at the distribution of a candidate summary as a function of \(R_0\).
• Ask whether the summary helps us distinguish different values of \(R_0\).

Does early growth rate distinguish different \(R_0\) values?

Yes — centers move strongly, and distributions overlap relatively little. Early growth rate is a useful summary.

Does final epidemic size distinguish different \(R_0\) values?

Section 3: Turn the simulator around

We have good evidence that our summaries distinguish different values of \(R_0\). Let’s use that.

Forward and inference problems

Forward: choose a parameter, simulate, summarize. Inference reverses the question — we observe the summary and ask what parameter could have produced it.

\[ \theta \xrightarrow{\ \text{simulate}\ } Y \xrightarrow{\ s(\cdot)\ } s(Y) \qquad\qquad s(Y_{\text{obs}}) \xrightarrow{\ \text{estimate}\ } \theta \]

For SIR:

\[ R_0 \xrightarrow{\ \text{simulate}\ } \text{epidemic} \xrightarrow{\ s(\cdot)\ } \text{early-growth statistic} \qquad\qquad \text{observed early-growth statistic} \xrightarrow{\ \text{estimate}\ } R_0 \]

The most naive algorithm

Choose a set \(\Theta\) of possible parameter values. For each \(\theta \in \Theta\):

1. Simulate many datasets.
2. Compute the same summary for each.
3. Average the simulated summaries.
4. Compare with the observed summary.

Choose the value that is closest to the observed summary.

The naive algorithm leaves a lot unsaid

• How many simulated datasets is enough, and does it matter?
• What exactly do we mean by “closest” — especially once there is more than one summary?
• How do we turn “try a bunch of \(\theta\)’s” into something we can actually code up and defend?

We need to pin this down precisely enough to implement.

Section 4: Simulated Method of Moments

The standard way to formalize this idea is Simulated Method of Moments (SMM):

\[ \hat\theta = \arg\min_\theta \; g(\theta)^\top W \, g(\theta), \qquad g(\theta) = s(Y_{\text{obs}}) - \overline{s(Y_\theta)}, \qquad \overline{s(Y_\theta)} = \frac1B\sum_{b=1}^B s(Y_\theta^{(b)}) \]

Note that for a single summary, \(W\) is a \(1 \times 1\) matrix (a single scalar) and the result doesn’t depend on the value it takes.

The components we need

• \(s(Y_{\text{obs}})\) — the observed summary.
• \(\overline{s(Y_\theta)}\) — the simulated average summary under a candidate \(\theta\).
• \(g(\theta)\) — how far apart they are.
• \(W\) — a weighting matrix: it specifies how to combine multiple summaries. We’ll use \(W = \hat\Omega^{-1}\), where \(\hat\Omega\) is our estimate (from the simulations) of \(\operatorname{Cov}(g(\theta_0))\)

We’ll build these up one at a time, starting with the simplest case: one parameter, one summary.

The one-summary case

With a single summary, \(g(\theta)\) and \(W\) are both scalars, and \(W\) drops out entirely:

\[ Q(\theta) = \big[s(Y_{\text{obs}}) - \overline{s(Y_\theta)}\big]^2, \qquad \hat\theta = \arg\min_\theta Q(\theta) \]

The algorithm is therefore: Pick the parameter whose simulations reproduce the observed summary most closely.

Minimizing \(Q(R_0)\)

\(Q(R_0)\), across a grid of candidate \(R_0\)’s:

What if we have more than one summary?

In the SIR example, we had two summaries: \[ s(Y) = \begin{pmatrix} s_1(Y) \\ s_2(Y) \end{pmatrix} = \begin{pmatrix} \text{growth rate} \\ \text{final size} \end{pmatrix}, \qquad g(\theta) = s(Y_{\text{obs}}) - \overline{s(Y_\theta)} \]

\(g(\theta)\) is now a vector.

SMM says that we should use \(W = \hat\Omega^{-1}\), where \(\hat \Omega\) is an estimate of \(\operatorname{Cov}(g(\theta_0))\).

Many problems here!

Suppose we took \(W = I\)

If our “estimate” of the covariance of the summary statistics is that they are independent and each have variance 1, that corresponds to \(W = I\).

Plugging that into the SMM objective function gives

\[ Q(\theta) = g_1(\theta)^2 + g_2(\theta)^2 \]

Why might this cause a problem?

Discrepancies should be judged relative to ordinary variability

The problem traces back to scale. Final size ranges over hundreds; growth rate ranges over tenths. A discrepancy of 10 is:

• huge, if the statistic normally varies by 1;
• trivial, if the statistic normally varies by 100.

Treating both raw discrepancies as equally important — as \(W=I\) does — silently lets the noisier, larger-scale summary dominate.

For independent summaries, the fix is to divide each by its own variance:

\[ Q(\theta) = \frac{g_1(\theta)^2}{\sigma_1^2} + \frac{g_2(\theta)^2}{\sigma_2^2} \]

Estimating the covariance \(\hat\Omega\) is exactly how we do this in general, and how we account for the summaries not being independent.

An example of equal weighting giving poor results

Repeated-estimation results, SIR example (true \(R_0=6\), where growth rate is the good summary and final size is much weaker):

\(W=I\) looked natural, but it is not generally a good idea: combining the two summaries naively bought us nothing — it just looks like the weaker summary.

Does weighting actually help?

Here we cheat and use a covariance estimate based on the true R_0.

There are more sophisticated ways of estimating that don’t rely on knowing the truth.

A weak/noisy summary can hurt if combined carelessly; weighting keeps it from having inappropriate influence — and here, legitimately combining two pieces of information beats using either alone.

There is theory behind this choice of \(W\)

Setting \(W = \operatorname{Cov}(g(\theta_0))^{-1}\) isn’t just a good heuristic we happened to verify empirically today.

Under standard conditions, it is the choice of \(W\) that minimizes the estimator’s variance among all valid weighting matrices — a classical result from generalized method of moments theory.

The derivation is beyond today’s scope. The takeaway: “weight by the inverse covariance” is the textbook default because it is provably the best you can do with these summaries, not an arbitrary convention.

Final takeaway

\[ \boxed{ \theta \to \text{simulate many datasets} \to \text{calculate summaries} \to \text{estimate expected summaries} } \]

\[ \boxed{ \text{choose } \theta \text{ so simulated summaries resemble observed summaries} } \]

SMM turns a forward simulator into a parameter-estimation method.