Identifiability: Can we actually estimate parameters?

Today

Last time, we were able to estimate one parameter in the SIR model using one or two hand-chosen summaries.

The goal for today is to scare you — sometimes this doesn’t work, and you don’t know in advance whether it will or not.

A model in which parameters are in principle estimable is an “identifiable” model, and a model in which parameters are not in principle estimable is not.

The plan

First, how to think about named distributions (should have done this earlier)

Then, lab + debrief

Section 1: Named distributions

We’ve been using named distributions (Bernoulli, exponential, binomial) without talking about how you should deal with those you haven’t seen before

We will talk about the mathematical definitions (PDFs and PMFs) later, when it is more directly necessary.

For now, the main things are just to know the following.

If \(Y\) is described as being a random variable drawn from some named distribution, you should understand:

• What values can \(Y\) take?
• What parameters does the distribution have?
• What does each parameter change?
• What are the allowable values and units of those parameters?
• How do I simulate from it in my software?
• Is the software using the same parameterization as the notation?

Applying the checklist to the Normal

\[ Y \sim N(\mu, \sigma^2) \]

• \(Y\) can take any real value.
• \(\mu\) controls location; same units as \(Y\).
• \(\sigma > 0\) controls spread; same units as \(Y\).
• \(\sigma^2\) is the variance — squared units.

Notation and software don’t always match

Mathematical notation for the Normal usually displays the variance, \(\sigma^2\). Both R and NumPy, though, ask for the standard deviation \(\sigma\) directly, not the variance.

rnorm(n, mean = mu, sd = sigma)

rng.normal(loc=mu, scale=sigma, size=n)

Passing a variance where the function expects a standard deviation can lead to incorrect estimates.

Some distributions are functions of simpler ones

Sometimes a random variable you don’t yet understand turns out to be a deterministic function of other random variables you already do understand.

When that’s true, it’s often a good way to understand how the random variable works.

The location-scale construction

For the normal model, the simpler piece is a standard normal distribution: \[ Z \sim N(0,1), \qquad Y = \mu + \sigma Z. \]

If you understand a standard normal random variable, this allows you to understand what a generic normal distribution is: \(\mu\) shifts the standard normal draw; \(\sigma\) stretches it.

Fixed \(\sigma\), varying \(\mu\)

Fixed \(\mu\), varying \(\sigma\)

Many random variables are constructed as deterministic functions of simpler ones

Location-scale is one instance of a pattern you’ll see many times.

For instance, the Binomial distribution is a sum of independent Bernoullis (random variables that take value either 0 or 1), and the Dirichlet distribution comes from drawing independent Gammas (random variables that take positive values) and scaling them so that they sum to 1.

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

\[ G_1,\ldots,G_K \overset{\text{ind}}{\sim} \operatorname{Gamma}(\alpha_k, 1), \qquad S = \sum_k G_k, \qquad \left(\frac{G_1}{S},\ldots,\frac{G_K}{S}\right) \sim \operatorname{Dirichlet}(\alpha_1,\ldots,\alpha_K). \]

Section 2: Estimating badly in a model we understand

Suppose our data model is \(Y_i \sim N(\mu, \sigma^2)\), and we are trying to estimate \(\mu\) and \(\sigma\) using SMM.

Per last time, we need to choose some summaries to use in our estimation procedure.

We will start off by using the mean and the median as our summaries.

Lab: Estimating \(\mu\) and \(\sigma\) in a normal model with SMM

Recall the setup from Lecture 4: choose a candidate \(\theta\), simulate from the model, and average the summaries of the simulated data:

\[ m(\theta) = \frac{1}{B} \sum_{b=1}^{B} s\big(Y_\theta^{(b)}\big), \qquad Y_\theta^{(1)},\ldots,Y_\theta^{(B)} \overset{\text{iid}}{\sim} \text{model at } \theta. \]

The general SMM objective is

\[ \hat\theta = \arg\min_\theta Q(\theta), \qquad Q(\theta) = g(\theta)^\top W\, g(\theta), \qquad g(\theta) = s(Y_{\text{obs}}) - m(\theta), \]

for some weighting matrix \(W\).

Today, \(\theta = (\mu,\sigma)\) and \(s(Y) = (\bar Y, \operatorname{median}(Y))\), and we’ll use \(W = I\) — a valid choice, though not an optimal one:

\[ Q(\mu,\sigma) = \big[s(Y_{\text{obs}}) - m(\mu,\sigma)\big]^\top \big[s(Y_{\text{obs}}) - m(\mu,\sigma)\big] \]

Lab Part A: mean and median

Open the lab notebook. You have an observed sample from a Normal distribution. Using \(s(Y) = (\bar Y, \operatorname{median}(Y))\):

• Evaluate \(Q(\mu,\sigma)\) over a grid of candidate \((\mu,\sigma)\) pairs.
• Plot the resulting surface.
• Find the minimum. Is it a point, or a ridge?
• If it’s a ridge: find two or three different values of \((\mu, \sigma)\) that give (essentially) the same objective value.
• Would running more simulations at each grid point fix this? Why or why not?

Lab Part B: weak identifiability

Now replace the median with the empirical 60th percentile. Recompute the objective surface.

