```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(ggplot2)
theme_set(theme_bw())
```

## A. Background and model

We have a dataset that we are willing to model as

$$
Y_i \overset{\text{iid}}{\sim} N(\mu, \sigma^2), \qquad i = 1,\ldots,n,
$$

and we want to use Simulated Method of Moments to estimate $(\mu,\sigma)$
from a chosen pair of summary statistics. Today's question is not "what
are the estimates" — it's "how well do different choices of summary
statistic actually pin down $(\mu,\sigma)$ at all?"

Here is the observed dataset we'll use for the whole lab:

```{r observed-data}
y_obs <- c(
  4.226, 5.571, 5.513, 6.598, 5.302, 2.529, 1.860, 4.149, 6.535, 6.294,
  2.851, 4.269, 5.017, 6.224, 3.615, 6.741, 5.948, 6.383, 2.505, 6.692,
  4.865, 7.873, 3.243, 4.627, 4.275, 4.902, 9.423, 5.388, 7.908, 7.538,
  3.688, 5.061, 6.453, 3.536, 4.859, 4.655, 4.866, 2.118, 2.790, 5.170,
  3.147, 5.565, 5.555, 4.817, 2.832, 10.738, 2.647, 4.595, -0.094, 3.852
)
length(y_obs)
```

As in Lecture 4, an SMM objective compares the observed summary to a
*simulated* average summary, $\overline{s(Y_\theta)}$, built by actually
simulating from the model $B$ times at a candidate $\theta=(\mu,\sigma)$,
computing the summary each time, and averaging.

## B. Provided: the simulator and a plotting function

For a candidate $(\mu,\sigma)$, `simulate_summary_avg()` simulates $B$
datasets of size $n$ from $N(\mu,\sigma^2)$, computes the sample mean
and the sample $p$-th percentile of each one, and returns the average
of each across the $B$ replicates — exactly $\overline{s(Y_\theta)}$
for $s(Y) = (\bar Y, q_p(Y))$.

```{r simulate-summary-avg}
simulate_summary_avg <- function(mu, sigma, n, B, p) {
  means  <- numeric(B)
  quants <- numeric(B)
  for (b in seq_len(B)) {
    y <- rnorm(n, mean = mu, sd = sigma)
    means[b]  <- mean(y)
    quants[b] <- quantile(y, p, type = 7)
  }
  c(mean = mean(means), quant = mean(quants))
}
```

This is the second-most expensive function in the whole course, after
the SIR simulator — every single grid point now costs real computer
time, unlike the instant, closed-form surfaces from the identification
debrief slides. Below is a plotting helper for a completed grid of
$Q(\mu,\sigma)$ values — you don't need to modify it, just call it once
your grid is ready.

```{r plot-surface-fn}
plot_objective_surface <- function(mu_grid, sigma_grid, Q_matrix, title = "") {
  df <- expand.grid(mu = mu_grid, sigma = sigma_grid)
  df$Q <- as.vector(Q_matrix)
  ggplot(df, aes(x = mu, y = sigma, z = Q)) +
    geom_raster(aes(fill = Q)) +
    geom_contour(color = "white", alpha = 0.5, bins = 15) +
    scale_fill_viridis_c(trans = "sqrt", name = "Q(mu,sigma)") +
    coord_cartesian(expand = FALSE) +
    labs(x = expression(mu), y = expression(sigma), title = title)
}
```

## C. Required exercise 1 (Part A): mean and median

Use $s(Y) = (\bar Y, \operatorname{median}(Y))$. The SMM objective is

$$
Q(\mu,\sigma) = \big[\bar y_{\text{obs}} - \overline{s_1(Y_{(\mu,\sigma)})}\big]^2
+ \big[\operatorname{median}(y_{\text{obs}}) - \overline{s_2(Y_{(\mu,\sigma)})}\big]^2,
$$

where the two barred quantities are exactly what
`simulate_summary_avg()` returns. Fill in the blanks below to compute
the observed summaries and evaluate $Q$ over a grid of candidate
$(\mu,\sigma)$ pairs. Notice the grid is much coarser than you might
expect ($25\times 25$), and $B$ is modest (100) — with genuine
simulation, every grid point costs real time, so we can't afford the
effectively-unlimited resolution we had when a formula did the work:

