Methods of Point Estimation: Method of Moments & Maximum Likelihood

The Method of Moments is historically the oldest systematic general technique for parameter estimation. Its underlying philosophy rests on the analogy principle: if a theoretical probability distribution $P_\theta$ depends on $k$ unknown parameters $\boldsymbol{\theta} = (\theta_1, \dots, \theta_k)^T$, and if the first $k$ population moments exist as functions of $\boldsymbol{\theta}$, then equating these theoretical moments to their empirical sample counterparts generates a system of $k$ simultaneous equations whose solution yields point estimates for $\boldsymbol{\theta}$.

Population Moments and Sample Moments

Let $X_1, X_2, \dots, X_n$ be an independent, identically distributed random sample from a distribution with probability density or mass function $f(x; \boldsymbol{\theta})$.

The $r$-th population raw moment is defined as:

assuming the integral converges absolutely ($\mathbb{E}[|X|^r] < \infty$). The $r$-th sample raw moment is the empirical arithmetic average:

Similarly, the $r$-th population central moment is $\mu_r(\boldsymbol{\theta}) = \mathbb{E}[(X - \mu)^r]$, and the $r$-th sample central moment is:

The MOM System of Equations

To estimate a $k$-dimensional parameter vector $\boldsymbol{\theta} = (\theta_1, \dots, \theta_k)^T$, we set up the first $k$ moment-matching equations:

Assuming the mapping $\boldsymbol{\alpha}: \Theta \to \mathbb{R}^k$ given by $\boldsymbol{\alpha}(\boldsymbol{\theta}) = (\alpha_1(\boldsymbol{\theta}), \dots, \alpha_k(\boldsymbol{\theta}))^T$ is one-to-one and continuously invertible on its range, the Method of Moments Estimator (MOME) $\tilde{\boldsymbol{\theta}}$ is the unique solution:

Advantages and Operational Limitations of MOM

The Method of Moments possesses several appealing practical virtues:

1. Computational Simplicity: In many classical distribution families (such as Gamma, Beta, Normal, Poisson, and Weibull), the moment equations reduce to straightforward algebraic formulas or low-degree polynomials that can be solved in closed form without iterative numerical optimization.

2. Distribution-Free Moments: Evaluating population moments often requires only basic calculus, avoiding the need to evaluate normalizing constants or complex likelihood normalizers.

However, the Method of Moments suffers from severe theoretical and practical deficiencies:

• Severe Inefficiency: MOM estimators often fail to attain the Cramér–Rao Lower Bound and discard substantial amounts of Fisher Information. In heavy-tailed distributions or mixtures, the asymptotic variance of MOM estimators can be orders of magnitude larger than that of maximum likelihood estimators.
• Inadmissibility and Boundary Violations: Because sample moments $m_r$ are unconstrained random variables, the resulting estimate $\tilde{\boldsymbol{\theta}} = \boldsymbol{\alpha}^{-1}(\mathbf{m})$ can fall outside the physically valid parameter space $\Theta$. For example, when estimating Beta parameters $(\alpha, \beta)$, a large sample variance $c_2 > m_1(1 - m_1)$ can produce negative parameter estimates $\tilde{\alpha} < 0$, which is physically impossible.
• Sensitivity to Outliers: High-order sample moments $m_3, m_4$ involve cubing and quadrupling observations, making them hypersensitive to extreme values and measurement errors.

Despite their finite-sample inefficiencies, Method of Moments estimators possess excellent asymptotic properties under broad regularity conditions.

Theorem: Strong Consistency of MOM Estimators

Let $X_1, \dots, X_n$ be an i.i.d. sample from $f(x; \boldsymbol{\theta})$ such that the first $k$ population moments $\alpha_1(\boldsymbol{\theta}), \dots, \alpha_k(\boldsymbol{\theta})$ exist and are finite. Let $\mathbf{g} = \boldsymbol{\alpha}^{-1}: \mathbb{R}^k \to \Theta$ be continuous at the true moment vector $\boldsymbol{\alpha}(\boldsymbol{\theta}_0)$. Then the MOM estimator $\tilde{\boldsymbol{\theta}}_n = \mathbf{g}(\mathbf{m}_n)$ is strongly consistent:

