1
$\begingroup$

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?

$\endgroup$
3
  • $\begingroup$ Could you be more specific in your question statement? Exactly what does "what's going on here" refer to? $\endgroup$
    – whuber
    Feb 13 at 18:57
  • $\begingroup$ Done (in the final paragraph) $\endgroup$
    – Eoin
    Feb 13 at 21:27
  • $\begingroup$ Thank you. You can find many answers by searching our site about (a) when and why to use $n$ instead of $n-1$ in the denominator of a variance estimate and, separately, (b) what an unbiased variance estimate is and how Jensen's Inequality implies the associated SD estimate is biased; and (c) why most MLEs are biased. The short answer is that nothing in MLE assures parameter estimates are unbiased. $\endgroup$
    – whuber
    Feb 13 at 21:29

1 Answer 1

3
+200
$\begingroup$

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.

$\endgroup$
3
  • $\begingroup$ This is excellent, thank you. Should $S_{corrected}$ in the last paragraph be uncorrected? $\endgroup$
    – Eoin
    Feb 14 at 9:00
  • $\begingroup$ It's beyond the scope of my original question, but any pointers to further reading on the $n - 1$ thing, or why maximum likelihood is biased in this regard, would be appreciated! Even a good textbook that covers this content would be useful. $\endgroup$
    – Eoin
    Feb 14 at 9:04
  • 1
    $\begingroup$ Re bias in ML Insight is afforded by a basic obvious fact of ML: because it finds its solutions by optimizing an objective function, the solutions survive any one-to-one transformation of the parameters. Jensen's Inequality then implies that if the MLE for a parameter $\theta$ is unbiased and $f$ is a convex transformation, then the MLE of $f(\theta)$ can't be unbiased. One example is the variance $v$ of a Normal distribution: if the MLE of $v$ even were unbiased (it's not), the MLE of $\sigma=\sqrt v$ would have to be biased. $\endgroup$
    – whuber
    Feb 14 at 14:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.