Finishing SMM + Parametric Bootstrap

Today

We have seen SMM + a couple of variants. First part of today will be wrapping up a couple of loose ends.

• How do we choose the weighting matrix \(W\)?
• How is this actually implemented?
• Asymptotic theory for SMM (not required knowledge for this course, but FYI)

Second half will be on the parametric bootstrap, which we can use for uncertainty quantification

How do we choose \(W\)?

Recall from Lecture 4: Ideally \(W=\Omega^{-1}\), where \(\Omega=\operatorname{Cov}_\theta(g(\theta_0))\), downweights noisy summaries and accounts for correlations between them.

The problem: \(\Omega\) depends on \(\theta_0\), which we don’t know.

Two-step estimation of \(W\)

1. Choose a simple initial weight matrix \(W_1\), e.g. the identity matrix, or a diagonal inverse-variance matrix.
2. Minimize the preliminary objective to get \(\hat\theta_1\).
3. Simulate many datasets at \(\hat\theta_1\), compute the summary vector for each, and estimate its covariance \(\hat\Omega(\hat\theta_1)\).
4. Set \(W_2=\hat\Omega(\hat\theta_1)^{-1}\).
5. Minimize again, holding \(W_2\) fixed, to get \(\hat\theta_2\).

In words: get into the right neighborhood with something ok (\(W\) is the identity or a reasonable diagonal matrix), learn the local scale and correlation of the summaries there, then refit using those values.

Iterative estimation of \(W\)

You can also keep refining the estimate of \(W\). After obtaining \(\hat\theta_k\):

1. Simulate many datasets at \(\hat\theta_k\) and estimate \(\hat\Omega(\hat\theta_k)\).
2. Set \[ W_{k+1}=\hat\Omega(\hat\theta_k)^{-1}. \]
3. Holding \(W_{k+1}\) fixed, compute \[ \hat\theta_{k+1} = \arg\min_\theta g(\theta)^\top W_{k+1}g(\theta). \]
4. Repeat until the estimates change very little or a maximum number of iterations is reached.

Potential issues

• If several summaries measure nearly the same thing, \(\hat\Omega\) may nearly singular, making its inverse unstable. You can combine or remove summaries, use only the diagonal, or regularize: \[ W=(\hat\Omega+\lambda I)^{-1}. \]
• Too few simulations: With \(q\) summaries, \(\hat\Omega\) contains \(q(q+1)/2\) distinct quantities to estimate. With \(B\leq q\) simulated datasets, the sample covariance matrix cannot be invertible, and \(B\) usually needs to be much larger than \(q\) for stable estimation.

Software exists

• Generic optimizers in R, Python, and Julia can minimize any user-written criterion, and will sometimes do ok (we’ll talk about fancier ideas for optimizing simulated objective functions later in the semester).
• There is a python package respy that appears to have a good interface for SMM. I do not believe R has anything similar.
• You still have to choose summaries, establish identifiability, etc.

Asymptotic normality of the SMM estimator (Not required)

Define the local sensitivity matrix

\[ G = \left.\frac{\partial E_\theta\{s(Y)\}}{\partial\theta^\top}\right|_{\theta=\theta_0}, \]

and let \(\Omega\) be the large-sample covariance of the discrepancy \(g(\theta_0)\). Then

\[ \sqrt n(\hat\theta-\theta_0) \ \dot\sim\ N(0, V_W), \qquad V_W = (G^\top WG)^{-1} G^\top W\Omega WG (G^\top WG)^{-1}. \]

Reading the formula

Under the optimal weighting from Section 1, \(W=\Omega^{-1}\), this collapses to

\[ V_{\text{opt}} = (G^\top \Omega^{-1} G)^{-1}, \qquad \operatorname{Var}(\hat\theta) \approx V_{\text{opt}}/n. \]

• \(G\): how strongly the expected summaries respond when parameters move — this is exactly the sensitivity we’ve been reading off objective-surface plots since Lecture 5.
• \(\Omega\): how much the summaries fluctuate, and how their fluctuations are related — the same covariance Section 1 just showed how to estimate.
• \(G\) rank-deficient: local nonidentification (Lecture 5’s ridges).
• \(G\) nearly rank-deficient: weak identification, large uncertainty.
• Larger \(n\): less sampling variability, as always.

A bit more on optimization

We want to find

\[ \hat\theta=\arg\min_\theta Q(\theta). \]

We have been mainly using grid search to find the minimum.

It is also possible to use some of the built-in optimization functions in R/python/etc on a simulated objective function, but there will generally be some issues.

Derivative-based methods

Derivative-based methods use the gradient