• Does the surface still look flat in \(\sigma\), or has something changed?
• Now try the 75th percentile instead. What changes about the minimum?

Keep your three surfaces (median, 60th, 75th percentile) — we’ll compare them directly in the debrief.

Section 4: Debrief

Was that just simulation noise?

The three surfaces in the lab were simulated — \(B=100\) per grid point, grainy, no two runs identical. That leaves a fair question: is the ridge you saw a real property of the model, or an artifact of a small \(B\)?

For the Normal specifically (and only because it’s simple enough that we have the math to work it out), we know the exact population values your simulations were estimating: \(E_{\mu,\sigma}(\bar Y) = \mu\) and \(q_p(\mu,\sigma) = \mu + \sigma z_p\).

Therefore, we can show what the objective surfaces would have looked like without simulation noise.

The three surfaces, side by side

Three words for what we saw

\[ m(\theta) = E_\theta\{s(Y)\} \]

\(m(\theta)\) is the population version of the simulated average summary from Lecture 4 — what \(\overline{s(Y_\theta)}\) converges to as \(B\to\infty\).

• Identifiable: \(m(\theta) = m(\theta_0)\) only when \(\theta = \theta_0\).
• Not identifiable: several different \(\theta\)’s give the same \(m(\theta)\).
• Weakly identifiable: \(m\) is one-to-one, but it changes very slowly in some direction.

The median panel was not identifiable in \(\sigma\). The 60th-percentile panel was weakly identifiable. The 75th-percentile panel was (comfortably) identifiable.

Section 5: Transfer to the SIR model

In the normal case, you probably had some intuition before we even did the lab that mean + median wasn’t going to work.

You also could have gone through the math that would have told you the same thing.

However, with more complicated models this won’t always be available, and so you have to just go by simulation.

Recall: SMM for the SIR model

Recall the SIR model is

\[ S \longrightarrow I \longrightarrow R \]

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

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

with \(\beta\) the transmission parameter, \(\gamma\) the recovery parameter, and \(R_0=\beta/\gamma\).

Last time, we assumed that we knew \(\gamma\) and tried to estimate \(\beta\) using final size and initial growth rate as summaries.

This time, we will try to use those summaries to estimate both parameters.

The first thing to do is to check whether the summary statistics change with the pairs of parameters.

The final-size surface over \((\beta,\gamma)\)

The white contours (simulated final size) run almost parallel to the dashed red lines — lines of constant \(R_0=\beta/\gamma\).

Two ways a summary map can fail to be one-to-one

• Normal example — insensitivity. \(\sigma\) is simply absent from \(m(\mu,\sigma)\); nothing about it leaks through, no matter which \((\mu,\sigma)\) you’re near.
• SIR example — confounding. Both \(\beta\) and \(\gamma\) enter the final-size map, but only through one combination of them, \(R_0=\beta/\gamma\). Neither parameter is absent — but only their ratio is visible.

Both are non-identifiability in the same formal sense (the map \(m\) is not one-to-one), but they are different failure geometries.

A second summary, a different shape

What do you predict?

Final size’s contours run along lines of constant \(R_0\). Growth rate’s contours run along a visibly different direction. If that’s right:

• Would an SMM objective built from final size alone have a well-localized minimum, or a ridge?
• Same question for growth rate alone.
• What do you expect from combining them — and does it matter whether you standardize each discrepancy by its ordinary variability first, the way we did in Lecture 4?

The actual objective surfaces

Final lab exercise: Competing failure modes

We are reliability engineers, looking at how a particular part fails. Each manufactured unit can fail two ways — wear or corrosion — whichever happens first:

\[ T_{\text{wear}} \sim \operatorname{Exponential}(\theta), \qquad T_{\text{corrosion}} \sim \operatorname{Exponential}(\kappa\theta). \]

Corrosion is \(\kappa\) times as hazardous as wear. \(\theta\) is a shared “how harsh are the operating conditions” factor — temperature, load, usage intensity, whatever — that scales both hazards up or down equally.

Units are inspected only when they fail. Repeated across many units, you only ever find out which cause got there first — not how long the unit lasted.

What is an exponential distribution?

\[ Y \sim \operatorname{Exponential}(\lambda) \]

• \(Y\) can take any nonnegative real value.
• \(\lambda\) controls the rate; units are events per unit time.
• \(E(Y) = 1/\lambda\) — a bigger rate means a shorter wait, on average.

Standard exponential, and varying \(\lambda\)

Lab Part C: competing failure modes

In this setup, we only have one number coming out of our simulator, the wear-failure fraction, so the whole objective is built from that single summary.

Your job is:

• Simulate many units at a candidate \((\theta,\kappa)\) and record the fraction that failed by wear.
• Build an SMM objective comparing that simulated wear-failure fraction to one fixed “observed” batch of units.
• Evaluate it over a grid of candidate \((\theta,\kappa)\) pairs, and plot the surface.
• What do you see? Are the parameters identifiable?

Things get much harder with more parameters

• The number of points in a grid scales poorly with the number of parameters.
• We can’t visualize the full parameter space anymore — only lower-dimensional slices of it.
• Identifying \(d\) parameters takes roughly \(d\) independently varying summaries, but hand-built summaries are hard to make independent once \(d\) gets larger.
• We’ll talk about some ideas for fixing these issues later.