When does maximum likelihood work and when it doesn't? I'm confused about the maximum likelihood method as compared to e.g. computing the arithmetic mean.
When and why does maximum likelihood produce "better" estimates than e.g. arithmetic mean? How is this verifiable?
 A: Let's interpret "computing the arithmetic mean" as estimation using the Method of Moments (MoM). I believe that's faithful to the original question since the method substitutes sample averages for theoretical ones. It also address @Xi'an's concern about an arbitrary parameter (from an arbitrary model).
If you're still with me, then I think a great place to go is Examples where method of moments can beat maximum likelihood in small samples? The question text points out that "Maximum likelihood estimators (MLE) are asymptotically efficient; we see the practical upshot in that they often do better than method of moments (MoM) estimates (when they differ)," and seeks specific cases where MoM estimators achieve a smaller mean squared error than its MLE counterpart. A few examples that are provided are in the context of linear regression, the two-parameter Inverse Gaussian distribution, and an asymmetric exponential power distribution.
This idea of "asymptotic efficiency" means that maximum likelihood estimators are probably close to using the data to its fullest potential (to estimate the parameter in question), a guarantee you don't get with method of moments in general. While maximum likelihood is not always "better" than working with averages, this efficiency property (if only in the limit) makes it a go-to method for most frequentists. Of course, the contrarian could argue that with the increasing size of data sets, if you're pointing at the right target with a function of averages, go with it.
A: There are several famous examples where maximum likelihood (ML) doesn't provide the best solution. See Lucien Le Cam's 1990 paper: "Maximum Likelihood: an introduction" [1], which is from his invited lectures at the Univ. of Maryland. 
The example that I like the most, because it's so straightforward, is this:
Consider two sequences of independent random variables $X_j$ and $Y_j$ indexed by $j = 1,...,n$. Let's assume that $X_j\sim N(\mu_j,\sigma^2)$ and $Y_j\sim N(\mu_j,\sigma^2)$. In other words, for each $j$ the pair $X_j$ and $Y_j$ are identically distributed with the same mean and variance and the mean is a function of $j$. 
What is the ML estimate of $\sigma^2$?
I won't ruin the fun by giving you the answer, but (no surprise) there are two ways to solve this using ML and they give different solutions. One is the "arithmetic mean" of the squared residuals (as one would expect), and the other is half the arithmetic mean. You can find the answer here on my Github page. 
A: While the arithmetic mean $\bar{x}$ may sound as the "natural" estimator, one could ask why it should be preferred to the MLE! The only sure property associated with the arithmetic mean is that it is an unbiased estimator of $\mathbb{E}[X]$ when this expectation is defined. (Think of the Cauchy distribution as a counter-example.) The later indeed enjoys a wide range of properties under regularity conditions on the likelihood function. To borrow from the wikipedia page, the MLE is


*

*consistent

*asymptotically normal

*efficient in that it achieves the minimum asymptotic variance

*invariant under bijective transforms

*within the parameter set even for constrained parameter sets


In comparison with the arithmetic mean, most of those properties are also satisfied for regular enough distributions. Except 4 and 5. In the case of exponential families, the MLE and the arithmetic mean are identical for estimating the parameter in the mean parameterisation (but not for other parameterisations). And the MLE exists for a sample from the Cauchy distribution.
However, when turning to finite sample optimality properties like minimaxity or admissibility, it may happen that the MLE is neither minimax nor admissible. For instance, the Stein effect shows there exist estimators with a smaller quadratic risk for all values of the parameter under some constraints on the distribution of the sample and the dimension of the parameter. This is the case when $x\sim\mathcal{N}_p(\theta,I_p)$ and $p\ge 3$.
