Week 2 Exercises

Complete solutions open Friday, September 4 at 12:00 PM

Do all the exercises below. The Extra Practice problems at the end are optional.

Maximum Likelihood

Exercise 1 A Function of What?

For the Bernoulli model, we wrote the likelihood as $L(\pi) = \pi^{k}(1 - \pi)^{(N - k)}$, where $k$ is the number of successes in $N$ trials. In this expression, which symbol is the variable, and which symbols are fixed? In one sentence describe the question that $L(\pi)$ answers.

Hint 1

In the pmf $f(x; \pi)$, we treat $\pi$ as fixed and $x$ as the variable. What changed when we relabeled $f(x; \pi)$ as $L(\pi)$?

Hint 2

Before the notes ever write $L(\pi)$, the Bernoulli section evaluates the same fixed toothpaste-cap data at several candidate values of $\pi$ (0.1, 0.45, 0.55, and 0.9), one at a time. Reread that paragraph.

Exercise 2 Why the Log? Predict, Then Check

Suppose you observe 2,000 Bernoulli trials with $\pi = 0.5$, so each observation has probability exactly 0.5. The likelihood is the product of these 2,000 probabilities.

Predict what the first line below returns before you run it, and write down your reasoning.

set.seed(42)
y <- rbinom(2000, 1, 0.5)

prod(dbinom(y, 1, 0.5))             # the likelihood, computed directly
sum(dbinom(y, 1, 0.5, log = TRUE))  # the log-likelihood

Now run both lines. Explain what happened, why it’s expected, and why it matters for ML. The log fixes this computational problem, but it also helps with the calculus. How?

Hint 1

Each of the 2,000 factors equals 0.5, so the product is $0.5^{2000}$. Figure out roughly how big that is as a power of ten. Compare that to the smallest positive number R can store (i.e., run .Machine$double.xmin).

A rule of thumb: $2^{10} = 1024 \approx 10^3$. Thus, ten halvings cost about three orders of magnitude.

Hint 2

The notes give two separate reasons for taking the log (i.e., where they define $L(\pi)$ in the Bernoulli example). One reason is what you just watched happen in R. The other reason is about calculus.

Exercise 3 What a Grid Can Miss

Suppose a colleague estimates a Bernoulli $\pi$ by computing the log-likelihood at the eleven values $0, 0.1, 0.2, \ldots, 1$ and reporting the best of the eleven as “the ML estimate.” In one or two sentences, in what sense is that answer wrong, and how would you tighten it? (If you’d like to actually run a grid search, that’s Extra Practice 1.)

Exercise 4 From Derivative to Estimator

Last week, you found $\frac{d\ell}{d\pi}$ for $\ell(\pi) = S \log(\pi) + (N - S)\log(1 - \pi)$. Set the derivative equal to zero and solve for $\hat{\pi}$.

Exercise 5 Reading a Contour Plot

The plot below shows the log-likelihood for the notes’ beta-distribution example. We have 100 observations (simulated from a beta distribution). I computed the log-likelihood $\log L(\alpha, \beta)$ over a grid of parameter pairs.

Roughly where is the ML estimate?
What does a single labeled curve on this plot mean?
For the Bernoulli model, we found $\hat{\pi}$ with pencil and paper. Why can’t we do that here?

For part (a)

Each labeled curve gives the log-likelihood value along it; higher numbers sit closer to the peak. Find the region where the labeled values are highest and read off its $(\alpha, \beta)$ coordinates.

For part (b)

The plot shows the height of $\log L(\alpha, \beta)$. Along one curve, what is constant?

For part (c)

Look back at what stopped us in the notes’ beta example.

Exercise 6 Diagnose the Error

A classmate wants the ML estimate of a Poisson rate $\lambda$ for the data below, and their code runs without an error or a warning:

y <- c(12, 7, 9, 12, 10)

ll_pois <- function(par, y) {
  sum(dpois(y, lambda = par, log = TRUE))
}

est <- optim(par = 1,
             fn = ll_pois,
             y = y,
             method = "Brent",
             lower = 0.001, upper = 100)
est$par

[1] 100

But the estimate is absurd. The average of y is 10, but $\hat{\lambda}$ comes back as 100. What went wrong? What’s the fix? What is the category of the mistake: a typo, an algebra slip, or something else?