Proof:

By Kolmogorov's Strong Law of Large Numbers (SLLN), each sample moment converges almost surely to its corresponding theoretical population moment:

In vector notation, $\mathbf{m}_n \xrightarrow{\text{a.s.}} \boldsymbol{\alpha}(\boldsymbol{\theta}_0)$. Because $\mathbf{g}$ is continuous at $\boldsymbol{\alpha}(\boldsymbol{\theta}_0)$, the Continuous Mapping Theorem immediately yields:

Asymptotic Normality via the Multivariate Delta Method

If the distribution possesses finite moments up to order $2k$ (so that $\mathbb{E}[X^{2k}] < \infty$), the multivariate Central Limit Theorem guarantees that the sample moment vector is asymptotically normal:

where the covariance matrix entries are given by:

If $\mathbf{g}(\mathbf{u})$ is continuously differentiable with Jacobian matrix $\mathbf{J}_{\mathbf{g}}(\mathbf{u}) = \nabla \mathbf{g}(\mathbf{u})$, the Multivariate Delta Method establishes that:

This provides closed-form asymptotic variance expressions for MOM estimators, enabling the construction of large-sample confidence intervals and hypothesis tests.

The foundational method of modern parametric inference is Ronald A. Fisher's Method of Maximum Likelihood. The core intuition of maximum likelihood estimation is that among all conceivable parameter values $\theta \in \Theta$, the preferred point estimate is the parameter value that assigns the highest probability (or probability density) to the realized, observed empirical data.

Definition of Maximum Likelihood Estimator

Let $\mathbf{x} = (x_1, \dots, x_n)$ be the realized sample from a model with joint density $f_n(\mathbf{x}; \boldsymbol{\theta})$. The likelihood function is $L(\boldsymbol{\theta}; \mathbf{x}) = f_n(\mathbf{x}; \boldsymbol{\theta})$. A Maximum Likelihood Estimator (MLE) $\hat{\boldsymbol{\theta}} = \hat{\boldsymbol{\theta}}(\mathbf{X})$ is any value in the parameter space $\Theta$ that globally maximizes the likelihood function:

Because the natural logarithm is a strictly increasing monotonic transformation on $(0, \infty)$, maximizing $L(\boldsymbol{\theta}; \mathbf{x})$ is mathematically identical to maximizing the log-likelihood function:

The Likelihood Equations

When $\Theta$ is an open subset of $\mathbb{R}^k$ and $\ell(\boldsymbol{\theta}; \mathbf{x})$ is continuously differentiable, any interior local maximum must satisfy Fermat's first-order stationary condition. The gradient of the log-likelihood is the Score Vector $\mathbf{S}(\boldsymbol{\theta}; \mathbf{x}) = \nabla_{\boldsymbol{\theta}} \ell(\boldsymbol{\theta}; \mathbf{x})$. Setting the score vector equal to zero generates the Likelihood Equations:

Second-Order Sufficiency Conditions and Hessian Curvature

To verify that a critical point $\hat{\boldsymbol{\theta}}$ satisfying the score equations is indeed a local maximum (rather than a saddle point or minimum), the Hessian matrix of second partial derivatives:

must be strictly negative-definite when evaluated at $\hat{\boldsymbol{\theta}}$:

Equivalently, the Observed Information Matrix $J(\hat{\boldsymbol{\theta}}; \mathbf{x}) = -\mathcal{H}(\hat{\boldsymbol{\theta}}; \mathbf{x})$ must be strictly positive-definite: all its eigenvalues must be strictly positive.

Solving the first-order likelihood equations $\nabla \ell = \mathbf{0}$ does not automatically identify the global MLE. Three distinct structural complications can arise:

1. Strictly Concave Log-Likelihoods (Global Uniqueness):

