Efron Hastie CASI exercise 4.3 This is an exercise from Computer Age Statistical Inference by Efron and Hastie, student edition of the book.

Draw a schematic graph of $\dot{l}_{x}(\theta)$ versus $\theta$. Use it to justify (4.25) $\left[\hat{\theta}\dot{=}\theta + \frac{\dot{l}_{x}(\theta)/n}{\ddot{l}_{x}(\theta)/n}\right]$

I am unsure of what is being asked here. The score $\dot{l}_{x}(\theta)$ has been defined in several ways across the chapter. For example, in equation 4.14 as $\frac{\partial}{\partial\theta}\log f_{\theta}(x)$, equation 4.16 with $\dot{l}_{x}(\theta)\sim (0, \mathcal{I}_{\theta})$ (no distribution in here), equation 4.26 with $\dot{l}_{x}(\theta)/n \dot{\sim} \mathcal{N}(0, \mathcal{I}_{\theta}/n)$, and in the specific case of the Gaussian distribution, equation 4.30 gives $\dot{l}_{x}(\theta)=\frac{1}{\sigma^{2}}\sum_{i=1}^{n}(x_{i}-\theta)$.
I made plots for equation 4.26 and equation 4.30 (with a single observation):


Additionally, what is there to justify equation 4.25? I think the point is that if we find the $\hat{\theta}$ so that $\dot{l}_{x}(\theta)=0$, then the estimate is unbiased. If we don't find the exact value, some form of spread on the bias is given in inverse proportion to $\ddot{l}_{x}(\theta)$. Is this interpretation correct? Is it the justification that the book is asking for?
 A: Rather than a schematic graph of the score function $\dot{l}_x(\theta)$ vs $\theta$, I plot $\dot{l}_x(\theta)$ for an assortment of distributions. Then a schematic can be generalized from those cases.
The setup is as follows: For a fixed parameter $\theta$, I take $n = 10$ draws from the distribution and plot the score function $\dot{l}_x(\theta) = \frac{\partial}{\partial\theta}\log L_x(\theta)$, the derivative of the log-likelihood $L_x(\theta)$ with respect to $\theta$. I repeat this "experiment" 20 times for each distribution, so that we can see the sampling variability of $\dot{l}_x(\theta)$.

We can see that:

*

*The score function is approximately linear in a neighborhood of the true parameter $\theta$. (And exactly linear for the Normal model.)

*[In the vertical direction] At the true parameter $\theta$, the score function $\dot{l}_x(\theta)$ varies about zero.

*[In the horizontal direction] The MLE $\hat{\theta}$ (the value of $\theta$ such that score is equal to 0) varies about $\theta$.

This basically describes what the schematic graph should look like. Now let's use these properties of the score function to justify the approximation (4.25):
$$
\begin{aligned}
\hat{\theta} = \theta + \frac{\dot{l}_x(\theta)/n}{-\ddot{l}_x(\theta)/n}
\end{aligned}
$$
For simplicity, let's start with a line $y = a + bx$ which crosses $y = 0$ at $x_0$. From this we can solve for the intercept $a = -bx_0$, so the line becomes $y = b(x - x_0)$ and we can rewrite this as $x_0 = x - y / b$. The first derivative of $y$ with respect to $x$ is the slope $b$.
The first-order Taylor series of the score function at the true parameter $\theta$ gives the linear approximation (ie. $y = a + bx$) and from it we get (4.25) by plugging in $\theta$ for $x$, $\hat{\theta}$ for $x_0$, $\dot{l}_x(\theta)$ for $y$ and the $\ddot{l}_x(\theta) = \frac{\partial}{\partial\theta}\dot{l}_x(\theta)$ for $b$. The relationship is only approximate.
PS. You seem to use the term "bias" to mean the difference $\hat{\theta}-\theta$ (for a particular dataset $x$); this is the error of the MLE. For any $x$, the error probably won't be exactly 0. But if the expected error is 0, then the MLE is unbiased.

R code to reproduce the figure.
# Number of successes for the negative binomial distribution
negbin.size <- 3

rbern <- function(n, prob) rbinom(n, 1, prob)
rnegbin <- function(n, prob) rnbinom(n, negbin.size, prob)

snorm <- function(x) {
  # Normal model with known variance
  function(mu) {
    n <- length(x)
    n * (mean(x) - mu)
  }
}
spois <- function(x) {
  function(lambda) {
    n <- length(x)
    n * (mean(x) / lambda - 1)
  }
}
sbern <- function(x) {
  function(p) {
    n <- length(x)
    n * (mean(x) / p - (1 - mean(x)) / (1 - p))
  }
}
scauchy <- function(x) {
  # Cauchy model with known scale
  function(x0) {
    score.x0 <- function(x0) {
      xt <- (x - x0)
      sum(2 * xt / (1 + xt^2))
    }
    sapply(x0, score.x0)
  }
}
snegbin <- function(x) {
  # Negative Binomial with known size
  function(p) {
    n <- length(x)
    xbar <- mean(x)
    n * (negbin.size / p - xbar / (1 - p))
  }
}
sexp <- function(x) {
  function(lambda) {
    n <- length(x)
    xbar <- mean(x)
    n * (1 / lambda - xbar)
  }
}

library("glue")
library("patchwork")
library("tidyverse")
iccmlr::theme_set_amazon()

plot_score <- function(theta, n, rdist, sdist, xlims,
                       reps = 20, title = waiver()) {
  score_function <-
    rep(theta, reps) %>%
    map(
      rdist,
      n = n
    ) %>%
    map(function(x) {
      geom_function(
        fun = sdist(x),
        alpha = 0.25
      )
    })

  ggplot() +
    geom_vline(
      xintercept = theta,
      color = "#DF536B"
    ) +
    geom_hline(
      yintercept = 0
    ) +
    score_function +
    labs(x = "θ", y = "score(θ)", title = title) +
    lims(x = xlims) +
    theme(
      plot.title = element_text(size = 10)
    )
}

set.seed(2023)

n <- 10

theta <- 10
theta.range <- c(5, 15)
p <- 0.5
p.range <- c(0.1, 0.9)

# continuous models
p1 <- plot_score(theta, n, rnorm, snorm,
  xlims = theta.range,
  title = glue("normal({theta},σ=1)")
)
p2 <- plot_score(theta, n, rexp, sexp,
  xlims = theta.range,
  title = glue("exponential({theta})")
)
p3 <- plot_score(theta, n, rcauchy, scauchy,
  xlims = theta.range,
  title = glue("cauchy({theta},γ=1)")
)

# discrete models
p4 <- plot_score(p, n, rbern, sbern,
  xlims = p.range,
  title = glue("bernoulli({p})")
)
p5 <- plot_score(p, n, rnegbin, snegbin,
  xlims = p.range,
  title = glue("negative-binomial(r = {negbin.size},{p})")
)
p6 <- plot_score(theta, n, rpois, spois,
  xlims = theta.range,
  title = glue("poisson({theta})")
)

(p1 + p2 + p3) / (p4 + p5 + p6)