Hint

When you hand optim() a function what does it do by default: maximize or minimize? See the Details section of ?optim.

Exercise 7 Trust It or Not?

Your optim() call returns $convergence = 1, and $par equals the upper bound you supplied. Do you trust the estimate? What do you check or change next?

Hint

?optim’s Value section documents every code $convergence can return, not only 0. Look up what 1 means there before you decide anything else.

Exercise 8 An Unfamiliar Density

Suppose a random variable has pdf $f(x; \theta) = \theta x^{\theta - 1}$ for $0 < x < 1$ and $\theta > 0$, and suppose we collect $N$ iid observations $x = \{x_1, x_2, \ldots, x_N\}$.¹ Find the ML estimator of $\theta$.

Hint 1

What steps did we follow for the Bernoulli and the Poisson? Start the same way here by writing $L(\theta)$ as a product over the observations.

Hint 2

Take the log and bring the exponent down. You should end up with two terms: one involving $\log \theta$ and one involving $\sum \log x_i$.

Self-check

Your formula and optim() must agree. For x <- c(0.2, 0.5, 0.9), evaluate your $\hat{\theta}$ formula, then run the code below and compare.

x <- c(0.2, 0.5, 0.9)

ll_fn <- function(theta, x) {
  sum(log(theta) + (theta - 1)*log(x))
}

est <- optim(par = 1,
             fn = ll_fn,
             x = x,
             control = list(fnscale = -1),
             method = "Brent",
             lower = 0, upper = 100)
est$par

The Invariance Property

Exercise 9 One Line of Invariance

Let $\hat{\theta}$ be the ML estimate of $\theta$. Suppose we are interested in $\psi = \theta^2$. What is the ML estimate of $\psi$?

Hint

This week’s notes has a chapter built entirely turning an ML estimate of one quantity into an ML estimate of a function of it, without writing down a new likelihood.

Exercise 10 The Exponential Model

The exponential distribution has pdf $f(t; \lambda) = \lambda e^{-\lambda t}$ for $t \geq 0$ and $\lambda > 0$. We might use the exponential distribution to model durations. A duration is the time until an event, such as the end of a coalition government.

Suppose we collect $N$ iid durations $t = \{t_1, t_2, \ldots, t_N\}$. Find the ML estimator of $\lambda$.
The parameter $\lambda$ is called the rate. The mean of the exponential distribution is $\frac{1}{\lambda}$.² Give the ML estimator of the mean.

For part (a)

This is very similar to the derivation of the ML estimate of $\lambda$ for the Poisson distribution, with one important difference. Be careful with your algebra!

Self-check

# simulate durations from a known rate and see whether your formulas
# recover it
set.seed(1)
lambda_true <- 2
t <- rexp(500, rate = lambda_true)

# fill in your part (a) estimator, as a function of t
lambda_hat <- NA

# a numerical check: maximize the exponential log-likelihood directly
ll_fn <- function(lambda, t) sum(dexp(t, rate = lambda, log = TRUE))
lambda_hat_numeric <- optim(par = 1, fn = ll_fn, t = t,
                             method = "Brent", lower = 1e-6, upper = 20,
                             control = list(fnscale = -1))$par

# these two should agree closely
lambda_hat
lambda_hat_numeric

# fill in your part (b) estimator, as a function of lambda_hat
mean_hat <- NA

# an independent numerical check: reparameterize the same likelihood in
# terms of the mean and maximize over that instead
ll_fn_mean <- function(mu, t) sum(dexp(t, rate = 1 / mu, log = TRUE))
mean_hat_numeric <- optim(par = 1, fn = ll_fn_mean, t = t,
                           method = "Brent", lower = 1e-6, upper = 20,
                           control = list(fnscale = -1))$par

# these two should agree closely
mean_hat
mean_hat_numeric

Exercise 11 Four Models, One Number

Suppose you observe the binary outcome y <- c(0, 1, 0, 1, 1, 1, 0) and you want to estimate the mean of the distribution that generated it. Naturally, you’d model these data as Bernoulli and estimate $\pi$. But suppose that instead you model these 0s and 1s with a normal distribution. What is your ML estimate of the mean? What if you use the Poisson? The exponential?