\[ \nabla Q(\theta) \]

to determine which direction is downhill.

If we do not provide derivatives, the optimizer may approximate them using nearby evaluations:

\[ \frac{Q(\theta+h)-Q(\theta)}{h}. \]

Examples include BFGS and L-BFGS-B.

These methods can be fast when the objective is smooth and its derivatives can be estimated accurately.

Derivative-free methods

Derivative-free methods compare objective values at several parameter settings and move toward the settings that look better.

For example, Nelder–Mead maintains a collection of points and moves, expands, or contracts them according to their objective values.

These methods are useful when derivatives are unavailable or the objective is not smooth.

However:

• they often require many objective evaluations;
• many are still local methods;
• they still rely on objective values being accurate enough to compare.

What should we do for now?

• Plot the objective, or low-dimensional slices of it, whenever possible.
• Use scientifically meaningful bounds or transformations.
• Try more than one starting value.
• Use a derivative-free method if the objective has no reliable gradient.
• Re-evaluate the best candidates with a larger, independent set of simulations.
• Compare simulations from the fitted model with the observed data.

Later in the course, we will return to methods designed more specifically for expensive, noisy, or rough simulation-based objectives.

Uncertainty quantification

We are in a statistics course, so we need to do this

We are going to start from zero

Some terminology

• An estimator is just a procedure \(h\) applied to data \(Y\): \(\hat \theta = h(Y)\). We usually interpret it as trying to estimate some property of a broader population, but the definition is just that it’s a function of the observed data.

• A sampling distribution is a thought experiment: If we were able to get new datasets from the same process, what would the resulting estimates look like?

• Formally, the sampling distribution of \(h\) is the distribution of \[ \hat\theta=h(Y) \qquad\text{when}\qquad Y\sim P_{\theta_0}. \]

How do we use sampling distributions?

The sampling distribution shows us how the estimate changes across repeated datasets.

We can use it to ask:

• How variable is the estimator? The standard deviation of the sampling distribution is called the standard error.

• Is the estimator systematically too high or too low? A systematic shift away from the true value is called bias.

• What does the uncertainty look like? The estimates might form a symmetric distribution, a skewed distribution, several clusters, or a long ridge.

A simple example

Suppose we want to estimate the probability \(p\) that a coin lands heads based on ten coin flips.

We use the proportion of heads out of our ten coin flips as the estimator

\[ \hat p=h(Y)=\frac{\text{number of heads}}{10}. \]

These are different datasets from the same process, so they give different estimates.

The sampling distribution

For this example, suppose we know the true value is \(p=0.6\). We can simulate the entire experiment many times:

1. Flip the coin 10 times.
2. Calculate \(\hat p\).
3. Save the estimate.
4. Repeat.

This histogram contains 3,000 estimates, not 3,000 individual coin flips.

Each estimate came from an entire dataset of 10 observations.

A different experiment, a different sampling distribution

For a dataset of \(n\) coin flips, our estimator is

\[ \hat p_n=h_n(Y_1,\ldots,Y_n) =\frac{1}{n}\sum_{i=1}^n Y_i. \]

Changing \(n\) gives us a different data-collection experiment and a different estimator:

• \(\hat p_{10}\) is computed from datasets containing 10 flips.
• \(\hat p_{50}\) is computed from datasets containing 50 flips.
• \(\hat p_{200}\) is computed from datasets containing 200 flips.

Each estimator therefore has its own sampling distribution.

In this example, experiments with more data produce sampling distributions that are more tightly concentrated around the true value.

What do we do if we don’t know the truth?

Here we could simulate the sampling distribution because we chose the true value \(p=0.6\) ourselves.

With real data, the true parameter is exactly what we do not know:

\[ Y\sim P_{\theta_0}, \qquad \theta_0\text{ unknown}. \]

The parametric bootstrap replaces the unknown \(\theta_0\) with our estimate \(\hat\theta\).

Parametric bootstrap

In the parametric bootstrap, we use \(P_{\hat\theta}\) as our stand-in for the unknown \(P_{\theta_0}\), and imitate the original experiment:

1. Fit the observed data: \(\hat\theta=h(Y_{\mathrm{obs}}).\)

2. Simulate many new datasets from the fitted model, using the same design and sample size as the original data: \(Y^{*(b)}\sim P_{\hat\theta}.\)

3. Refit the same estimator to every simulated dataset: \(\hat\theta^{*(b)}=h\left(Y^{*(b)}\right).\)

4. Use the distribution of the bootstrap estimates \(\hat\theta^{*(1)},\ldots,\hat\theta^{*(B)}\) to approximate the unknown sampling distribution of \(\hat\theta\).

