What should we match?

Today

Last time, we saw what happened when we used poor summaries or not enough summaries.

This time, we’ll look at adding more summaries and we’ll see other general strategies for summaries to pick.

Specifically:

1. Look at what happens when we add more and more summaries.
2. Compare the entire empirical distributions instead of hand-picked summaries (maximum mean discrepancy).
3. What to do when we can’t do that — fit a small auxiliary model and use its parameters as summaries (indirect inference).

Section 0: Why don’t we just use the best summary?

• Stat theory has the idea of “sufficiency” of a a summary for a particular model. Idea is there is no loss of information going from full data to the summary.
• Therefore, we should use sufficient statistics!
• But if you are in a situation where you know what the sufficient statistics are, you probably don’t want to use SBI.

Section 1: Why not keep adding summaries?

Summaries can behave in various ways:

• Informative: changes appreciably relative to its baseline fluctuations as that parameter direction moves.
• Redundant: moves too, but tracks information another summary already gave us.
• Weak: barely moves at all, relative to its baseline fluctuations.

Roughly,

\[ \hat\theta - \theta_0 \;\approx\; \frac{\text{summary noise}}{\text{sensitivity of the expected summary at } \theta_0}. \]

A summary helps to the extent its numerator is small and its denominator is large.

Why not just use a lot of summaries?

If any single summary might be redundant or weak, the obvious thing is to just add everything you can think of.

For instance:

• means, variances, skewnesses, kurtoses, higher-order moments
• quantiles;
• histogram-bin frequencies;
• autocorrelations;
• regression coefficients fit to the data;
• selected points along a trajectory;
• slopes or curvature of a fitted curve.

Nothing stops you from computing all of these and stacking them into one long \(s(Y)\).

More summaries, in principle

Simulated method of moments comes from a non-simulation-based method called GMM (generalized method of moments).

If you’re not simulating, and you use the optimal weighting matrix, you can just add more and more summaries and get a better estimator.

The problem is that if you’re simulating, adding non-informative/weakly-informative summaries also adds noise, and so summaries have to be good compared to the amount of noise they contribute.

A worked example

Recall estimating \((\beta,\gamma)\) in the SIR model.

When we used final size and early growth rate, we got an ok estimate of both \(\beta\) and \(\gamma\).

What if we add a third summary: \(I_2\), the raw case count on day 2?

Is \(I_2\) informative, redundant, or weak?

• correlation with growth rate: 0.41; with final size: 0.04, so it’s not entirely redundant with either
• sensitivity-to-noise ratio of \(I_2\) for \(\beta\) is \(\approx 8.8\), versus growth rate’s \(\approx 23.5\) (matching the formula from Section 1: a finite-difference approximation to \(\frac{dE_\theta[s]}{d\beta}\) at the true \(\beta\) — bump \(\beta\) slightly, divide the change in the summary’s mean by the size of that bump — divided by the summary’s own noise, its SD)

Including \(I_2\) degrades performance

Using \(\hat\Omega^{-1}\) estimated from \(B=10\) simulated replicates at the true \((\beta,\gamma)\):

Section 2: A different idea — compare distributions directly

Stepping back

SMM reduces \(Y\) to a short, hand-picked vector \(s(Y)\), and compares vectors. Section 1 just showed how easy it is to get that choice wrong — even a genuinely valid, non-redundant summary can hurt once you have to estimate how to weight it.

What if we didn’t have to hand-pick a summary vector at all? What if we asked a more direct question:

How far apart are the whole simulated and observed distributions?

We’ll build this up from the simplest possible case: outcomes that fall into a small number of categories.

The household outbreak model

Suppose we have a bunch of independent SIR-generated ’epidemics,” each one within a household of size \(K\).

We assume that each household under consideration is seeded with one index case, we think transmission within the household follows the SIR model, and we get to observe the total number of people infected in each household.

\[ S_i(0) = K-1, \qquad I_i(0) = 1, \qquad R_i(0) = 0, \]