Complete the table below: for each model, the ML estimates of the parameters, then the estimate of the mean via the invariance property. What do you notice? Explain why it happens.

Distribution	Parameter(s)	ML estimates	Estimate of the mean
Bernoulli	$\pi$	??	??
Normal	$\mu$, $\sigma^2$	??	??
Poisson	$\lambda$	??	??
Exponential	$\lambda$	??	??

Hint 1

For each model, what is $E(Y)$ in terms of that model’s parameter? (Three of the four means are the parameter. The fourth is a transformation of the parameter.)

Hint 2

Estimate each model’s parameters by ML (all four estimators appear in the notes or in Exercise 10). If needed, transform the estimates with the invariance property.

Exercise 12 N or N − 1?

R’s var(y) is not the ML estimate of $\sigma^2$ for a normal model. What’s the difference between the two? Which one is bigger? And does the difference matter at $N = 10{,}000$? Check your answers numerically on any vector you like.

Hint

The notes derive the ML estimator for $\sigma^2$ by solving $\frac{d \log L}{d \sigma^2} = 0$. Find that derivation. What is the denominator? Compare it to the denominator for var(), which is explained deep in the Details of ?var.

Self-check

For the numeric check, build the sum of squared deviations from the mean. Then use the two alternative denominators. Compare each result to var().

# pick any vector, e.g.:
y <- rnorm(1000)

var(y)

# your ML formula for sigma^2, using y directly (not from var(y))

Exercise 13 The SD the Model Implies

Holland (2015) measures the number of enforcement operations (i.e., a “count”) against street vendors across districts in three Latin American cities. The data are in crdata::holland2015.³ The notes estimated the Poisson rate $\lambda$ separately for each city. Now let’s estimate the standard deviation of the number of operations (rather than the rate or mean).

What does the Poisson model assume about the SD of the data, in terms of $\lambda$?
Use the invariance property to estimate the SD of the operations counts in each city from $\hat{\lambda}$.
Compute the sample SD in each city and compare. What do you find, and what does it tell you about the model?

For parts (a) and (b)

For a Poisson distribution, how is the variance related to the mean? Write the SD as a function of $\lambda$, then estimate it by plugging in $\hat{\lambda}$.

The distributions appendix has the answer.

For part (c)

You already grouped the operations counts by city to compute $\hat{\lambda}$ in part (b). Apply that same grouping here, and compute each city’s sample SD with sd().

When the Recipe Fails

Exercise 14 The Discrete Uniform

Suppose a discrete uniform distribution on $\{0, 1, \ldots, K\}$, with pmf $f(x; K) = \frac{1}{K + 1}$ for $x \in \{0, 1, \ldots, K\}$. Suppose we observe a sample of size 3: 276, 159, and 912.

Find the ML estimate of $K$. Hint: The likelihood is discontinuous in $K$, so differentiate-and-solve will mislead you. But the maximum is immediately apparent once you write out the likelihood.
Find the method of moments estimate of $K$. Hint: The mean of this distribution is $\tfrac{K}{2}$. Set the sample mean equal to the model mean and solve.
For these data, the method of moments estimate falls below the largest observation. This is unsatisfying because it’s an estimate of $K$ that the data themselves rule out. How bad can this get? Construct a three-observation dataset that makes the ratio $\hat{K}_{MM} / \max(x)$ as small as you can. What is the smallest value the ratio can take with three observations?
In a sentence or two, say what these two failures teach us about choosing an estimator.

For part (c)

The ratio you’re minimizing is $2 \cdot \text{avg}(x) / \max(x)$ (i.e., part (b)’s formula divided by the largest observation). With $\max(x)$ fixed, which direction should you push the other two observations to shrink the average? Work that out before you pick numbers.

For part (c), a second nudge

Fact: this distribution never produces a value below 0, so that’s the floor for any observation you choose. Combined with part (b)’s formula, how far can you push $\text{avg}(x)$ down without changing $\max(x)$?

Self-check

Plug your candidate dataset from part (c) into the code below and confirm the ratio matches what you computed by hand.

x <- c(a, b, c)               # your three-observation dataset from part (c)
K_hat_mm <- 2 * mean(x)       # method of moments estimate, from part (b)
ratio <- K_hat_mm / max(x)
ratio

Checking the Model

Exercise 15 The Remaining Wait

