Instructor and TA copy: complete solutions, always open
Suppose \(y\) follows a normal model with mean \(\mu\) and SD \(\sigma\). Then \(\operatorname{avg}(y)\) is the best unbiased estimator of \(\mu\). But \(\operatorname{median}(y)\) is also an unbiased estimator of \(\mu\). Use a Monte Carlo simulation to show how \(\operatorname{avg}(y)\) is a better estimator of \(\mu\) than \(\operatorname{median}(y)\).
The average is better because the SE of the average estimator is smaller than the SE of the median estimator.
# parameters
mu <- 0
sigma <- 1
N <- 10 # sample size of y
n_mc_sims <- 10000 # monte carlo repetitions
# containers
mean_hat <- numeric(n_mc_sims)
median_hat <- numeric(n_mc_sims)
# simulate many data sets y; compute and store estimates for each
for (i in 1:n_mc_sims) {
y <- rnorm(N, mean = mu, sd = sigma)
mean_hat[i] <- mean(y)
median_hat[i] <- median(y)
}
sd(mean_hat) # se of avg estimator[1] 0.318158sd(median_hat) # se of median estimator[1] 0.374932For a large number of repetitions, do the following:
Compare the mean and SD of the simulated point estimates to the theoretical mean and SD of the asymptotic sampling distribution. We care, of course, about the actual sampling distribution, but we have to use an asymptotic approximation in practice. Discuss the differences and whether they seem important or unimportant.
Optional: Go beyond the mean and SD. Use a table of quantiles of the simulated point estimates and theoretical distribution and/or plots of the CDFs to compare the simulated (i.e., actual) and theoretical distributions.
For the exponential model, the asymptotic sampling distribution of \(\hat{\lambda}\) is \(\hat{\lambda} \overset{a}{\sim} N\left(\lambda,\; \frac{\lambda^2}{N}\right)\).
Solution
# load packages
library(tinytable)
options(tinytable_tt_digits = 3)
options(tinytable_latex_placement = "H")
# set seed for reproducibility
set.seed(1234)
# parameters
lambda <- 2
N <- 50 # sample size
R <- 10000 # number of mc repetitions
# loop
lambda_hat <- numeric(R) # container to store results
for (i in 1:R) {
y <- rexp(N, rate = lambda) # simulate data set
lambda_hat[i] <- 1/mean(y) # compute and store estimate
}
# empirical summaries
emp_mean <- mean(lambda_hat)
emp_sd <- sd(lambda_hat)
# asymptotic
asy_mean <- lambda
asy_sd <- sqrt(lambda^2 / N)
# data frame of results
res <- data.frame(
Statistic = c("Mean", "SD"),
Empirical = c(emp_mean, emp_sd),
Asymptotic = c(asy_mean, asy_sd)
)
# print table
tt(
res,
digits = 3,
align = "lrr",
) | Statistic | Empirical | Asymptotic |
|---|---|---|
| Mean | 2.04 | 2 |
| SD | 0.293 | 0.283 |
# quantile comparison table
probs <- c(0.01, 0.05, 0.10, 0.25, 0.50, 0.75, 0.90, 0.95, 0.99)
emp_q <- quantile(lambda_hat, probs = probs)
asy_q <- qnorm(probs, mean = asy_mean, sd = asy_sd)
# combine into data frame
quant <- data.frame(
p = probs,
empirical = emp_q,
asymptotic = asy_q,
difference = emp_q - asy_q
)
# make table
tt(
quant,
digits = 3, # 3 decimals for numbers
align = "lrr",
)| p | empirical | asymptotic | difference |
|---|---|---|---|
| 0.01 | 1.47 | 1.34 | 0.1302 |
| 0.05 | 1.61 | 1.53 | 0.073 |
| 0.1 | 1.69 | 1.64 | 0.0503 |
| 0.25 | 1.83 | 1.81 | 0.0242 |
| 0.5 | 2.01 | 2 | 0.013 |
| 0.75 | 2.22 | 2.19 | 0.028 |
| 0.9 | 2.43 | 2.36 | 0.0628 |
| 0.95 | 2.56 | 2.47 | 0.0995 |
| 0.99 | 2.83 | 2.66 | 0.176 |
# load packages
library(tidyverse)
# plot ecdf and asymptotic distribution
gg_data <- data.frame(lambda_hat = lambda_hat)
ggplot(gg_data, aes(x = lambda_hat)) +
stat_ecdf() +
stat_function(fun = function(x) pnorm(x,
mean = asy_mean,
sd = asy_sd),
color = "red") +
labs(
title = "ECDF of Simulated ML Estimates vs Asymptotic Sampling Distribution",
x = expression(hat(lambda)),
y = "Cumulative Probability"
)
Discussion
Suppose a discrete uniform distribution from 0 to \(K\). The pmf is \(f(x; K) = \frac{1}{K + 1}, \quad x \in \{0, 1, \ldots, K\}\). The ML estimator is \(\hat K = \max(y)\). The method of moments estimator is \(\hat K = 2 \times \operatorname{avg}(y)\). Use a Monte Carlo simulation to evaluate bias and variance of each estimator for a small (\(N = 3\)), medium (\(N = 25\)), and large sample size (\(N = 1,000\)). Which estimator do you recommend?
Note: This model violates several of the regularity conditions, so the usual asymptotic guarantees do not hold.
Solution
# load packages
library(tinytable)
library(dplyr)
options(tinytable_tt_digits = 3)
options(tinytable_latex_placement = "H")
# set seed for reproducibility
set.seed(1234)
# parameters
K <- 100 # true K
Ns <- c(3, 25, 1000) # sample sizes
R <- 10000 # number of mc repetitions
# container: list of data frames
res_list <- list()
# loop over sample sizes
for (N in Ns) {
K_hat_mle <- numeric(R) # container for MLE
K_hat_mom <- numeric(R) # container for MOM
for (i in 1:R) {
y <- sample(0:K, N, replace = TRUE) # simulate data set
K_hat_mle[i] <- max(y) # MLE
K_hat_mom[i] <- 2 * mean(y) # MOM
}
# empirical summaries
emp_mean_mle <- mean(K_hat_mle)
emp_sd_mle <- sd(K_hat_mle)
emp_mean_mom <- mean(K_hat_mom)
emp_sd_mom <- sd(K_hat_mom)
# store results for this N
res_list[[as.character(N)]] <- data.frame(
N = N,
Estimator = c("MLE", "MOM"),
Mean = c(emp_mean_mle, emp_mean_mom),
SD = c(emp_sd_mle, emp_sd_mom),
Bias = c(emp_mean_mle - K, emp_mean_mom - K),
MSE = c(mean((K_hat_mle - K)^2), mean((K_hat_mom - K)^2))
)
}
# combine all into one data frame
res <- bind_rows(res_list)
# print table
tt(
res,
digits = 5,
align = "lrrrrr",
)| N | Estimator | Mean | SD | Bias | MSE |
|---|---|---|---|---|---|
| 3 | MLE | 75.055 | 19.7364 | -24.9454 | 1011.7588 |
| 3 | MOM | 99.952 | 33.59 | -0.048467 | 1128.1803 |
| 25 | MLE | 96.665 | 3.6737 | -3.3353 | 24.6191 |
| 25 | MOM | 100.191 | 11.7466 | 0.191144 | 138.0057 |
| 1000 | MLE | 100 | 0 | 0 | 0 |
| 1000 | MOM | 99.993 | 1.8362 | -0.007073 | 3.3712 |
Notice that the ML estimator is severely biased downward in small samples. But it also has a much smaller variance. It’s not clear how to trade off these two. However, it’s worth noting that the severely biased ML estimator has a smaller mean-squared error.
Suppose the usual exponential model. The ML estimate of the rate \(\lambda\) is \(\hat{\lambda} = \frac{1}{\text{avg}(y)}\). Using the invariance property, the ML estimate of the mean \(\frac{1}{\lambda}\) is \(\hat{\mu} = \frac{1}{\hat{\lambda}}\). Show analytically that the ML estimate of the mean is unbiased, while the ML estimate of the rate is biased.
Jensen’s Inequality. If \(g\) is a convex function and \(X\) is a random variable with \(\mathbb E[X]\) finite, then \(\mathbb E\big[g(X)\big] \ge g \big(\mathbb E[X]\big)\), with strict inequality unless \(X\) is degenerate (i.e., constant almost surely).
First show that \(\hat{\mu}\) is unbiased. Then use Jensen’s inequality to show that \(\hat{\lambda}\) is biased upward. You do not need to find the size of the bias, just use Jensen’s inequality to show that the ML estimate of the rate is biased. To use Jensen’s inequality, let \(g(x)=1/x\), and notice that \(g(x)\) is convex on \((0,\infty)\) (since \(g''(x)=2/x^3>0\)).
Solution
1) Unbiasedness of the ML estimate of the mean
\(\hat{\mu} = \text{avg}(y) = \frac{1}{N} \sum_{i = 1}^N y_i\) is the sample mean of iid observations. \(\mathbb E[y_i]=\frac{1}{\lambda}\). Then
\[ \begin{aligned} \mathbb E[\hat\mu] &= \mathbb E \left[ \frac{1}{N} \sum_{i = 1}^N y_i \right] && \text{(definition of $\hat\mu$)}\\ &= \frac{1}{N} \cdot \mathbb E \left[\sum_{i = 1}^N y_i \right] && \text{(linearity)}\\ &= \frac{1}{N} \cdot \left[\sum_{i = 1}^N \mathbb E (y_i) \right] && \text{(sum of expectations)}\\ &= \frac{1}{N} \cdot \left[\sum_{i = 1}^N \frac{1}{\lambda} \right] && \text{(each has mean $1/\lambda$)}\\ &= \frac{1}{N} \cdot \left[N \cdot \frac{1}{\lambda} \right] && \text{(sum $N$ identical terms)}\\ &= \frac{1}{\lambda} && \text{(simplify)}\\ &= \mu && \text{(definition of $\mu$)} \end{aligned} \]
Because \(\mathbb E(\hat\mu) = \mu\), \(\hat{\mu}\) is an unbiased estimate of the mean.
2) Bias of the ML estimate of the rate
\(\hat{\lambda} = \dfrac{1}{\text{avg}(y)} = \left(\dfrac{1}{N}\sum_{i=1}^N y_i\right)^{-1}\). Let \(g(x)=1/x\), and notice that \(g(x)\) is convex on \((0,\infty)\) (since \(g''(x)=2/x^3>0\)). Then
\[ \begin{aligned} \mathbb E[\hat{\lambda}] &= \mathbb E\left[ \frac{1}{\text{avg}(y)} \right] && \text{(definition of $\hat{\lambda}$)}\\ &= \mathbb E\left[g\big(\text{avg}(y)\big)\right] && \text{(definition of $g$ above)}\\ &> g\big(\mathbb E[\text{avg}(y)]\big) && \text{(Jensen's inequality, $g$ convex)}\\ &= \frac{1}{\mathbb E[\text{avg}(y)]} && \text{(apply $g$)}\\ &= \frac{1}{\frac{1}{N}\sum_{i=1}^N \mathbb E(y_i)} && \text{(linearity of expectation)}\\ &= \frac{1}{\frac{1}{N}\sum_{i=1}^N \frac{1}{\lambda}} && \text{(each $y_i$ has mean $1/\lambda$)}\\ &= \frac{1}{\frac{1}{N}\cdot N \cdot \frac{1}{\lambda}} && \text{(sum of $N$ identical terms)}\\ &= \lambda && \text{(simplify)} \end{aligned} \]
Since the inequality is strict for finite \(N\) (the distribution of \(\text{avg}(y)\) is non-degenerate), \(\mathbb E[\hat{\lambda}] > \lambda\). Thus, the ML estimate of the rate is biased upward.
Suppose a Poisson model. The ML estimate of the mean parameter \(\lambda\) is \(\hat\lambda = \text{avg}(y)\). Find an estimate of the SE using the observed Fisher information.
Solution
The log-likelihood is \(\ell(\lambda) = \sum_{i=1}^N \big[y_i \log \lambda - \lambda - \log(y_i!)\big]\).
First, find the second derivative with respect to \(\lambda\).
\[ \begin{aligned} \ell(\lambda) &= \sum_{i=1}^N \big[y_i \log \lambda - \lambda - \log(y_i!)\big] && \text{Poisson log-likelihood} \\[6pt] \frac{\partial \ell}{\partial \lambda} &= \sum_{i=1}^N \left( y_i \cdot \frac{1}{\lambda} - 1 - 0 \right) && \text{use } \frac{d}{d\lambda}\log \lambda = \frac{1}{\lambda},\; \frac{d}{d\lambda}(-\lambda) = -1,\; \log(y_i!) \text{ is constant} \\[6pt] &= \sum_{i=1}^N \left( \frac{y_i}{\lambda} - 1 \right) = \frac{1}{\lambda} \sum_{i=1}^N y_i - N && \text{collect terms} \\[10pt] \frac{\partial^2 \ell}{\partial \lambda^2} &= \frac{\partial}{\partial \lambda} \left( \frac{1}{\lambda} \sum_{i=1}^N y_i - N \right) && \text{differentiate again} \\[6pt] &= \sum_{i=1}^N \left( - \frac{y_i}{\lambda^2} \right) - 0 && \text{use } \frac{d}{d\lambda}\left(\frac{1}{\lambda}\right) = -\frac{1}{\lambda^2} \text{, } N \text{ is constant} \\[6pt] &= - \sum_{i=1}^N \frac{y_i}{\lambda^2}. \end{aligned} \]
The observed Fisher information is
\[ \mathcal I_{\text{obs}}(\lambda) = -\frac{\partial^2 \ell}{\partial \lambda^2} = \sum_{i=1}^N \frac{y_i}{\lambda^2}. \]
Asymptotically, \(\operatorname{Var}(\hat\lambda) \approx \mathcal I_{\text{obs}}(\lambda)^{-1}\).1
To find the estimate of the SE, replace \(\lambda\) with \(\hat\lambda\):
\[ \widehat{\operatorname{SE}}(\hat\lambda) = \sqrt{\mathcal I_{\text{obs}}(\hat\lambda)^{-1}} = \sqrt{\frac{\hat\lambda^2}{\sum_{i=1}^N y_i}} = \sqrt{\frac{\hat\lambda^2}{N \hat\lambda}} = \sqrt{\frac{\hat\lambda}{N}} = \sqrt{\frac{\operatorname{avg}(y)}{N}}, \]
since \(\sum_{i=1}^N y_i = N \hat\lambda\).
operationsThe code below loads the number of enforcement operations for each district in Bogota from Holland (2015) from the {crdata} package. See ?crdata::holland2015 for details.
# load holland2015 data set from crdata package
holland2015 <- crdata::holland2015
# get operations in bogota only
ops <- holland2015$operations[holland2015$city == "bogota"]
# print
ops [1] 15 10 15 20 8 4 8 16 4 28 2 8 4 4 4 4 8 3 4A Poisson model.
Task: Using the Poisson estimators from the previous question, find the ML estimate and SE estimate of \(\lambda\) for a Poisson model.
But before moving on, compare the SD of the predictive distribution (no need to simulate, it depends only and directly on \(\lambda\)) and the SD of the data. Why might this be a problem for the SE estimate?
Solution.
First, fit the Poisson model and compute the ML estimate and SE of \(\lambda\) without using optim(), following the observed-information derivation above.
# load packages
library(tidyverse)
library(tinytable)
# load holland2015 data set from crdata package
holland2015 <- crdata::holland2015
# get operations in bogota only
ops <- holland2015$operations[holland2015$city == "bogota"]
# estimate lambda with ml
N <- length(ops)
lambda_hat_pois <- mean(ops); lambda_hat_pois[1] 8.894737se_lambda_hat_pois <- sqrt(lambda_hat_pois / N); se_lambda_hat_pois[1] 0.6842105# compare predictive SD to data SD
sd_pred_pois <- sqrt(lambda_hat_pois) # predictive SD under Poisson
sd_data <- sd(ops) # empirical SD
# summarize results in a small table
res_pois <- data.frame(
Statistic = c("Predictive SD ", "Data SD"),
Value = c(sd_pred_pois, sd_data)
)
tt(res_pois, digits = 3, align = "lr")| Statistic | Value |
|---|---|
| Predictive SD | 2.98 |
| Data SD | 6.94 |
Suppose an exponential model of \(y\). We know that the ML estimate of the rate \(\lambda\) is \(\hat \lambda = \frac{1}{\operatorname{avg}(y)}\). By the invariance property, the ML estimate of the mean \(\mu = \frac{1}{\lambda}\) is \(\hat \mu = \frac{1}{\hat \lambda}\).
For an exponential model with rate \(\lambda\), the density is \(f(y; \lambda) = \lambda e^{-\lambda y}\), so \(\ell(\lambda) = \sum \big(\log \lambda - \lambda y_i\big)\).
Solution
The log-likelihood is \[ \ell(\lambda) \;=\; \sum_{i=1}^N \big[\log \lambda \;-\; \lambda y_i\big]. \]
First, find the second derivative with respect to \(\lambda\).
\[ \begin{aligned} \ell(\lambda) &= \sum_{i=1}^N \big[\log \lambda - \lambda y_i\big] && \text{Exponential log-likelihood} \\[6pt] \frac{\partial \ell}{\partial \lambda} &= \sum_{i=1}^N \left(\frac{1}{\lambda} - y_i\right) && \text{use } \frac{d}{d\lambda}\log \lambda = \frac{1}{\lambda},\; \frac{d}{d\lambda}(-\lambda y_i) = -y_i \\[6pt] &= \frac{N}{\lambda} \;-\; \sum_{i=1}^N y_i && \text{collect terms} \\[10pt] \frac{\partial^2 \ell}{\partial \lambda^2} &= \frac{\partial}{\partial \lambda}\!\left(\frac{N}{\lambda} - \sum_{i=1}^N y_i\right) && \text{differentiate again} \\[6pt] &= -\,\frac{N}{\lambda^2} \;-\; 0 && \text{use } \frac{d}{d\lambda}\!\left(\frac{1}{\lambda}\right) = -\frac{1}{\lambda^2},\; \sum y_i \text{ is constant}. \end{aligned} \]
The observed Fisher information is \[ \mathcal I_{\text{obs}}(\lambda) = -\frac{\partial^2 \ell}{\partial \lambda^2} = \frac{N}{\lambda^2}. \]
To find the estimate of the SE, replace \(\lambda\) with \(\hat\lambda\): \[ \widehat{\operatorname{SE}}(\hat\lambda) = \sqrt{\mathcal I_{\text{obs}}(\hat\lambda)^{-1}} = \sqrt{\frac{\hat\lambda^{2}}{N}} = \frac{\hat\lambda}{\sqrt{N}} = \frac{1}{\sqrt{N}\,\operatorname{avg}(y)}. \]
Now suppose we use the invariance property to obtain \(\hat\mu = \dfrac{1}{\hat\lambda}\). Then we can use the delta method to estimate the SE with \(\tau(\lambda) = \dfrac{1}{\lambda}\).
The first derivative is \(\tau'(\lambda) = -\frac{1}{\lambda^2}\).
Plugging in, we have \[ \begin{aligned} \widehat{\operatorname{Var}}(\hat\mu) &\approx \big[\tau'(\hat\lambda)\big]^2 \cdot \widehat{\operatorname{Var}}(\hat\lambda) \;=\; \left(\frac{1}{\hat\lambda^2}\right)^2 \cdot \frac{\hat\lambda^{2}}{N} \\[6pt] &= \frac{1}{N\,\hat\lambda^{2}}. \end{aligned} \]
And the standard error is \[ \widehat{\mathrm{SE}}(\hat\mu) \;\approx\; \sqrt{\frac{1}{N\,\hat\lambda^{2}}} \;=\; \frac{1}{\sqrt{N}\,\hat\lambda} \;=\; \frac{\operatorname{avg}(y)}{\sqrt{N}}. \]
Suppose you obtain the ML estimate \(\hat\theta=(\hat a,\hat b)=(2,\,1)\) with an estimated covariance matrix of \[ \widehat{\mathrm{Var}}(\hat\theta) = \widehat{\mathrm{Var}}(\hat a,\hat b) = \begin{bmatrix} 1 & \tfrac{1}{2}\\[2pt] \tfrac{1}{2} & 1 \end{bmatrix}. \]
However, you want to estimate the quantity of interest \(q = \tau(a,b)=\dfrac{a}{b}\). Find \(\hat q\) using the invariance property and estimate the SE of \(\hat q\) using the delta method. (Find the gradient analytically.) Confirm your work with R.
Solution
First, find the gradient. For \(\tau(a,b)=a/b\),
\[ \nabla\tau(a,b)= \begin{bmatrix} \frac{\partial}{\partial a}(a/b)\\[6pt] \frac{\partial}{\partial b}(a/b) \end{bmatrix} = \begin{bmatrix} 1/b\\[4pt] -\,a/b^2 \end{bmatrix}. \] Plugging in \((\hat a,\hat b)=(2,1)\), we have \(\nabla\tau(\hat a,\hat b)=\begin{bmatrix} 1\\[2pt] -2 \end{bmatrix}\).
Now do the required matrix multiplication.
\[ \begin{aligned} \widehat{\mathrm{Var}}[\tau(\hat\theta)] &\approx \overbrace{\begin{bmatrix}1 & -2\end{bmatrix}}^{\nabla\tau(\hat\theta)^{\!\top}} \Biggl( \underbrace{\begin{bmatrix} 1 & \tfrac{1}{2}\\[2pt] \tfrac{1}{2} & 1 \end{bmatrix}}_{\widehat{\mathrm{Var}}(\hat\theta)} \underbrace{\begin{bmatrix}1\\ -2\end{bmatrix}}_{\nabla\tau(\hat\theta)} \Biggr) \\[10pt] &= \begin{bmatrix}1 & -2\end{bmatrix} \begin{bmatrix} 1\cdot 1 + \tfrac{1}{2}\cdot(-2)\\[4pt] \tfrac{1}{2}\cdot 1 + 1\cdot(-2) \end{bmatrix} = \begin{bmatrix}1 & -2\end{bmatrix} \begin{bmatrix} 0\\[2pt] -1.5 \end{bmatrix} = 0 + 3 = 3. \end{aligned} \] Thus \(\widehat{\mathrm{Var}}[\widehat{\tau}]=3\) and \(\widehat{\mathrm{SE}}(\widehat{\tau})=\sqrt{3}\).
We can confirm our work with R.
# given values
a_hat <- 2
b_hat <- 1
V_hat <- matrix(c(1, 1/2,
1/2, 1), nrow = 2, byrow = TRUE)
# gradient at (a_hat, b_hat) for tau(a,b) = a/b
grad <- c(1 / b_hat, -a_hat / b_hat^2) # c(1, -2)
# delta method
var_tau_hat <- grad %*% V_hat %*% grad
se_tau_hat <- sqrt(var_tau_hat)
var_tau_hat # 3 [,1]
[1,] 3se_tau_hat # sqrt(3) [,1]
[1,] 1.732051Load a tibble version of the 2024 turnout data from the UF Election Lab into R.
Model the proportion of the voting-eligible population that voted using a beta distribution. Estimate the parameters of the beta distribution \(\alpha\) and \(\beta\) as well as the mean \(\mu = \frac{\alpha}{\alpha + \beta}\) and the SD \(\sigma = \sqrt{\frac{\alpha \beta}{(\alpha + \beta)^2 (\alpha + \beta + 1)}}\). Use a parametric bootstrap to estimate the SE of each.
Solution
# load packages
library(tidyverse)
library(tinytable)
# load data as tibble
turnout_2024 <- tribble(
~state, ~vep_turnout,
"Alabama", 0.5893,
"Alaska", 0.6378,
"Arizona", 0.6360,
"Arkansas", 0.5348,
"California", 0.6206,
"Colorado", 0.7314,
"Connecticut", 0.6708,
"Delaware", 0.6702,
"District of Columbia", 0.6357,
"Florida", 0.6671,
"Georgia", 0.6826,
"Hawaii", 0.5027,
"Idaho", 0.6344,
"Illinois", 0.6325,
"Indiana", 0.5869,
"Iowa", 0.7078,
"Kansas", 0.6318,
"Kentucky", 0.6219,
"Louisiana", 0.6077,
"Maine", 0.7424,
"Maryland", 0.6930,
"Massachusetts", 0.6803,
"Michigan", 0.7464,
"Minnesota", 0.7635,
"Mississippi", 0.5745,
"Missouri", 0.6434,
"Montana", 0.6820,
"Nebraska", 0.6796,
"Nevada", 0.6580,
"New Hampshire", 0.7405,
"New Jersey", 0.6724,
"New Mexico", 0.5957,
"New York", 0.6044,
"North Carolina", 0.7032,
"North Dakota", 0.6310,
"Ohio", 0.6539,
"Oklahoma", 0.5328,
"Oregon", 0.7194,
"Pennsylvania", 0.7143,
"Rhode Island", 0.6328,
"South Carolina", 0.6214,
"South Dakota", 0.6400,
"Tennessee", 0.5761,
"Texas", 0.5657,
"Utah", 0.6415,
"Vermont", 0.7089,
"Virginia", 0.7119,
"Washington", 0.7017,
"West Virginia", 0.5546,
"Wisconsin", 0.7664,
"Wyoming", 0.6128
)
# beta log-likelihood
beta_ll_fn <- function(theta, y) {
alpha <- theta[1]
beta <- theta[2]
ll <- sum(dbeta(y, shape1 = alpha, shape2 = beta, log = TRUE))
return(ll)
}
# create a function to fit the beta model
est_beta <- function(y) {
est <- optim(par = c(3, 3), fn = beta_ll_fn, y = y,
control = list(fnscale = -1),
method = "BFGS") # for >1d problems
if (est$convergence != 0) print("Model did not converge!")
res <- list(est = est$par)
return(res)
}
# ml estimates
ml_est <- est_beta(turnout_2024$vep_turnout)
a <- ml_est$est[1]; a # alpha[1] 39.59784b <- ml_est$est[2]; b # beta[1] 21.29332a / (a + b) # mean; using a and b for readability[1] 0.6503053sqrt( (a * b) / ((a + b)^2 * (a + b + 1)) ) # SD[1] 0.06061623# parametric bootstrap
n_bs <- 100
bs_alpha <-
bs_beta <-
bs_mu <-
bs_sigma <-
numeric(n_bs) # containers for the estimates
for (i in 1:n_bs) {
# bootstrap sample
bs_y <- rbeta(length(turnout_2024$vep_turnout), shape1 = ml_est$est[1], shape2 = ml_est$est[2])
# fit model to bootstrapped sample
bs_est <- est_beta(bs_y)
# store ml estimates of alpha and beta for readability below
bs_a <- bs_est$est[1]
bs_b <- bs_est$est[2]
# compute and store all qis
bs_alpha[i] <- bs_a
bs_beta[i] <- bs_b
bs_mu[i] <- bs_a / (bs_a + bs_b)
bs_sigma[i] <- sqrt( (bs_a * bs_b) / ((bs_a + bs_b)^2 * (bs_a + bs_b + 1)) )
}[1] "Model did not converge!"# bootstrap SEs
se_alpha <- sd(bs_alpha)
se_beta <- sd(bs_beta)
se_mu <- sd(bs_mu)
se_sigma <- sd(bs_sigma)
# tinytable (avoid underscores in visible labels)
tab_df <- tibble::tibble(
Parameter = c("$\\alpha$", "$\\beta$", "$\\mu$", "$\\sigma$"),
`ML estimate` = c(a_hat = a,
b_hat = b,
mu_hat = a / (a + b),
sigma_hat = sqrt( (a * b) / ((a + b)^2 * (a + b + 1)) )),
`Bootstrap SE`= c(se_alpha, se_beta, se_mu, se_sigma)
)
tt(
tab_df,
digits = 3,
caption = "Beta model: ML estimates and parametric-bootstrap standard errors.",
)| Parameter | ML estimate | Bootstrap SE |
|---|---|---|
| $\alpha$ | 39.5978 | 7.99754 |
| $\beta$ | 21.2933 | 4.36365 |
| $\mu$ | 0.6503 | 0.00827 |
| $\sigma$ | 0.0606 | 0.00547 |
Consider a simple Bernoulli model for \(y \in \{0,1\}\) with parameter \(\pi \in [0,1]\). Create a small sample with a rare event: take \(N = 20\) trials with exactly one success.
Part 1. Wald confidence interval.
Compute the ML estimate \(\hat \pi = \operatorname{avg}(y)\), the Wald SE \(\widehat{\mathrm{SE}}(\hat \pi) = \sqrt{\hat \pi(1-\hat \pi)/N}\), and the usual \(95\%\) Wald confidence interval \(\hat \pi \pm 1.96 \cdot\widehat{\mathrm{SE}}(\hat \pi)\). What do you notice about the interval relative to the parameter space \([0,1]\)?
Part 2. Parametric bootstrap confidence interval.
Use a parametric bootstrap to obtain a \(95\%\) percentile confidence interval for \(\pi\) based on the fitted model \(\text{Bernoulli}(\hat \pi)\). Compare the bootstrap interval with the Wald interval. What stands out? Briefly explain.
For the bootstrap, repeatedly simulate \(y^{*(b)}_1,\dots,y^{*(b)}_N \stackrel{iid}{\sim} \text{Bernoulli}(\hat \pi)\), recompute \(\hat \pi^{*(b)} = \operatorname{avg}(y^{*(b)})\), and then take the empirical \(2.5^\text{th}\) and \(97.5^\text{th}\) percentiles of \(\{\hat \pi^{*(b)}\}\).
Solution.
set.seed(1)
# data: N = 20 with exactly one success
N <- 20
y <- c(1, rep(0, N - 1)) # rare event sample
y <- sample(y) # shuffle (irrelevant, but tidy)
# ml estimate and wald 95%ci
pi_hat <- mean(y)
se <- sqrt(pi_hat * (1 - pi_hat) / N)
wald_ci <- pi_hat + c(-1, 1) * 1.96 * se
pi_hat[1] 0.05se[1] 0.04873397wald_ci[1] -0.04551858 0.14551858The Wald interval extends below 0, which is outside the parameter space \([0,1]\).
# parametric bootstrap w/ percentile 95% ci
n_bs <- 2000
bs_est <- numeric(n_bs) # a container for the estimates
for (i in 1:n_bs) {
bs_y <- rbinom(20, size = 1, prob = pi_hat)
bs_est[i] <- mean(bs_y)
}
print(quantile(bs_est, probs = c(0.025, 0.975)), digits = 2) # 95% ci 2.5% 97.5%
0.00 0.15 The Wald interval uses a normal approximation centered at \(\hat \pi\) with variance estimated by Fisher information. With extreme \(\hat \pi\) near zero and small \(N\), this approximation is poor and the lower bound of the CI can dip below 0. In contrast, the parametric bootstrap respects the model’s support by construction (every \(\hat \pi^{*(b)} \in [0,1]\)), so the percentile interval stays within \([0,1]\) and typically provides more reasonable coverage in this setting.
Rather than “approximately equals” (\(\approx\)), it’s more technically correct to write “asymptotically equals.” However, we treat “asymptotically equals” to mean “approximately” in practice. Given the practical orientation of the course, it seems better to write “approximately” here.↩︎