```{r part-a-objective, eval=FALSE}
obs_mean   <- ___
obs_median <- ___

mu_grid    <- seq(0, 10, length.out = 25)
sigma_grid <- seq(0.1, 6, length.out = 25)
B <- 100

Q_median <- matrix(NA, nrow = length(mu_grid), ncol = length(sigma_grid))
for (i in seq_along(mu_grid)) {
  for (j in seq_along(sigma_grid)) {
    sim <- simulate_summary_avg(mu_grid[i], sigma_grid[j],
                                 n = length(y_obs), B = B, p = ___)
    Q_median[i, j] <- (obs_mean - sim["mean"])^2 + (obs_median - sim["quant"])^2
  }
}

plot_objective_surface(mu_grid, sigma_grid, Q_median, title = "mean + median")
```

This will take on the order of ten seconds — that's genuine simulation
time, not a typo, and it's the whole reason last section's grid was
this coarse.

Questions:

1. Describe the shape of the surface. Is the minimum a single point, or
   something else?
2. Pick a value of $\mu$ near the bottom of the valley. Read off the
   objective value at that $\mu$ for three quite different values of
   $\sigma$ (e.g. $\sigma=0.5$, $\sigma=2$, $\sigma=5$). Are they
   exactly equal? Should you expect them to be?

## D. Required exercise 2 (Part B): weak identification

Now replace the median with the **60th percentile**. Adapt your code
from Part A (same `mu_grid`, `sigma_grid`, and `B`):

```{r part-b-p60, eval=FALSE}
obs_p60 <- as.numeric(quantile(y_obs, ___, type = 7))

Q_p60 <- matrix(NA, nrow = length(mu_grid), ncol = length(sigma_grid))
for (i in seq_along(mu_grid)) {
  for (j in seq_along(sigma_grid)) {
    sim <- simulate_summary_avg(mu_grid[i], sigma_grid[j],
                                 n = length(y_obs), B = B, p = ___)
    Q_p60[i, j] <- (obs_mean - sim["mean"])^2 + (obs_p60 - sim["quant"])^2
  }
}

plot_objective_surface(mu_grid, sigma_grid, Q_p60, title = "mean + 60th percentile")
```

Now do the same with the **75th percentile**:

```{r part-b-p75, eval=FALSE}
obs_p75 <- as.numeric(quantile(y_obs, ___, type = 7))

Q_p75 <- matrix(NA, nrow = length(mu_grid), ncol = length(sigma_grid))
for (i in seq_along(mu_grid)) {
  for (j in seq_along(sigma_grid)) {
    sim <- simulate_summary_avg(mu_grid[i], sigma_grid[j],
                                 n = length(y_obs), B = B, p = ___)
    Q_p75[i, j] <- ___
  }
}

plot_objective_surface(mu_grid, sigma_grid, Q_p75, title = "mean + 75th percentile")
```

Question:

1. Compare the median, 60th-percentile, and 75th-percentile surfaces.
   The 60th-percentile surface should look visibly different from the
   flat median surface — but is its minimum as sharply localized as
   the 75th-percentile surface's, or still noticeably more spread out?
   What might this suggest about how precise the estimate is?

## E. Required exercise 3 (Part C): competing failure modes
Parts A and B relied on genuine simulation, one grid point at a time.

**The model.** 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 story you like — that scales *both* hazards
up or down equally. Units are inspected only when they fail, and only
the **cause of failure** is recorded (wear or corrosion) — not how long
the unit lasted. Repeated across many independently manufactured units,
the only thing you ever get is which cause fired first for each one.
(This is the standard *competing risks* setup from reliability
engineering and survival analysis.)


### The simulator

```{r failure-simulator}
## simulate n independent units, return the fraction that failed by wear
wear_fraction <- function(n, theta, kappa) {
  T_wear      <- rexp(n, rate = theta)
  T_corrosion <- rexp(n, rate = kappa * theta)
  mean(T_wear < T_corrosion)
}
```

A whole batch of `n` units is one vectorized call — no loop over grid
points needed, unlike the SIR simulator. Try it a few times at the same
$(\theta,\kappa)$ to see how much the wear-failure fraction bounces
around, and compare two very different $\theta$'s at the same $\kappa$:

```{r failure-try-it}
wear_fraction(n = 200, theta = 5, kappa = 2)
wear_fraction(n = 200, theta = 5, kappa = 2)
wear_fraction(n = 200, theta = 1, kappa = 2)
```

### An observed batch of failures

```{r failure-observed, eval=FALSE}
set.seed(2026)
theta0 <- ___   # pick any true value you like -- you won't need it until the end
kappa0 <- ___   # pick any true value you like -- you won't need it until the end
n_obs  <- 300
obs_wear_fraction <- wear_fraction(n_obs, theta0, kappa0)
obs_wear_fraction
```

(In a real problem you wouldn't know `theta0`/`kappa0` either — setting
them yourself here just lets you check your answer at the end.)

### Build the SMM objective

For a candidate $(\theta,\kappa)$, simulate $B$ units and compare the
simulated wear-failure fraction to `obs_wear_fraction`:

```{r failure-objective, eval=FALSE}
theta_grid <- seq(0.5, 10, length.out = 40)
kappa_grid <- seq(0.2, 5, length.out = 40)
grid <- expand.grid(theta = theta_grid, kappa = kappa_grid)

B <- 500  # simulated units per grid point
grid$Q <- NA
for (i in seq_len(nrow(grid))) {
  sim_wear_fraction <- wear_fraction(___, ___, ___)
  grid$Q[i] <- (obs_wear_fraction - sim_wear_fraction)^2
}

ggplot(grid, aes(x = theta, y = kappa, z = Q)) +
  geom_raster(aes(fill = Q)) +
  geom_contour(color = "white", alpha = 0.5, bins = 15) +
  scale_fill_viridis_c(trans = "sqrt", name = "Q") +
  coord_cartesian(expand = FALSE) +
  labs(x = expression(theta), y = expression(kappa))
```

Questions:

1. Pick a value of $\kappa$ near the bottom of the ridge. Read off $Q$
   at that $\kappa$ for a few very different values of $\theta$. Does
   it change?
2. Go back and check on `theta0` and `kappa0` from the `failure-observed` chunk
   above. Does the ridge pass through $\kappa_0$? Does it single out
   $\theta_0$, or does it run through every $\theta$ you tried?
3. Compare how noisy this surface looks to the Normal surfaces from
   Parts A/B and to the SIR surfaces from class — all three are
   built from simulation, but they don't look equally grainy.
   What affects how noisy an SMM objective surface looks: the number
   of simulations $B$, the number of grid points, how variable the
   underlying summary statistic is, or some combination? (A win/loss
   fraction is a proportion — think about how its variability compares
   to a sample mean's.)
4. Suppose you increased `B` a lot — many more units per grid point.
   Would that ever turn the ridge into a point? Compare your answer
   here to the same question in Part A.
5. In reliability terms: $\kappa$ is the *relative hazard* of corrosion
   versus wear, and it's exactly what you can learn from failure-cause
   records alone. What additional piece of data collection — beyond
   recording *which* cause of failure occurred — would let you also
   recover $\theta$, the absolute failure rate?

## F. Optional extensions

These are separate from the required work above — come back to them if
you have time.

- **Variance instead of a quantile.** Write a version of
  `simulate_summary_avg()` that tracks the sample variance instead of a
  percentile, and build the corresponding objective surface by
  simulation, the same way as Parts A/B. What does the resulting
  surface look like compared to the mean+median surface? (Hint: think
  about which parameter each summary is blind to, and what that
  predicts about the orientation of the ridge, before you plot
  anything.)
- **Three summaries at once.** Build an objective using $s(Y) = (\bar Y,
  \operatorname{median}(Y), q_{0.75}(Y))$ — three summaries, two
  parameters. Does adding the median back in, alongside the 75th
  percentile, change where the minimum is or how sharp it is?
- **Push $B$ around.** Redo the mean+median surface (Part A) with a
  much smaller $B$ (say 20) and a much larger one (say 1000), on a
  coarser grid if 1000 is too slow. Does a small $B$ ever produce a
  surface with an accidental, spurious minimum along the $\sigma$
  direction, rather than a clean ridge?