\[ \Delta I_i(t) \mid S_i(t), I_i(t) \sim \operatorname{Binomial}\!\Big(S_i(t),\, 1-(1-\beta/K)^{I_i(t)}\Big), \qquad \Delta R_i(t) \mid I_i(t) \sim \operatorname{Binomial}(I_i(t), \gamma), \] \[ X_i \;=\; \big[K - S_i(\infty)\big] \;\in\; \{1,\ldots,K\}, \]

We assume households are independent (unreasonable), so \(X_1,\ldots,X_n\) are iid draws from one distribution determined by \(\theta=(\beta,\gamma)\).

Comparing observed and simulated data using histograms

From \(n\) observed households we get empirical proportions \(\hat p_j\) (fraction with \(j\) people infected); from \(n\) simulated households at some candidate \(\theta\), empirical proportions \(\hat q_j\):

One answer

Treat \(\hat{\mathbf p}\) and \(\hat{\mathbf q}\) as two points in \(\mathbb R^8\) and take squared Euclidean distance:

\[ d^2(\hat{\mathbf p}, \hat{\mathbf q}) = \sum_{j=1}^{8} (\hat p_j - \hat q_j)^2. \]

For the histograms above (observed: \(R_0=2.5\); candidate \(\theta\): \(R_0=1.5\), from \(n=300\) households each),

\[ d^2 = \sum_j (\hat p_j - \hat q_j)^2 \approx 0.108. \]

(For comparison: two independent samples of 300 households at the same \(\theta\) as the observed data give \(d^2\approx 0.009\) — noise alone, about a twelfth as large. The candidate above is genuinely far off, not just noisy.)

MMD: A metric for distributions

For any kernel \(k\), and observed/simulated samples \(x_1,\ldots,x_n\) and \(y_1,\ldots,y_m\),

\[ \widehat{\operatorname{MMD}}^2 = \underbrace{\frac{1}{n^2}\sum_{i,i'} k(x_i,x_{i'})}_{\text{average within-observed similarity}} \;+\; \underbrace{\frac{1}{m^2}\sum_{j,j'} k(y_j,y_{j'})}_{\text{average within-simulated similarity}} \;-\; \underbrace{\frac{2}{nm}\sum_{i,j} k(x_i,y_j)}_{\text{average between-sample similarity}}. \]

We think of a kernel as being a function that takes two values and measures how similar they are.

We will start with the equality kernel, which says that values have similarity 1 if they are equal and similarity 0 otherwise.

We will then generalize to kernels that don’t require exact matches.

Checking the definition against the histogram case

Define the equality kernel on our outcome space:

\[ k(x,x') = \mathbf 1\{x = x'\}. \]