If the log-likelihood function $\ell(\boldsymbol{\theta}; \mathbf{x})$ is strictly concave over a convex parameter space $\Theta$—meaning its Hessian $\mathcal{H}(\boldsymbol{\theta}; \mathbf{x})$ is strictly negative-definite everywhere on $\Theta$—then the log-likelihood surface is strictly unimodal. Any stationary point satisfying $\nabla \ell(\hat{\boldsymbol{\theta}}) = \mathbf{0}$ is guaranteed to be the unique global maximum. Example: In natural exponential families with open convex parameter spaces (such as Gaussian, Poisson, Binomial, and Gamma), the log-likelihood is strictly concave, ensuring that the MLE is unique whenever a solution exists.

2. Multimodal Likelihoods (Local Maxima and Spurious Roots):

In complex models—including finite mixture distributions, neural networks, Gaussian processes, and Cauchy distributions—the likelihood surface features multiple local peaks, ridges, and saddle points. Standard root-finding algorithms (such as Newton–Raphson) can converge to an inferior local maximum or saddle point. Finding the true global MLE requires multiple random restarts, global annealing, or grid searches.

3. Boundary Solutions:

If the maximum of $L(\boldsymbol{\theta})$ occurs on the topological boundary $\partial \Theta$ of the parameter space, the score vector at the maximum may not be zero ($\nabla \ell(\hat{\boldsymbol{\theta}}) \ne \mathbf{0}$). In constrained optimization (such as enforcing non-negative variances $\sigma^2 \ge 0$ in random-effects models), the MLE must be identified using Karush–Kuhn–Tucker (KKT) conditions.

When the support of the probability density depends on the unknown parameter $\theta$, standard calculus techniques collapse. The likelihood function is not differentiable across the support boundary, and the MLE is typically determined by order statistics rather than score equations.

Canonical Example: Uniform Distribution $\text{Uniform}(0, \theta)$

Let $X_1, \dots, X_n \sim \text{Uniform}(0, \theta)$ where $\theta > 0$. The joint likelihood function is:

where $x_{(n)} = \max(x_1, \dots, x_n)$. Notice the analytical behavior of $L(\theta)$ as a function of $\theta$:

• If $\theta < x_{(n)}$, at least one observation falls outside the interval $[0, \theta]$, so $L(\theta) = 0$.
• If $\theta \ge x_{(n)}$, all observations are inside $[0, \theta]$, so $L(\theta) = \frac{1}{\theta^n}$.

On the domain $\theta \in [x_{(n)}, \infty)$, the function $L(\theta) = \theta^{-n}$ is a strictly decreasing function of $\theta$. Its derivative $\frac{dL}{d\theta} = -n \theta^{-(n+1)} < 0$ is strictly negative everywhere, so setting the derivative to zero yields no solution! To maximize a strictly decreasing function over the interval $[x_{(n)}, \infty)$, we must choose the smallest permissible value in the domain. Therefore, the Maximum Likelihood Estimator is:

Notice that $\hat{\theta}$ is an order statistic located precisely on the boundary of the feasible domain. Furthermore, because $X_i < \theta$ with probability $1$, we have $X_{(n)} < \theta$ with probability $1$. Thus, the MLE systematically underestimates $\theta$:

The MLE is biased, but asymptotically unbiased as $n \to \infty$.

Shifted Exponential and Pareto Models

Similarly, for a shifted exponential distribution $f(x; \theta, \lambda) = \lambda e^{-\lambda (x - \theta)} \mathbf{1}_{\{x \ge \theta\}}$, the likelihood is strictly increasing with respect to $\theta$ up to the threshold $x_{(1)} = \min(x_1, \dots, x_n)$, making the MLE:

In non-regular models, maximum likelihood estimators frequently achieve super-efficient convergence rates of order $\mathcal{O}(n)$ rather than the standard $\mathcal{O}(\sqrt{n})$ Gaussian rate.

One of the most powerful and celebrated practical advantages of maximum likelihood estimation over all other estimation frameworks (including unbiasedness and Method of Moments) is its functional invariance property, originally proved in full generality by Zeena (1966).