The exponential distribution from Exercise 10 is our first model for durations. What does the model actually claim about waiting? This turns out to be really interesting and important!

Show that the cdf of the exponential distribution is $F(t; \lambda) = \Pr(T \leq t) = 1 - e^{-\lambda t}$. Then define the survival function $S(t; \lambda) = \Pr(T > t) = 1 - F(t; \lambda)$ and interpret it (i.e., for an input $t$, what does $S$ return)?
Suppose you’ve already waited $s$ units of time without the event occurring. Using the definition of conditional probability and the survival function, find $\Pr(T > t + s \mid T > s)$. Compare your answer to $\Pr(T > t)$. What do you notice? Write it down in one sentence before you start part (c).
Now check your sentence against a simulation. The code below draws 10,000 durations and computes, for the draws that lasted past time 2, how much longer they lasted.
```
set.seed(123)
durations <- rexp(10000, rate = 1)
remaining <- durations[durations > 2] - 2
```
Compare the distribution of remaining to the distribution of the original durations (plotting both ECDFs in one panel works well). Do the plots confirm or refute your sentence from (b)?
Name one political process that your finding might describe well, and one it clearly can’t. Say why.

For part (a)

The cdf accumulates the density: $F(t; \lambda) = \int_0^t \lambda e^{-\lambda u}\, du$. Use integration rules from Week 1.

For part (b)

What is the event “$T > t + s$ and $T > s$”? Notice that one of the two conditions implies the other.

For part (b), a second nudge

Write the conditional probability as a ratio of survival functions and substitute $S(t; \lambda) = e^{-\lambda t}$.

Exercise 16 Herron’s Hockey Data

Herron’s hockey data set records the time between hits⁴ in the regulation periods of all 82 regular-season games for the Chicago Blackhawks.

Before looking at the data, do you expect an exponential model to fit the times between hits? Write down your prediction and your reasoning. Notice that your derivation in Exercise 15 can’t settle this one. Whether hits actually behave this way is an empirical question. It could be that hits cluster (a hit sparks retaliation). Or maybe a long lull makes the next hit imminent.

Model seconds_btw_hits as exponential, estimate the rate and the mean, and use the predictive distribution to evaluate the fit. Were you right?

# load data directly from the web
hockey <- read_csv("https://gist.githubusercontent.com/carlislerainey/0bc3018cd2377022fd045e1c932110a2/raw/fd6dcc28a7c0df456d779e9f0d7a82f15b9b5844/herron-hockey.csv")

# quick look
glimpse(hockey)

Rows: 1,175
Columns: 5
$ game_id                  <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, …
$ period_id                <dbl> 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 1, 1, 1, 1, 2, …
$ time_of_hit              <chr> "25S", "2M 8S", "2M 58S", "3M 8S", "4M 13S", …
$ seconds_played_in_season <dbl> 25, 128, 178, 188, 253, 560, 643, 1407, 1825,…
$ seconds_btw_hits         <dbl> NA, 103, 50, 10, 65, 307, 83, 764, 418, 535, …

Hint 1

The general recipe for the predictive distribution is in the notes and it works for any fitted distribution. You’ve already derived a formula for $\hat{\lambda}$ in Exercise 10. Get that number for the hockey data first. Then simulate.

Hint 2

The notes offer a three-step procedure: estimate $\hat{\lambda}$ (you derived the estimator in Exercise 10), simulate several fake data sets from the fitted exponential distribution, and compare them to the observed data with histograms or ECDFs.

Exercise 17 Heavy Tails: GDP Growth and the Location-Scale t

The location-scale t distribution is a flexible model for data with heavier-than-normal tails (i.e., more extreme observations than a normal model expects). Its pdf is

\[ f(y; \mu, \sigma, \nu) = \frac{\Gamma\left(\frac{\nu+1}{2}\right)}{\Gamma\left(\frac{\nu}{2}\right) \sqrt{\nu \pi} \, \sigma} \left[1 + \frac{1}{\nu} \left( \frac{y - \mu}{\sigma} \right)^2 \right]^{-\frac{\nu+1}{2}}, \]