Two observations are similarity 1 if they land in the same category, similarity 0 otherwise. Plug this into the general definition. Among the \(n^2\) pairs \((x_i,x_{i'})\), the fraction landing in the same category \(j\) is \(\hat p_j^2\), so

\[ \frac{1}{n^2}\sum_{i,i'} k(x_i,x_{i'}) = \sum_j \hat p_j^2, \qquad \frac{1}{nm}\sum_{i,j} k(x_i,y_j) = \sum_j \hat p_j \hat q_j, \]

so

\[ \widehat{\operatorname{MMD}}^2 = \sum_j \hat p_j^2 + \sum_j \hat q_j^2 - 2\sum_j \hat p_j \hat q_j = \sum_j (\hat p_j-\hat q_j)^2. \]

SMM on the category frequencies

Our dataset is \(X_1,\ldots,X_n\) (here, \(n\) households). Define one summary per possible value of \(X_i\):

\[ s_j(X_1,\ldots,X_n) = \frac{1}{n}\sum_{i=1}^n \mathbf 1\{X_i=j\} = \hat p_j, \qquad j=1,\ldots,8. \]

Now run ordinary SMM with this eight-dimensional summary vector and the simplest possible weighting, \(W=I\):

\[ Q(\theta) = g(\theta)^\top W g(\theta) = \sum_j(\hat p_j-\hat q_j)^2, \qquad g(\theta) = \hat{\mathbf p}-\hat{\mathbf q}(\theta), \quad W=I. \]

Then \(Q(\theta)\) is exactly MMD between observed data and data simulated with parameter value \(\theta\).

Pros and cons

Cons:

• If we interpret this as SMM with a particular set of summaries and \(W = I\), our previous observations about using \(W = I\) go through. Higher-probability categories are going to have more weight, summaries are not going to be combined optimally.
• We are kind of using a lot of summaries, and we have seen that adding a lot of extra summaries isn’t always useful.

Pros:

• We don’t have to hand-build summaries.
• We will be able to generalize it to continuous data and less stringent sorts of kernels.

Section 3: Smooth kernels

The equality kernel is very stringent/doesn’t work for continuous data

• \(k(x,x')=\mathbf 1\{x=x'\}\) says: same category, similarity 1; anything else, similarity 0
• for our household example, a household with 2 people infected is exactly as “different” from one with 3 as it is from one with 8.
• maybe we think a household with 2 people infected is more similar to a household with 3 than it is to one with 8.
• the equality kernel won’t work at all for continuous data, where every pair of points will have similarity 0.

Therefore, we generalize to kernels that are less stringent, for instance the Gaussian kernel:

\[ k_h(x,x') = \exp\left\{-\frac{(x-x')^2}{2h^2}\right\}. \]