A concrete example: bootstrapping the coin flip

Go back to Dataset 3 from a few slides ago: H H H H T H H T H H, 8 heads out of 10, so \(\hat p = 0.8\).

In real life we would stop here — we don’t get to see \(p=0.6\) again. The parametric bootstrap asks: if \(\hat p=0.8\) were the truth, how much would \(\hat p\) bounce around across repeated 10-flip experiments?

1. Fit: \(\hat p = h(Y_{\text{obs}}) = 0.8\).
2. Simulate \(B=3000\) new datasets of 10 flips from \(\operatorname{Bernoulli}(0.8)\).
3. Refit \(\hat p^{*(b)}\) on each — here, just the new proportion of heads.
4. Look at the distribution of the \(\hat p^{*(b)}\)’s.

Applying this to SMM

Let \(R\) be the number of simulations used to evaluate the SMM objective, and let \(B\) be the number of bootstrap datasets.

First, fit the observed data:

\[ \hat m_R(\theta) = \frac{1}{R}\sum_{r=1}^R s\left(Y_\theta^{(r)}\right), \qquad Y_\theta^{(r)}\sim P_\theta, \]

\[ \hat\theta = \arg\min_\theta \left[s(Y_{\mathrm{obs}})-\hat m_R(\theta)\right]^\top \hat W \left[s(Y_{\mathrm{obs}})-\hat m_R(\theta)\right]. \]

Then, for each bootstrap replicate \(b=1,\ldots,B\):

1. Simulate a new dataset from the fitted model:

   \[ Y^{*(b)}\sim P_{\hat\theta}. \]

2. Compute its summaries:

   \[ s^{*(b)}=s\left(Y^{*(b)}\right). \]

3. For each candidate \(\theta\), simulate from \(P_\theta\) and estimate its expected summaries:

   \[ \hat m_R^{*(b)}(\theta) = \frac{1}{R}\sum_{r=1}^R s\left(Y_\theta^{*(b,r)}\right), \qquad Y_\theta^{*(b,r)}\sim P_\theta. \]

4. Refit the SMM estimator:

   \[ \hat\theta^{*(b)} = \arg\min_\theta \left[s^{*(b)}-\hat m_R^{*(b)}(\theta)\right]^\top W^{*(b)} \left[s^{*(b)}-\hat m_R^{*(b)}(\theta)\right]. \]

What do we do once we have this bootstrap distribution?

Suppose we have bootstrap estimates

\[ \hat\theta^{*(1)},\ldots,\hat\theta^{*(B)}. \]

We can:

• Plot them. For a vector of parameters, we plot the estimates together to the extent possible to see correlations among the estimates.

• Estimate the standard error. This is just the standard deviation of the bootstrap estimates.

• Look for bias. Compare the center of the bootstrap distribution with the original estimate.

• Construct a simple confidence interval. The 2.5th and 97.5th percentiles of the bootstrap estimates give a 95% percentile interval.

Section 7: Worked demonstration

We return to the Normal location-scale model from Lecture 5:

\[ Y_i\sim N(\mu,\sigma^2), \qquad \theta=(\mu,\sigma). \]

We will fit it using the identifiable summary pair

\[ s(Y)= \begin{pmatrix} \bar Y\\ \hat q_{.75}(Y) \end{pmatrix}. \]

SMM is not the sensible way to fit a Normal model. We are using this model because we understand it well enough to check whether the bootstrap is doing what we intend.

True values (not known for real analyses, but known for pedagogical reasons here) \(\mu_0=5\), \(\sigma_0=2\), and \(n=40\).

The estimator we will use

For a candidate \(\theta=(\mu,\sigma)\), simulate \(R\) datasets of size \(n=40\):

\[ Y_\theta^{(1)},\ldots,Y_\theta^{(R)} \sim N(\mu,\sigma^2), \]

and estimate the expected summaries:

\[ \hat m_R(\theta) = \frac{1}{R}\sum_{r=1}^R s\left(Y_\theta^{(r)}\right). \]

Then minimize

\[ Q_Y(\theta) = \left[s(Y)-\hat m_R(\theta)\right]^\top W \left[s(Y)-\hat m_R(\theta)\right]. \]

Our estimator is the complete numerical procedure

\[ h(Y)=\arg\min_{\mu,\sigma}Q_Y(\mu,\sigma). \]

Step 0: The oracle sampling distribution (not part of the bootstrap!)

This step only exists because we’re in a teaching simulation and happen to know \(\theta_0=(\mu_0,\sigma_0)=(5,2)\). It is not something the parametric bootstrap does, and not something you could do with real data.