Theorem: The Invariance Property of MLEs

Let $\hat{\boldsymbol{\theta}}$ be the Maximum Likelihood Estimator of a parameter vector $\boldsymbol{\theta} \in \Theta$. Let $g: \Theta \to \Omega$ be any arbitrary function (which need not be one-to-one or invertible). Then the Maximum Likelihood Estimator of the transformed parameter $\tau = g(\boldsymbol{\theta})$ is:

Proof:

If $g$ is a one-to-one (bijective) mapping, the result is trivial: reparameterizing the likelihood $L^*(\tau) = L(g^{-1}(\tau))$ preserves the supremum, so $\hat{\tau} = g(\hat{\boldsymbol{\theta}})$. Now consider the general case where $g$ is not one-to-one (many-to-one). To define the likelihood function for $\tau$, we utilize the concept of the profile likelihood (induced likelihood). For any $\tau \in \Omega$, define:

The maximum likelihood estimator of $\tau$ is defined as the value $\hat{\tau}$ that maximizes $L_g(\tau; \mathbf{x})$ over $\Omega$:

Now, evaluate the induced likelihood at the point $\tau^* = g(\hat{\boldsymbol{\theta}})$:

Because $\hat{\boldsymbol{\theta}}$ itself belongs to the set $\{\boldsymbol{\theta} \in \Theta : g(\boldsymbol{\theta}) = g(\hat{\boldsymbol{\theta}})\}$, the supremum over that set is at least as large as $L(\hat{\boldsymbol{\theta}}; \mathbf{x})$:

Combining the two inequalities:

Therefore, $\hat{\tau} = g(\hat{\boldsymbol{\theta}})$ achieves the global maximum of the induced likelihood function. This completes the proof.

Significance of Invariance:

Invariance gives maximum likelihood estimators immense practical flexibility:

• If $\hat{\sigma}^2$ is the MLE of variance, then the MLE of standard deviation is immediately $\hat{\sigma} = \sqrt{\hat{\sigma}^2}$.
• If $\hat{p}$ is the MLE of a binomial success probability, the MLE of the odds ratio $\frac{p}{1-p}$ is simply $\frac{\hat{p}}{1 - \hat{p}}$.
• If $\hat{\lambda}$ is the MLE of a Poisson rate, the MLE of the zero-event probability $e^{-\lambda}$ is $e^{-\hat{\lambda}}$.

Recall from Unit 1 that unbiased estimators never possess this property: $g(T_{\text{UMVUE}})$ is almost never unbiased for $g(\theta)$. Maximum likelihood estimators liberate the practitioner from re-deriving estimators under every non-linear reparameterization.

In multiparameter models, the parameter vector often partitions into structural parameters of interest $\boldsymbol{\psi}$ and nuisance parameters $\boldsymbol{\lambda}$:

We wish to perform inference on $\boldsymbol{\psi}$ without being impeded by uncertainty in $\boldsymbol{\lambda}$.

Definition of the Profile Likelihood

For each fixed candidate value of $\boldsymbol{\psi}$, define the constrained maximum likelihood estimate of the nuisance parameter:

The Profile Log-Likelihood Function for $\boldsymbol{\psi}$ is defined as:

The overall global MLE of the structural parameter $\hat{\boldsymbol{\psi}}$ maximizes the profile likelihood:

Curvature of the Profile Likelihood and the Schur Complement

A celebrated theorem in likelihood inference connects the curvature of the profile likelihood directly to the full Fisher Information matrix. Let $\mathcal{H}(\boldsymbol{\theta})$ be the full Hessian matrix partitioned as:

By the multivariable chain rule, the second derivative of the profile log-likelihood evaluated at the MLE is:

This is the Schur complement of $\mathcal{H}_{\lambda\lambda}$ in the full Hessian matrix! Inverting this expression reveals that:

The curvature of the profile log-likelihood automatically incorporates the exact variance inflation caused by estimating the unknown nuisance parameters! This enables the construction of valid profile likelihood confidence intervals:

In most non-trivial statistical models, the likelihood equations $\nabla \ell(\boldsymbol{\theta}) = \mathbf{0}$ cannot be solved analytically. Numerical optimization algorithms are required to compute the MLE.

1. The Newton–Raphson Algorithm

The Newton–Raphson method approximates the log-likelihood locally around the current iterate $\boldsymbol{\theta}^{(t)}$ via a second-order Taylor series expansion:

Differentiating with respect to $\boldsymbol{\theta}$ and setting the gradient to zero yields the celebrated iterative update rule:

where $J(\boldsymbol{\theta}^{(t)}) = -\mathcal{H}(\boldsymbol{\theta}^{(t)})$ is the observed information matrix. Near the global maximum, Newton–Raphson exhibits quadratic convergence (doubling the number of correct digits at each iteration). However, if the starting point $\boldsymbol{\theta}^{(0)}$ is far from the peak, the observed Hessian may not be negative-definite, causing the algorithm to diverge.

2. Fisher's Method of Scoring

Ronald Fisher proposed replacing the observed information matrix $J(\boldsymbol{\theta}^{(t)})$ with the Expected Fisher Information Matrix $\mathbf{I}(\boldsymbol{\theta}^{(t)}) = \mathbb{E}[J(\boldsymbol{\theta}^{(t)})]$:

Fisher scoring has significant advantages:

• Because $\mathbf{I}(\boldsymbol{\theta})$ is the covariance matrix of the score vector, it is guaranteed to be positive semi-definite everywhere in $\Theta$, ensuring that each step moves in an ascent direction.
• In generalized linear models (GLMs) with canonical link functions, the observed and expected information matrices coincide ($J \equiv \mathbf{I}$), making Fisher scoring and Newton–Raphson identical (known as Iteratively Reweighted Least Squares / IRLS).

Quasi-Newton Optimization: The BFGS Algorithm in Likelihood Maximization

In high-dimensional parametric modeling (such as multinomial logistic regression, neural networks, or spatial autoregressions with $p > 100$ parameters), calculating the full $p \times p$ analytical Hessian matrix $\mathcal{H}(\boldsymbol{\theta})$ at each iteration is computationally prohibitive or analytically intractable. The Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm is the most popular quasi-Newton optimization method. Instead of computing the true Hessian, BFGS builds a sequence of symmetric positive-definite approximations $\mathbf{B}_t \approx \mathcal{H}(\boldsymbol{\theta}^{(t)})^{-1}$ using only gradient (score) evaluations:

where:

• $\mathbf{s}_t = \boldsymbol{\theta}^{(t+1)} - \boldsymbol{\theta}^{(t)}$ is the parameter step vector.
• $\mathbf{y}_t = \nabla \ell(\boldsymbol{\theta}^{(t+1)}) - \nabla \ell(\boldsymbol{\theta}^{(t)})$ is the change in the score gradient.
• $\rho_t = \frac{1}{\mathbf{y}_t^T \mathbf{s}_t}$ is a normalization scalar.

The BFGS algorithm guarantees that the search direction $\mathbf{d}_t = \mathbf{B}_t \nabla \ell(\boldsymbol{\theta}^{(t)})$ is always an ascent direction, avoids matrix inversions entirely via rank-2 updates, and achieves superlinear convergence to the Maximum Likelihood Estimator.

When data contain missing values, latent cluster assignments, or censored observations, the marginal likelihood function is an intractable sum or integral over all possible hidden configurations. In 1977, A. P. Dempster, N. M. Laird, and D. B. Rubin formalized the Expectation–Maximization (EM) algorithm as a general iterative framework for maximum likelihood estimation in incomplete-data problems.

Complete Data vs Observed Data

Let $\mathbf{Y}$ denote the observed data, and let $\mathbf{Z}$ denote the unobserved latent variables or missing data. The combination $\mathbf{X} = (\mathbf{Y}, \mathbf{Z})$ is the complete data. The observed data likelihood is the marginal integral:

Evaluating this integral directly is often computationally impossible. The EM algorithm avoids integration by operating iteratively on the complete-data log-likelihood $\ln f(\mathbf{Y}, \mathbf{Z}; \boldsymbol{\theta})$ through two alternating steps:

The Two Steps of the EM Algorithm:

1. Expectation Step (E-step):

Given the current parameter iterate $\boldsymbol{\theta}^{(t)}$, compute the conditional expectation of the complete-data log-likelihood with respect to the latent variables $\mathbf{Z}$, conditioned on the observed data $\mathbf{Y}$ and current parameters $\boldsymbol{\theta}^{(t)}$:

2. Maximization Step (M-step):

Update the parameter estimate by globally maximizing the surrogate function $Q(\boldsymbol{\theta} \mid \boldsymbol{\theta}^{(t)})$ over $\Theta$:

Monotonicity and Convergence (Dempster–Laird–Rubin Theorem)

By Jensen's Inequality, every EM iteration is guaranteed to never decrease the observed log-likelihood:

The observed likelihood sequence $\{\ell(\boldsymbol{\theta}^{(t)})\}$ converges monotonically to a local maximum or stationary point of the likelihood surface, providing exceptional numerical stability.

While Maximum Likelihood Estimators are asymptotically efficient under a perfectly specified model, they can be catastrophically sensitive to minor model misspecifications and heavy-tailed contamination (outliers).

To overcome this brittleness, Peter J. Huber (1964) developed the theory of M-Estimators (generalized maximum likelihood estimators). Instead of minimizing $-\sum \ln f(x_i; \theta)$, an M-estimator minimizes a general objective function $\rho$:

Differentiating with respect to $\theta$ yields the estimating equation:

where $\psi(u) = \rho'(u)$ is the score function.

The Influence Function and Huber's Robust Criterion

Frank Hampel introduced the Influence Function $\text{IF}(x; T, F)$, which measures the asymptotic bias induced in an estimator $T$ by an infinitesimal point contamination at location $x$:

For an M-estimator, the influence function is directly proportional to $\psi(x)$:

If $\psi(x)$ is unbounded (as in classical sample mean where $\rho(u) = u^2 / 2$ and $\psi(u) = u$), a single catastrophic outlier at $x = \infty$ can drag the estimate to infinity! The breakdown point of the sample mean is $0\%$. To bound the influence function and protect against outliers, Huber introduced the Huber Loss Function:

The Huber $\psi$-function acts like the mean for small deviations ($|u| \le k$) and like the median for large outliers ($|u| > k$), producing estimators that achieve $95\%$ Gaussian efficiency while remaining completely bounded and robust against severe data contamination.

Asymptotic Efficiency and Bahadur Efficiency of Maximum Likelihood Estimators

Beyond classical first-order efficiency defined by the Cramér–Rao lower bound and asymptotic variance, modern mathematical statistics evaluates estimators through higher-order asymptotic frameworks and large deviation theory. In 1960, R. R. Bahadur introduced an exact measure of asymptotic efficiency based on the exponential rate of decay of tail probabilities.

Let $T_n$ be a consistent estimator of a scalar parameter $\theta$. The Bahadur efficiency assesses how rapidly the probability of substantial estimation error $\mathbb{P}_\theta(|T_n - \theta| > \epsilon)$ tends to zero as the sample size $n \to \infty$. By the Sanov and Cramér large deviation theorems, under regular exponential families, this probability decays exponentially:

where $c(T, \epsilon)$ represents the Bahadur slope of the estimator. Bahadur established that for any consistent estimator, the slope is bounded from above by the Kullback–Leibler divergence between the true parameter and the perturbed parameter:

Remarkably, under standard regularity conditions, the maximum likelihood estimator achieves this theoretical upper bound exactly:

Consequently, the maximum likelihood estimator is optimal not only in terms of its localized Gaussian asymptotic variance around the true parameter, but also globally across the entire sample space in penalizing large, catastrophic deviations. This profound result provides a deep mathematical justification for selecting maximum likelihood estimators in safety-critical scientific and engineering applications.