where $\mu$ is the location parameter (it shifts the distribution, like the normal’s mean), $\sigma$ is the scale parameter (it spreads the distribution, like the normal’s SD), and $\nu$ controls the heaviness of the tails. The $\mu$ and $\sigma$ work similarly to a normal model. Small $\nu$ means heavy tails, and as $\nu \to \infty$ the distribution converges to the normal. (For $\nu > 10$ or so, the $t$ and the normal distributions are hard to tell apart, and $\nu = 1$ is the Cauchy.)

In R, metRology::dt.scaled() computes this density.⁵ Be careful with the argument names in metRology::dt.scaled(). mean is the location $\mu$ (though $\mu$ isn’t always a mean), sd is the scale $\sigma$ (though $\sigma$ isn’t the SD), and df is $\nu$.

We’ll model cross-national GDP growth, which mixes a tight cluster of ordinary economies with a few extreme performers. The code below downloads percentage GDP growth for 2022 from the World Bank.⁶

# load package
library(WDI)

# get annual % gdp growth for 2022
# - "NY.GDP.MKTP.KD.ZG" is percentage gdp growth
#   see https://data.worldbank.org/indicator/NY.GDP.MKTP.KD.ZG
df <- WDI(indicator = "NY.GDP.MKTP.KD.ZG",
          start = 2022,
          end = 2022,
          extra = TRUE) %>%
  # drop aggregates (e.g., European Union)
  filter(region != "Aggregates") %>%
  mutate(pct_gdp_growth = NY.GDP.MKTP.KD.ZG) %>%
  select(country, year, region, pct_gdp_growth) %>%
  na.omit()

# plot histogram
ggplot(df, aes(x = pct_gdp_growth)) +
  geom_histogram(bins = 40)

As a baseline, model pct_gdp_growth with a normal distribution. Estimate $\mu$ and $\sigma$ by ML (both have closed forms), then use the predictive distribution to assess the fit. What does the normal model miss?

Now fit the location-scale t, estimating $\mu$, $\sigma$, and $\nu$ with optim(). Complete the skeleton below by filling in the two blanked lines.

ll_t <- function(par, y) {
  # unpack the parameter vector
  mu    <- par[1]   # location
  sigma <- par[2]   # scale
  nu    <- par[3]   # degrees of freedom (tail heaviness)

  # 1. guard: impossible parameter values should return -Inf
  #    ... your code here ...

  # 2. compute the log-likelihood with metRology::dt.scaled()
  #    ll <- ... your code here ...

  return(ll)
}

Reasonable starting values are c(median(y), sd(y), 10).

Use the predictive distribution to compare the t model to the normal model. What is the ML estimate of $\nu$, and what does it give us?

For part (b)

The density is undefined when the scale or the tail parameter is zero or negative. Decide what your function should return for those impossible values, so that optim() doesn’t consider them.

For part (b), the code

If sigma or nu is not positive, return -Inf; otherwise return sum(dt.scaled(y, mean = mu, sd = sigma, df = nu, log = TRUE)).

Self-check

Your normal-model estimates in part (a) must match the closed forms exactly so that $\hat{\mu}$ = mean(y) and $\hat{\sigma}$ = sqrt(sum((y - mean(y))^2)/length(y)). If your predictive simulations look off, check these two numbers before debugging anything else.

Reflection

Exercise 18 What Should a Model Match?

The exercises above repeatedly compare a fitted model against observed data, and the models often disagree with the data about everything except the mean. Three questions, a few sentences each:

Is it important that a model mimic features of the data beyond the mean?
When might the SD, the tails, or the behavior of waits matter substantively, for the political question itself?
When might those features matter statistically, even if the mean is all you care about?

For part (a)

Exercise 11 earlier on this page fit several different models to the same data. Reread it and compare what varied across the fits to what didn’t, then think about what that contrast implies for when matching the mean is enough.

For part (b)

Two of the exercises above already put this question in a concrete setting: Exercise 17 compared a normal model to a heavier-tailed model for a country’s GDP growth, and Exercise 15 asked whether a wait so far should change your forecast of the wait still to come. Reread whichever is closer to where you’re stuck, and ask what would go wrong for the actual political question, not just for the histogram, if you’d picked the wrong model.

Extra Practice

Extra Practice 1 Bernoulli Grid Search

Suppose you design a Bernoulli experiment that generates successes and failures with an unknown probability $\pi$. You want to estimate $\pi$, so you run the experiment three times and get the outcomes y <- c(0, 1, 0), where 0 is a failure and 1 is a success.

