Relationships between sample/population standard deviation, standard error, and maximum likelihood

Question

I'm interested in using Fisher information to calculate standard errors for maximum likelihood estimates, as discussed, e.g. here and here.

To make sure everything was implemented correctly, I was comparing my results for a simple model, $x \sim \text{Normal}(\mu, \sigma)$, to the standard error of the mean calculated in the standard way, $\text{SEM} = \frac{s_\text{Corrected}}{\sqrt{n}}$, where $s_\text{Corrected}$ is the corrected sample standard deviation (that is, dividing by $n - 1$).

library(dplyr)
library(purrr)

# Simulate data
n = 20
ncol = 5
X = matrix(NA, n, ncol)
for(i in 1:ncol){
  X[,i] = rnorm(n, mean = i, sd = i)
}

means = apply(X, 2, mean)
(sample_sds   = apply(X, 2, sd))
## [1] 1.032571 2.100124 2.595767 4.144215 5.024134
(sems_from_sample_sd = sample_sds / sqrt(n))
## [1] 0.2308898 0.4696021 0.5804311 0.9266746 1.1234306

For reasons that will become clear shortly, I also calculated something similar, $\text{SEM}_\text{Alt} = \frac{s_\text{Uncorrected}}{\sqrt{n}}$, where $s_\text{Uncorrected}$ is the uncorrected sample standard deviation, AKA the population standard deviation (dividing by $n$).

#' Function to calculate population (rather than sample) sd
sd_population = function(x){
  n = length(x)
  sqrt((n-1)/n) * sd(x)
}

(population_sds = apply(X, 2, sd_population))
## [1] 1.006425 2.046948 2.530040 4.039281 4.896920
(sems_from_population_sds = population_sds / sqrt(n))
## [1] 0.2250436 0.4577115 0.5657342 0.9032107 1.0949847

Finally, I tried estimating the same thing by maximising the likelihood for parameters $\mu$ and $\sigma$, and calculating the Fisher information, following the code linked above.

I find first that the maximum likelihood estimate for $\sigma$ matches the uncorrected standard deviation estimate from above. I think I understand why this is the case (because estimation of standard deviations is complicated).

Second, and more confusingly, the standard errors estimated in this way correspond to the standard errors calculated using the population standard deviation statistic. I find this very hard to explain.

# Estimate parameters and SEs with maximum-likelihood and Fisher information

loglik_function = function(pars, x){
  sum(dnorm(x, pars[1], pars[2], log = T))
}

fit_model = function(x){
  starting_values = c(1, 1)
  fit = optim(par = starting_values, fn = loglik_function,
              control = list(fnscale = -1), x = x, hessian = T)
  fisher_information = solve(-fit$hessian)
  se = sqrt(diag(fisher_information))
  data.frame(term = c('mu', 'sigma'),
         estimate = fit$par,
             se = se)
}

model_fits = apply(X, 2, fit_model)

sd_estimates = map_dbl(model_fits, function(df){
  filter(df, term == 'sigma')$estimate
})
sd_estimates # Matches population SD
## [1] 1.006358 2.046619 2.529634 4.038780 4.897565

mean_ses = map_dbl(model_fits, function(df){
  filter(df, term == 'mu')$se
})
mean_ses # Matches SEs calculated using population SD
## [1] 0.2250285 0.4576378 0.5656434 0.9030987 1.0951288

Can anyone explain, or point me towards an explanation of, what's going on here. Specifically, why does the standard error of the maximum-likelihood estimate correspond to $\frac{s_\text{Uncorrected}}{\sqrt{n}}$, rather than $\frac{s_\text{Corrected}}{\sqrt{n}}$, (which I've always understood to be) the textbook definition of standard error of the mean?

Related to this, are any reasons we should prefer one definition of the standard error over the other?

Answer 1

The standard error of the mean is the standard deviation of the estimator of the mean, so $\sigma/\sqrt{n}$. We don't know the standard error of the mean, so we need to estimate it. As you say, there are (at least) two plausible estimators

• $s_{\text corrected}/\sqrt{n}$
• $s_{\text uncorrected}/\sqrt{n}$

One advantage of $s_{\text corrected}/\sqrt{n}$ is that its square, $s^2_{\text corrected}/n$ is unbiased ($s_{\text corrected}/\sqrt{n}$ itself is not unbiased). Another is that we usually estimate $\sigma$ with $s_{\text corrected}$, so the relationship between the estimator of $\sigma$ and the estimator of $\sigma/\sqrt{n}$ is simpler.

An advantage of $s_{\text uncorrected}/\sqrt{n}$ is that its square $s^2_{\text uncorrected}/n$ is the observed Fisher information (as you found), which is something that can be computed for a very wide range of models, not just for $N(\mu,\sigma^2)$. In contrast, there isn't a simple analogue of $s_{\text corrected}/\sqrt{n}$ for standard errors of most maximum likelihood estimators -- there aren't unbiased estimators of most things. All that maximum likelihood estimator guarantees (under certain conditions) is that the bias is small compared to the standard error when the sample size is large.

Another advantage of $s_{\text uncorrected}/\sqrt{n}$ is that $s_{\text uncorrected}/\sqrt{n}$ is the estimator of the standard error obtained by plugging the maximum likelihood estimator of $\sigma$ into the formula $\sigma/\sqrt{n}$, an approach that generalises to other estimators than the mean.

So, there isn't a preferred estimator on a principled basis; there are reasons to prefer either option. The corrected estimator is more commonly used and is a simple function of the unbiased estimator of variance; the uncorrected estimator generalises to maximum likelihood estimation in many different models, not just $N(\mu, \sigma^2)$

I think it's probably more common to use $s_{\text corrected}/\sqrt{n}$ under the name "standard error of the mean", and I'd guess that the main reason is that you get it by plugging the usual estimate $s_{\text corrected}$ of $\sigma$ into the standard error formula $\sigma/\sqrt{n}$. Fortunately, it rarely matters much which estimator you use (though you should know which one you're using).

Some people do define the standard error of the mean as $s_{\text corrected}/\sqrt{n}$, but it's not a good definition and they probably shouldn't. I'm glad to see that Wikipedia gets it right.