• \(k_h(x,x')\) is close to 1 when \(x,x'\) are close;
• it decays toward 0 as \(|x-x'|\) grows;
• \(h\) (the bandwidth) sets the scale of “close.”

What the two kernels look like

The same formula, a richer kernel

For a finite set of outcomes, any kernel can be written as a similarity matrix \(K\), and

\[ \operatorname{MMD}^2(P,Q) = (\mathbf p - \mathbf q)^\top K (\mathbf p - \mathbf q). \]

• With \(K=I\) (the equality kernel), this collapses back to \(\sum_j(p_j-q_j)^2\) (this is just a different way of writing the same thing)
• With a smooth \(K\), a mismatch that moves mass to a neighboring category is weighted differently than one that moves it far away.
• This vector-and-matrix picture only makes sense for a finite set of categories. But the underlying MMD formula — sums of \(k(x_i,x_{i'})\) over pairs of observations, a few slides back — never actually needed a finite outcome space; it’s computable directly from any samples, continuous ones included.

Section 5: A contrived example showing MMD with continuous data

Let

\[ S_i \sim \operatorname{Categorical}(\{-1,1\}), \quad P(S_i=-1)=P(S_i=1)=\tfrac12, \qquad \varepsilon_i \sim N(0, 1-a^2), \]

\[ Y_i = aS_i + \varepsilon_i, \qquad 0 \le a < 1. \]

This is constructed so that the mean and variance don’t change with \(a\).

What sorts of data does the model produce for different values of \(a\)?

Two different methods for inferring \(a\)

Suppose the true value is some \(a_0>0\), unknown to us, and we build two SMM-style objectives from one fixed observed sample:

• \(Q_{\text{mean,var}}(a)\): matches simulated mean and variance to observed mean and variance.
• \(Q_{\text{MMD}}(a)\): the Gaussian-kernel MMD between observed and simulated samples.

The actual objectives

So have we completely solved the problem?

The MMD curve you just saw used one particular bandwidth. What happens at other values of \(h\)?

Reading the three panels

• Too small (\(h=0.03\)): only nearly-identical pairs register as similar at all. You can see the true signal, but it’s noisier than it is for the larger bandwidth.
• Good (\(h=0.3\)): less noisy than the smaller bandwidth, with a more convincing minimum.
• Too large (\(h=3\)): the kernel defines every pair of points to be almost equally similar to each other, which washes out the signal

Section 7: Indirect inference

What if you only have one object?

For MMD to make sense, your dataset needs to consist of repeated samples from the same distribution.

What if you only have one sample, like a macroeconomic time series or a single epidemic trajectory?

Indirect inference, formally

A structural model \(Y\sim P_\theta\) (SIR, in our case) has the parameter \(\theta\) we actually want. An auxiliary model is any fitting procedure \(\hat\alpha(\cdot)\) that turns a dataset into a vector of fitted numbers (usually the parameters of the model).

1. Fit the auxiliary model to the observed data: \(\hat\alpha_{\text{obs}} = \hat\alpha(Y_{\text{obs}})\).
2. At a candidate \(\theta\), simulate \(Y_\theta^{(1)},\ldots,Y_\theta^{(B)}\sim P_\theta\) and fit the same auxiliary model to each: \(\hat\alpha_\theta^{(b)} = \hat\alpha\big(Y_\theta^{(b)}\big)\).
3. Average: \(\overline{\hat\alpha_\theta} = \frac{1}{B}\sum_{b=1}^B \hat\alpha_\theta^{(b)}\).
4. Choose

\[ \hat\theta = \arg\min_\theta \; g(\theta)^\top W\, g(\theta), \qquad g(\theta) = \hat\alpha_{\text{obs}} - \overline{\hat\alpha_\theta}. \]

That’s exactly SMM (Lecture 4) with \(s(Y) = \hat\alpha(Y)\): an auxiliary model’s fitted coefficients standing in for a hand-picked summary vector.

Structural vs. auxiliary models

• The structural model, concretely: SIR, with \(\theta=(\beta,\gamma)\).
• The auxiliary model: any convenient, easy-to-fit descriptive model. Its only job is to produce a handful of numbers that summarize the shape of a dataset well.

For an epidemic curve, a natural auxiliary model is a quadratic fit to the log case counts:

\[ \log(I_t + 1) \;\approx\; \alpha_0 + \alpha_1 t + \alpha_2 t^2. \]

\(\alpha_0\) tracks the starting level, \(\alpha_1\) the early growth rate, \(\alpha_2\) the curvature — whether growth is accelerating or already bending over.

The auxiliary model is not supposed to describe how we actually think the world behaves, it’s just supposed to be a useful summary of the dataset.

What the auxiliary fit looks like

One simulated epidemic (points), and the quadratic fit to \(\log(I_t+1)\) (line).

\((\hat\alpha_0,\hat\alpha_1,\hat\alpha_2)\) would then be used as summaries for this particular simulation.

The algorithm, for this example

1. \(\hat\alpha_{\text{obs}}\): fit the quadratic to the one observed epidemic curve.
2. At a candidate \((\beta,\gamma)\), simulate \(B\) epidemics and fit the quadratic to each, getting \(\hat\alpha_{(\beta,\gamma)}^{(1)},\ldots,\hat\alpha_{(\beta,\gamma)}^{(B)}\).
3. Average the \(B\) fitted coefficient vectors.
4. Move \((\beta,\gamma)\) until that average resembles \(\hat\alpha_{\text{obs}}\).

The auxiliary model as a measuring instrument

A good auxiliary model should:

• respond to the structural parameters you actually care about;
• describe features that vary systematically across simulated datasets, not just noise;
• be stable enough to fit the same way, every time, without babysitting;
• separate parameter directions that need to be told apart.

Indirect inference doesn’t make summary choice disappear

The auxiliary model is still a choice, and it can still go wrong in exactly the ways Lecture 5 warned about:

• a badly chosen auxiliary model can be redundant with itself (its own fitted coefficients highly correlated, adding little beyond the first);
• it can be weak for a structural parameter that never shows up in its fitted output;
• it systematizes how you build \(s(Y)\), it does not guarantee that \(s(Y)\) is informative.

Summing up

Now we have a bunch of different ways to answer the question we asked in the first class — does the simulation resemble the data? — by picking a representation \(T(Y)\) and a discrepancy \(d\):

\[ \hat\theta = \arg\min_\theta \; d\big(T(Y_{\text{obs}}),\, T(Y_\theta)\big). \]

Some are more automatic than others, but they all require choices.