Use ML to estimate $\pi$. But don’t find the maximum analytically or with a hill-climbing algorithm. Instead, use a grid search: use seq() to create ten to twenty candidate values of $\pi$, compute the log-likelihood for each candidate, and locate the candidate that produces the largest log-likelihood. Report your results in a figure. (A table works too, if you’d like both.)

Hint

You’ll need the Bernoulli log-likelihood as a function you can evaluate. That is, write $\log L(\pi) = k \log(\pi) + (N - k)\log(1 - \pi)$, with $k$ successes in $N$ trials, as an R function of $\pi$. Then you can give this R function the whole grid at once to compute the log-likelihoods for the entire grid.

Self-check

Your grid search and the closed-form Bernoulli estimate from your notes should land on nearly the same $\pi$. The grid is just more coarse. Compute both, store them, and check that they differ by less than your grid’s spacing.

# best_grid_pi: the pi your grid search returned as the best candidate
# closed_form_pi: the closed-form Bernoulli estimate from your notes
best_grid_pi <- NA
closed_form_pi <- mean(y)
abs(best_grid_pi - closed_form_pi)

Extra Practice 2 `optim()` Meets the Closed Form

We found that the sample average is the ML estimator of the parameter $\lambda$ of the Poisson distribution. For the data set y <- c(12, 7, 9, 12, 10), show that maximizing the Poisson log-likelihood with optim() produces the same answer as the closed-form solution.

Hint 1

You need a Poisson log-likelihood function before you can call optim(). See the example for the beta distribution in the notes. Then write the Poisson version yourself.

Hint 2

Once you have a log-likelihood function, the beta example in the notes walks through every argument the optim() call needs, one by one, in the numbered list right after the call.

Self-check

Before calling optim(), check that the log-likelihood function itself works as expected. It should return a larger value at the closed-form estimate than at a nearby value that isn’t the estimate.

# evaluate your log-likelihood function at the closed-form estimate
ll_fn(par = mean(y), y = y)

# and at a nearby value that is NOT the estimate
ll_fn(par = mean(y) + 1, y = y)

# if ll_fn is correct, the first number is larger

Extra Practice 3 Corrupted Data

Take the WDI GDP-growth data from Exercise 17 and corrupt them by replacing one observation with a severe data-entry error. Perhaps replace one value of pct_gdp_growth with 10,000. Re-fit the normal and t models. How did each model’s estimate of the location change? Why? Is resistance to corruption a desirable property?

Footnotes

This is Exercise 9 on p. 425 of DeGroot and Schervish’s Probability and Statistics.↩︎
Showing this takes integration by parts; here we’ll take it as known.↩︎
crdata isn’t on CRAN. Install it once with remotes::install_github("carlislerainey/crdata").↩︎
From Wikipedia, a hit is: “Intentionally initiated contact with the player possessing the puck that causes that player to lose possession of the puck.”↩︎
metRology is on CRAN: install.packages("metRology").↩︎
The WDI package (CRAN) fetches the data live, so this chunk needs an internet connection.↩︎

Maximum Likelihood

Exercise 1 A Function of What?

Exercise 2 Why the Log? Predict, Then Check

Exercise 3 What a Grid Can Miss

Exercise 4 From Derivative to Estimator

Exercise 5 Reading a Contour Plot

Exercise 6 Diagnose the Error

Exercise 7 Trust It or Not?

Exercise 8 An Unfamiliar Density

The Invariance Property

Exercise 9 One Line of Invariance

Exercise 10 The Exponential Model

Exercise 11 Four Models, One Number

Exercise 12 N or N − 1?

Exercise 13 The SD the Model Implies

When the Recipe Fails

Exercise 14 The Discrete Uniform

Checking the Model

Exercise 15 The Remaining Wait

Exercise 16 Herron’s Hockey Data

Exercise 17 Heavy Tails: GDP Growth and the Location-Scale t

Reflection

Exercise 18 What Should a Model Match?

Extra Practice

Extra Practice 1 Bernoulli Grid Search

Extra Practice 2 optim() Meets the Closed Form

Extra Practice 3 Corrupted Data

Footnotes

Extra Practice 2 `optim()` Meets the Closed Form