Repeat the complete experiment many times:

1. Simulate a dataset: \(Y^{(a)}\sim P_{\theta_0}\).
2. Fit it using the simulation-based SMM estimator: \(\hat\theta^{(a)}=h\left(Y^{(a)}\right)\).
3. Record the estimate: \(\hat\theta^{(a)}\).

The distribution of these estimates, \(\hat\theta^{(1)},\ldots,\hat\theta^{(A)}\), is the sampling distribution of \(\hat\theta\) — the one we would like to know in a real analysis, but ordinarily cannot calculate, because it requires knowing \(\theta_0\).

Step 1: Fit one observed dataset

This is where a real analysis actually starts.

We start by simulating the “observed” data: \(Y_{\mathrm{obs}}\sim P_{\theta_0}\). If you were doing a real analysis, the observed data would be given to you; you wouldn’t simulate it.

Apply the same simulation-based SMM procedure: \(\hat\theta=h(Y_{\mathrm{obs}})=(\hat\mu,\hat\sigma)\).

In this case, imagine we get

\[ \hat\mu=4.918, \qquad \hat\sigma=2.135. \]

Step 2: Construct the bootstrap distribution

For each bootstrap replicate \(b=1,\ldots,B\):

1. Simulate a new dataset from the fitted model: \(Y^{*(b)}\sim P_{\hat\theta}\).
2. Refit it using exactly the same SMM procedure: \(\hat\theta^{*(b)}=h\left(Y^{*(b)}\right)\).
3. Record the estimate: \(\hat\theta^{*(b)}\).

The distribution of \(\hat\theta^{*(1)},\ldots,\hat\theta^{*(B)}\) is the bootstrap distribution.

Comparing the oracle and bootstrap distributions

The two panels show

\[ \underbrace{\hat\theta^{(1)},\ldots,\hat\theta^{(A)}}_{\text{oracle sampling distribution: repeat at }\theta_0} \qquad\text{and}\qquad \underbrace{\hat\theta^{*(1)},\ldots,\hat\theta^{*(B)}}_{\text{bootstrap distribution: repeat at }\hat\theta}. \]

What we calculate from the bootstrap fits

From

\[ \hat\theta^{*(1)},\ldots,\hat\theta^{*(B)}, \]

calculate

\[ \widehat{\operatorname{SE}}_{\mathrm{boot}}(\hat\theta_j) = \operatorname{sd} \left( \hat\theta_j^{*(1)},\ldots,\hat\theta_j^{*(B)} \right), \]

and

\[ \widehat{\operatorname{bias}}_{\mathrm{boot},j} = \operatorname{mean} \left( \hat\theta_j^{*(1)},\ldots,\hat\theta_j^{*(B)} \right) -\hat\theta_j. \]

A simple 95% percentile interval uses the 2.5th and 97.5th percentiles of the bootstrap estimates.

95% bootstrap percentile interval: \(\hat\mu\in[4.23,5.56]\), \(\hat\sigma\in[1.33,3.08]\).

The oracle column is available only because we created the data ourselves. In practice, we would have only the bootstrap column — and here it lands close to the oracle’s, which is the whole point: a genuinely re-optimized bootstrap, with no closed-form shortcuts anywhere in the loop, still recovers essentially the same answer the impossible oracle experiment would have given.

The summary choice still matters

Now change only one part of the estimator:

\[ s_{.75}(Y)= \begin{pmatrix} \bar Y\\ \hat q_{.75}(Y) \end{pmatrix} \qquad\longrightarrow\qquad s_{.55}(Y)= \begin{pmatrix} \bar Y\\ \hat q_{.55}(Y) \end{pmatrix}. \]

For each choice, repeat the same procedure:

1. simulate the summaries at candidate parameter values;
2. minimize the SMM objective;
3. repeat the entire fit across many datasets.

Weak identification reappears

We saw before that the 55th percentile changes very little as \(\sigma\) changes, relative to its sampling variability.

Consequently, when we use the 55th percentile:

• the objective is much flatter in the \(\sigma\) direction;
• estimates of \(\sigma\) vary much more across datasets;
• optimization becomes less stable; and
• constrained estimates may accumulate at the boundary of the parameter space.

Weak identification is visible both as large statistical uncertainty and optimization difficulties.

What we saw

• Repeated simulation + fitting at the known truth produced the oracle sampling distribution.
• Repeated simulation + fitting at the fitted value produced the parametric bootstrap distribution.
• Comparing the two showed what the bootstrap is trying to approximate/showed it does ok.
• Changing the summaries changed the sampling distribution, even though the generative model and sample size stayed the same.