Maximum Likelihood Estimation

Suppose we watch twenty visitors arrive at a website and twelve of them click a button. We want to estimate the underlying click probability. Or suppose a sensor returns five readings around 10 volts, and we want to estimate the true voltage and the measurement noise. In both cases we have data and a probability model with unknown parameters. Maximum likelihood estimation gives us a single, general recipe for turning those data into estimates. You choose the parameter values that make the observed data as plausible as possible under the model.

The likelihood is not the probability that a parameter is correct. Instead, we hold the observed data fixed and vary the parameter. For each candidate parameter value, we ask how likely the data would have been if that value had generated them. The maximum likelihood estimate, usually written $\hat\theta$, is the parameter value with the largest likelihood. This reverses our usual perspective. A probability distribution varies the data while its parameter is fixed; a likelihood varies the parameter while the observed data are fixed. Keeping that distinction clear prevents many later mistakes, especially when we compare likelihoods with posterior distributions or confidence intervals.

Taking a logarithm turns a product of many probabilities into a sum of many terms, which is easier to differentiate and far less prone to numerical underflow. The resulting log-likelihood is the central object in this chapter. It connects naturally to loss functions: under independent Gaussian noise with a common positive variance, maximizing the likelihood over the mean parameters is equivalent to minimizing the sum of squared errors we studied in Sum of Squared Errors (SSE). You will derive the main estimators by hand, inspect the curvature that produces standard errors, and then implement the same calculations with NumPy and SciPy. The ideas here also prepare us for the next chapter on gradient descent and numerical optimization, where we study iterative parameter fitting and numerical checks.

Maximum likelihood is a general construction rather than a single formula. Once we can write down a probability model for the data, we can form a likelihood, take its logarithm, and optimize. The same recipe produces the sample proportion for Bernoulli data, the sample mean and variance for Gaussian data, logistic regression coefficients for binary outcomes, and many survival and time-series models. This consistency is valuable because you do not need a new estimation principle for every model family.

The method also gives us more than a point estimate. The first derivative of the log-likelihood, called the score, tells us where the likelihood is flat. The second derivative, or curvature, tells us how sharply the likelihood peaks. That curvature produces approximate standard errors and confidence intervals. In regular models, the maximum likelihood estimator is asymptotically normal and asymptotically efficient. Under the regular asymptotic theory, it attains the inverse-information variance bound within the class of regular estimators. You should read that as a large-sample statement, not a promise that every finite-sample estimate is unbiased or precise.

Finally, maximum likelihood connects statistical modeling to machine learning losses. Gaussian likelihood leads to squared error, Bernoulli likelihood leads to log loss, and Poisson likelihood leads to a deviance loss. This makes likelihood a bridge between probability models, optimization objectives, and predictive evaluation. When we see a loss function in a machine learning model, it often has a likelihood interpretation behind it.

Maximum likelihood inherits every assumption of the probability model. If we choose a Gaussian model for data with heavy outliers, the resulting mean can be pulled strongly toward those outliers. If correlated Bernoulli observations are treated as independent, the standard errors can be wrong. Positive covariance terms increase the variance of the sample proportion; negative terms can decrease it. With a misspecified model, the estimator may converge to a pseudo-true parameter that is closest to the data-generating process in a Kullback-Leibler sense. That parameter may not answer the scientific or business question you care about.

Maximum likelihood is also not guaranteed to be unbiased in finite samples. The Gaussian variance estimator is a standard example: it divides the sum of squared errors by $n$, not $n-1$, so it is biased downward on average. At boundaries, such as a Bernoulli probability estimated from all successes, the maximum may sit at the edge of the parameter space, and the usual normal approximation can break down. Some models are not identifiable, meaning different parameter values produce exactly the same distribution for the data, so no amount of data can separate them.

Numerical difficulties are common as well. The likelihood can be flat, multimodal, or unbounded in unconstrained parameterizations. Optimization may stop at a local maximum, fail because of scaling, or wander toward infinity when a parameter approaches a boundary. These are not reasons to avoid likelihood, but they are reasons to check convergence, inspect curvature, and report uncertainty rather than treating a single number as certain.

We can build the likelihood in three steps. First, hold the observed data fixed; they are no longer random quantities to be predicted. Second, choose a probability model that assigns a probability or density to each observation for any candidate parameter value. Third, compare candidate parameter values by asking which one gives the entire observed sample the highest combined plausibility. The formula below is simply the bookkeeping for those three steps.

Suppose we observe independent and identically distributed data $y_1,\dots,y_n$. Let $f(y_i;\theta)$ be the probability mass function for discrete data or the probability density function for continuous data, where $\theta$ is a parameter or a vector of parameters. The likelihood is

$$L(\theta; y_1,\dots,y_n) = \prod_{i=1}^{n} f(y_i;\theta).$$

where:

• $L(\theta; y_1,\dots,y_n)$: the likelihood of the parameter $\theta$ given the full observed sample
• $n$: the number of observations in the sample
• $y_i$: the $i$-th observed data value, held fixed as $\theta$ varies
• $f(y_i;\theta)$: the probability mass (discrete case) or probability density (continuous case) assigned to $y_i$ by the model with parameter $\theta$
• $\prod_{i=1}^{n}$: the product over all observations from $i=1$ to $i=n$

The product symbol $\prod_{i=1}^{n}$ means that we multiply the terms for $i=1$ through $i=n$. That multiplication is not an arbitrary convention. Independence is exactly what makes the joint probability or joint density of the sample factor into one term per observation. If the observations were dependent, the product of marginal terms would not be the joint distribution, and the likelihood would have to be built from the joint model instead. The semicolon in $f(y_i;\theta)$ separates the observed data $y_i$ from the parameter $\theta$; it reminds us that the data are fixed and $\theta$ is the input we are varying. For discrete data, each $f(y_i;\theta)$ is a probability, so the product is the probability of the observed sequence under independence. For continuous data, each term is a density, so the product is a likelihood density rather than the probability of observing an exact value.

The log-likelihood is

$$\begin{aligned} \ell(\theta) &= \log L(\theta; y_1,\dots,y_n) \\ &= \log \prod_{i=1}^{n} f(y_i;\theta) && \text{(substitute the likelihood)} \\ &= \sum_{i=1}^{n} \log f(y_i;\theta) && \text{(log of a product is a sum)}. \end{aligned}$$

The symbols in this expression are:

• $\ell(\theta)$: the log-likelihood, the natural logarithm of the likelihood
• $\log$: the natural logarithm
• $\sum_{i=1}^{n}$: the sum over all observations from $i=1$ to $i=n$
• the remaining symbols are defined as above.

The logarithm is monotone increasing. If $L(\theta_1) > L(\theta_2)$, then $\ell(\theta_1) > \ell(\theta_2)$. Therefore the parameter that maximizes $L$ also maximizes $\ell$. Logs turn products into sums, which are easier to differentiate and much less likely to underflow when $n$ is large. It is common to see negative log-likelihood minimized instead, because optimization software usually minimizes rather than maximizes. These three facts, factorization under independence, monotone transformation, and numerical stability, explain nearly every manipulation that follows.

The likelihood is a function of the parameter, not a probability distribution over the parameter. This point deserves emphasis. If we compute $L(0.3)=0.002$ and $L(0.6)=0.010$, we can say that the observed data are five times more likely under $0.6$ than under $0.3$ in the likelihood ratio sense. We cannot say that there is a 1% probability that the parameter equals $0.6$. To obtain a probability distribution over $\theta$, we would need a prior distribution and a normalization step, which leads to Bayesian inference. Here we remain in the frequentist likelihood framework.

The independence assumption is what allows us to multiply the individual terms. If observations are correlated, the product of marginal probabilities is not the joint probability, and the likelihood must be built from the joint distribution instead. Identically distributed means that every observation has the same distribution. Independent observations may instead have different densities, all indexed by a shared parameter; their likelihood is still the product of those individual densities. Dependence requires a joint model.

For discrete data, the likelihood is a probability of the observed data. For continuous data, it is a density. A density can be greater than one when the distribution is concentrated over a narrow range, so a likelihood value above one is not an error. For estimation, what matters is the relative size of the likelihood at different parameter values, not whether the number looks like a probability.

To turn this function into an estimate, we need a way to locate its peak. The score is the first derivative of the log-likelihood:

$$s(\theta) = \frac{\partial \ell(\theta)}{\partial \theta}.$$

where:

• $s(\theta)$: the score, the first derivative of the log-likelihood with respect to the parameter
• $\theta$: the parameter (or parameter vector) being estimated
• $\partial \ell(\theta)/\partial \theta$: the partial derivative of the log-likelihood with respect to $\theta$

At a differentiable maximum in the interior of the parameter space, the score is zero:

$$s(\hat\theta)=0.$$

This equation is called the score equation or the likelihood equation. In simple models we can solve it algebraically. In more complex models we solve it numerically. The score tells us the local slope of the log-likelihood. A positive score means increasing $\theta$ increases the log-likelihood, and a negative score means decreasing $\theta$ increases it. A zero score means we are at a flat point, which could be a maximum, a minimum, or a saddle point. A negative second derivative gives a strict local maximum in one dimension; a positive one gives a strict local minimum. A zero second derivative is inconclusive. For a vector parameter, negative definiteness of the Hessian is a sufficient local maximum test.

The general formula becomes concrete only after we choose a probability model. You now derive two estimators that appear constantly in practice. The first is the Bernoulli probability, which gives us the sample proportion. The second is the Gaussian mean and variance, which gives us the sample mean and the maximum likelihood variance. These derivations show how the general likelihood formula becomes a concrete estimate.

For Bernoulli data, each observation is either $1$ or $0$. Assume $n>0$, and let $p\in[0,1]$ be the probability of a success. For $0<p<1$, the probability mass is

$$f(y_i;p)=p^{y_i}(1-p)^{1-y_i}.$$

where:

• $y_i$: the $i$-th observation, taking the value $1$ for a success or $0$ for a failure
• $p$: the probability of a success, the unknown parameter
• $1-p$: the probability of a failure

When $y_i=1$, this becomes $p^1(1-p)^0=p$. When $y_i=0$, it becomes $p^0(1-p)^1=1-p$. If we observe $k$ successes and $n-k$ failures, then the product over all observations is

$$L(p)=p^{k}(1-p)^{n-k}.$$

where:

• $k$: the number of observed successes
• $n-k$: the number of observed failures
• $p$: the success probability

The log-likelihood is

$$\ell(p)=k\log p+(n-k)\log(1-p).$$

where:

• $\ell(p)$: the log-likelihood for the Bernoulli success probability $p$
• $k$: the number of observed successes
• $n-k$: the number of observed failures
• $p$: the success probability

Differentiating with respect to $p$ gives

$$\frac{d\ell}{dp}=\frac{k}{p}-\frac{n-k}{1-p}.$$

where:

• $\frac{d\ell}{dp}$: the derivative of the Bernoulli log-likelihood with respect to $p$
• $\frac{k}{p}$: the term contributed by the successes
• $-\frac{n-k}{1-p}$: the derivative term contributed by the failures

Setting this equal to zero, we multiply both sides by $p(1-p)$:

$$k(1-p)-p(n-k)=0.$$

where:

• $k$: the number of observed successes
• $1-p$: the candidate failure probability
• $n-k$: the number of observed failures
• $p$: the success probability

Expanding gives $k-kp-pn+pk=0$, and the $-kp+pk$ terms cancel, leaving $k-pn=0$. Therefore

$$\hat p=\frac{k}{n}.$$

The second derivative is

$$\frac{d^2\ell}{dp^2}=-\frac{k}{p^2}-\frac{n-k}{(1-p)^2}<0,$$

where:

• $\frac{d^2\ell}{dp^2}$: the second derivative of the Bernoulli log-likelihood
• $\frac{k}{p^2}$: the magnitude of the negative curvature contribution from the successes
• $\frac{n-k}{(1-p)^2}$: the magnitude of the negative curvature contribution from the failures

This derivative is negative for $0<p<1$. With $0<k<n$, the stationary point is the unique global maximum. On the closed parameter space $[0,1]$, $k=0$ or $k=n$ puts the maximum at $p=0$ or $p=1$, using the limiting convention $0\log 0=0$. The plug-in Wald standard error is zero there, but that does not establish zero uncertainty. The interior information formula below is not valid at those endpoints.

If we only record the count $k$ rather than the exact sequence of successes and failures, the binomial likelihood is $\binom{n}{k}p^k(1-p)^{n-k}$. The binomial coefficient does not depend on $p$, so it does not change the maximizing value. That is why we can work with the simpler Bernoulli form.

For Gaussian data, suppose $y_i\sim N(\mu,\sigma^2)$. The density is

$$f(y_i;\mu,\sigma^2)=\frac{1}{\sqrt{2\pi\sigma^2}} \exp\left[-\frac{(y_i-\mu)^2}{2\sigma^2}\right].$$

where:

• $\mu$: the mean of the normal distribution
• $\sigma^2$: the variance of the normal distribution
• $\sigma$: the standard deviation, the positive square root of the variance
• $\exp[\cdot]$: the exponential function

The log-likelihood for all $n$ observations is

$$\ell(\mu,\sigma^2) = -\frac{n}{2}\log(2\pi) -\frac{n}{2}\log(\sigma^2) -\frac{1}{2\sigma^2}\sum_{i=1}^{n}(y_i-\mu)^2.$$

where:

• $\ell(\mu,\sigma^2)$: the Gaussian log-likelihood for all $n$ observations
• $n$: the number of observations
• $\log(2\pi)$: the constant log term from the normalizing factor
• $\log(\sigma^2)$: the log of the variance
• $\frac{1}{2\sigma^2}\sum_{i=1}^{n}(y_i-\mu)^2$: the sum of squared deviations divided by twice the variance

To estimate $\mu$, we differentiate with respect to $\mu$. The first two terms do not depend on $\mu$, and the derivative of the squared term brings down a factor of $-2(y_i-\mu)$:

$$\frac{\partial \ell}{\partial \mu} = \frac{1}{\sigma^2}\sum_{i=1}^{n}(y_i-\mu).$$

where:

• $\frac{\partial \ell}{\partial \mu}$: the score with respect to the Gaussian mean
• $\sigma^2$: the variance of the normal distribution
• $y_i-\mu$: the deviation of the $i$-th observation from the candidate mean

Setting this to zero gives

$$\sum_{i=1}^{n}(y_i-\mu)=0 \quad\Longrightarrow\quad \sum_{i=1}^{n}y_i-n\mu=0 \quad\Longrightarrow\quad \hat\mu=\frac{1}{n}\sum_{i=1}^{n}y_i=\bar y.$$

where:

• $\sum_{i=1}^{n}(y_i-\mu)=0$: the score set to zero
• $\sum_{i=1}^{n}y_i$: the sum of all observations
• $n\mu$: the sum if every observation were equal to $\mu$
• $\hat\mu$: the maximum likelihood estimate of the mean
• $\bar y$: the sample mean

The maximum likelihood estimate of the mean is the sample mean. This result does not depend on the value of $\sigma^2$.

To estimate the variance, let $v=\sigma^2$ and write the sum of squared deviations as

$$\mathrm{SSE}=\sum_{i=1}^{n}(y_i-\bar y)^2.$$

After substituting $\mu=\bar y$, the profile log-likelihood for $v>0$ becomes

$$\ell(\bar y,v) = -\frac{n}{2}\log(2\pi) -\frac{n}{2}\log v -\frac{\mathrm{SSE}}{2v}.$$

where:

• $\ell(\bar y,v)$: the Gaussian log-likelihood with the mean fixed at its fitted value and candidate variance $v$
• $n$: the number of observations
• $v$: the candidate variance $\sigma^2$
• $\mathrm{SSE}$: the sum of squared deviations from the sample mean

Differentiating with respect to $v$ gives

$$\frac{\partial \ell}{\partial v} = -\frac{n}{2v}+\frac{\mathrm{SSE}}{2v^2}.$$

Setting this to zero, we multiply by $2v^2$:

$$-nv+\mathrm{SSE}=0 \quad\Longrightarrow\quad \hat v=\frac{\mathrm{SSE}}{n}.$$

where:

• $n$: the number of observations
• $v$: the candidate variance
• $\mathrm{SSE}$: the sum of squared deviations from the sample mean
• $\hat v$: the maximum likelihood variance estimate

If $\mathrm{SSE}>0$, the profile likelihood rises for $v<\mathrm{SSE}/n$ and falls for $v>\mathrm{SSE}/n$, so this is the global maximum over $v>0$. If all observations are identical, $\mathrm{SSE}=0$ and the Gaussian likelihood is unbounded as $v$ tends to zero: no positive-variance MLE is attained. For nondegenerate i.i.d. Gaussian samples, the variance estimator is biased downward because

$$\mathbb{E}\left[\frac{\mathrm{SSE}}{n}\right] =\frac{n-1}{n}\sigma^2.$$

where:

• $\mathbb{E}[\cdot]$: expectation over repeated samples
• $\mathrm{SSE}$: the sum of squared deviations from the sample mean
• $n$: the number of observations
• $\sigma^2$: the true variance

The bias shrinks as $n$ grows, and the estimator is consistent, but it is not unbiased. This is a concrete reminder that maximum likelihood and unbiasedness are different properties.

The Gaussian likelihood also explains why squared error appears so often in regression. If we model $y=X\beta+\varepsilon$ with $\varepsilon\sim N(0,\sigma^2 I)$, the log-likelihood is

$$\ell(\beta,\sigma^2) = -\frac{n}{2}\log(2\pi\sigma^2) -\frac{1}{2\sigma^2}(y-X\beta)^\top(y-X\beta).$$

Here $y$ is a column vector of responses, $X$ is a matrix with one row per observation and one column per predictor, $\beta$ is a column vector of coefficients, and $I$ is the identity matrix. The expression $(y-X\beta)^\top(y-X\beta)$ is exactly the sum of squared errors written in matrix notation. For a fixed $\sigma^2$, maximizing the likelihood over $\beta$ is the same as minimizing the sum of squared errors. We will see the full regression derivation in the later chapters on ordinary least squares, but the likelihood connection is already visible here.

The mean controls the center of the distribution, and the standard deviation controls its width. When you evaluate a likelihood, you are not plotting a distribution over $\mu$ or $\sigma$. You are comparing candidate distributions over the observed data $y$.

The score equation $s(\hat\theta)=0$ is the first-order condition for a maximum. To confirm a maximum and to measure uncertainty, we use the second derivative. Intuitively, the first derivative tells us where the likelihood is flat, while the second derivative tells us how quickly it curves away from that flat point. A sharp peak is informative. A broad, nearly flat peak is not. The observed information is

$$J(\hat\theta)=-\left.\frac{d^2\ell(\theta)}{d\theta^2}\right|_{\theta=\hat\theta}.$$

where:

• $J(\hat\theta)$: the observed information, evaluated at the estimate $\hat\theta$
• $\hat\theta$: the maximum likelihood estimate of the parameter
• $d^2\ell(\theta)/d\theta^2$: the second derivative of the log-likelihood with respect to $\theta$
• the vertical bar with subscript $\theta=\hat\theta$: evaluation of the second derivative at the estimate

For a vector parameter, the second derivative becomes a Hessian matrix, and the observed information is the negative Hessian:

$$J(\hat\theta)=-\nabla^2\ell(\hat\theta).$$

The Fisher information is the expected value of the observed information under the model:

$$I_n(\theta) = \mathbb{E}_{\theta}\left[J(\theta)\right] = \mathbb{E}_{\theta}\left[-\nabla^2\ell(\theta)\right],$$

where:

• $I_n(\theta)$: the Fisher information in the full sample for the parameter $\theta$
• $\mathbb{E}_{\theta}[\cdot]$: the expectation taken over repeated samples generated with the true parameter equal to $\theta$
• $\nabla^2\ell(\theta)$: the Hessian of the log-likelihood (the matrix of second partial derivatives)

For a correctly specified regular i.i.d. model, $I_n(\theta)=nI_1(\theta)$, where $I_1$ is per-observation information. Under the usual consistency and derivative-convergence conditions, $J(\hat\theta)/n$ and $I_n(\hat\theta)/n$ converge to the same nonsingular limit. This statement concerns information per observation; their unnormalized difference need not tend to zero.

At a regular interior estimate with positive information, a large-sample standard error for a scalar parameter is

$$\widehat{\mathrm{se}}(\hat\theta) \approx \sqrt{J(\hat\theta)^{-1}} = \sqrt{\frac{1}{-\ell''(\hat\theta)}}.$$

where:

• $\widehat{\mathrm{se}}(\hat\theta)$: the approximate standard error of the estimate
• $J(\hat\theta)^{-1}$: the reciprocal of the observed information (the inverse information matrix in the vector case)
• $\ell''(\hat\theta)$: the second derivative of the log-likelihood evaluated at the estimate

For a regular interior vector estimate with nonsingular positive-definite observed information, its inverse approximates the covariance matrix:

$$\widehat{\mathrm{Cov}}(\hat\theta) \approx \left[-\nabla^2\ell(\hat\theta)\right]^{-1}.$$

The diagonal entries give approximate variances, and their square roots give approximate standard errors. This is where the shape of the likelihood becomes uncertainty: a sharply curved peak has large information and small standard error, while a flat peak has small information and large standard error.

The quadratic approximation below shows why curvature is used as a local summary of information. Near the peak, the log-likelihood and its second-order approximation almost coincide; farther away, the approximation drifts from the exact curve.

For the Bernoulli model with $0<k<n$, the observed information at $\hat p=k/n$ is

$$J(\hat p)=\frac{n}{\hat p(1-\hat p)},$$

where:

• $J(\hat p)$: the observed information for the Bernoulli probability
• $n$: the number of trials
• $\hat p$: the estimated success probability
• $1-\hat p$: the estimated failure probability

so the approximate standard error is

$$\widehat{\mathrm{se}}(\hat p) = \sqrt{\frac{\hat p(1-\hat p)}{n}}.$$

where:

• $\widehat{\mathrm{se}}(\hat p)$: the approximate standard error of the Bernoulli estimate
• $\hat p$: the estimated success probability
• $1-\hat p$: the estimated failure probability
• $n$: the number of trials

This is the familiar standard error of a sample proportion. For the Gaussian mean with known variance, the observed information is $n/\sigma^2$, so the standard error is $\sigma/\sqrt{n}$. For the Gaussian variance, the observed information at $\hat v=\mathrm{SSE}/n$ is $n/(2\hat v^2)$. The approximate standard error is therefore $\hat v\sqrt{2/n}$, where $\hat v=\hat\sigma^2$ is the variance estimate.

When no closed-form solution exists, you can maximize the log-likelihood numerically. A simple one-dimensional approach is a grid search: evaluate the log-likelihood on a fine grid of parameter values and choose the grid point with the largest value. Grid search is easy to reason about but becomes impractical in more than a few dimensions.

A more scalable approach is Newton-Raphson. It comes from a second-order Taylor expansion of the objective around the current iterate, and it approximately solves the score equation $\nabla\ell(\theta)=0$ at each step. To maximize $\ell(\theta)$, you update

$$\theta_{t+1} = \theta_t - \left[\nabla^2\ell(\theta_t)\right]^{-1} \nabla\ell(\theta_t).$$

where:

• $\theta_t$: the parameter value at iteration $t$
• $\theta_{t+1}$: the updated parameter value
• $\nabla\ell(\theta_t)$: the gradient (score) of the log-likelihood at $\theta_t$
• $\nabla^2\ell(\theta_t)$: the Hessian matrix of second derivatives at $\theta_t$
• $[\cdot]^{-1}$: the matrix inverse

The gradient points in the direction of increase, and the inverse Hessian rescales that direction. The update is equivalent to solving $\nabla^2\ell(\theta_t)(\theta_{t+1}-\theta_t) = -\nabla\ell(\theta_t)$, which makes the score approximately zero after the step. When the Hessian is negative definite and the score is nonzero, the Newton direction is an ascent direction, but the full step can overshoot and decrease the likelihood. A line search can shorten the step. If the starting value is poor or the Hessian is not negative definite, Newton-Raphson can diverge. Practical implementations often minimize the negative log-likelihood with quasi-Newton methods, add bounds, or reparameterize positive quantities using logarithms. For example, instead of optimizing $\sigma^2$ directly, you can optimize $\log\sigma$ and exponentiate the result. This keeps the standard deviation positive.

A model is identifiable if different parameter values produce different distributions for the data. Formally, if

$$P_{\theta_1}(A)=P_{\theta_2}(A) \quad\text{for every measurable event }A$$

implies $\theta_1=\theta_2$, then the parameter is identifiable. Here $P_\theta(A)$ is the model probability of an event $A$. If distinct parameters give the same distribution, no amount of data from that model can distinguish those equivalent values. A continuous equivalence can create a flat ridge, while a discrete label symmetry can create separate equivalent peaks. A mixture model with unlabeled components is a classic example: swapping the labels of two components gives the same likelihood. In linear regression, perfectly collinear predictor columns make $X^\top X$ singular, so different coefficient vectors can produce the same fitted values. We discussed the practical consequences of that in Multicollinearity.

Identifiability is a property of the model and the data design, not a property of the optimizer. A numerical optimizer may return a point on a flat ridge, but the data do not distinguish that point from the other equivalent points. Before trusting an estimate, ask whether the data contain enough variation to distinguish the parameters. If not, we need to simplify the model, collect more informative data, or impose a constraint.

Sampling uncertainty is different from identifiability. Even in an identifiable model, the maximum likelihood estimate is a random variable because it depends on a random sample. For a correctly specified identifiable model with an interior true parameter, sufficient smoothness, nonsingular Fisher information, a consistent maximizing solution, and the derivative-convergence conditions used in the asymptotic argument, we obtain:

$$\sqrt{n}\left(\hat\theta-\theta_0\right) \xrightarrow{d} N\left(0,I_1(\theta_0)^{-1}\right),$$

where:

• $\sqrt{n}$: the square root of sample size, which scales the centered estimation error to a nondegenerate limit under the stated regularity conditions
• $\hat\theta$: the maximum likelihood estimate
• $\theta_0$: the true parameter value
• $\xrightarrow{d}$: convergence in distribution as $n\to\infty$
• $N(0,I_1(\theta_0)^{-1})$: the limiting normal distribution with mean zero and covariance $I_1(\theta_0)^{-1}$
• $I_1(\theta_0)$: the per-observation Fisher information evaluated at the true parameter $\theta_0$

This is the basis for the curvature-based standard errors above. The approximation can be poor in small samples, near boundaries, or when the model is misspecified. Alternatives include profile likelihood, a bootstrap appropriate to the sampling design, and Bayesian inference with an explicit prior. Each requires its own assumptions and coverage or sensitivity checks; none automatically repairs a boundary or a misspecified model.

Maximum likelihood also has failure modes that are worth recognizing. The likelihood can be unbounded, as in a Gaussian mixture when a component variance collapses to zero around a single point. In logistic regression, complete separation can push coefficients toward infinity because the likelihood keeps increasing as the separating boundary becomes sharper. Non-concave likelihoods can have multiple local maxima, so different starting values may produce different answers. And if the model is wrong, the estimate may be precise but precisely wrong.

Maximum likelihood has several properties that make it useful beyond the basic examples. The first is invariance under reparameterization. If $\hat\theta$ is the maximum likelihood estimate of $\theta$, and $g$ is a one-to-one function, then $g(\hat\theta)$ is the maximum likelihood estimate of $g(\theta)$. For example, if $0<\hat p<1$, then $\hat p/(1-\hat p)$ is the finite maximum likelihood estimate of the odds. This property is convenient, but it does not mean that standard errors transform by simply plugging in the transformed estimate. Uncertainty usually requires the delta method or a profile-likelihood calculation.

The second property is additivity for independent observations. The log-likelihood is a sum of per-observation log-density or log-mass contributions, so each new independent observation adds one term. This makes the likelihood scalable and explains why the curvature typically grows with sample size. In regular models, information adds across independent observations. For a correctly specified regular i.i.d. model with fixed, finite, positive definite per-observation information, total information is proportional to $n$, and the usual large-sample standard errors shrink at rate $1/\sqrt{n}$.

Under regularity conditions, maximum likelihood estimators are consistent and asymptotically efficient. Consistency means the estimate converges to the true parameter as the sample size grows. Asymptotic efficiency means attaining the inverse-Fisher-information bound within the regular estimator class. It does not compare every possible consistent estimator at every parameter value. These are large-sample guarantees. They do not imply that the estimate is unbiased in a small sample, that the model is correct, or that the likelihood is well behaved in every finite dataset.

Finally, maximum likelihood is not the same as least squares in general. For independent Gaussian errors and any fixed common positive variance, maximizing over the mean parameters is equivalent to minimizing squared error. Bernoulli, Poisson, and other non-Gaussian likelihoods generally give different objectives, although their maximizing estimates can sometimes coincide with least-squares estimates. This is why the likelihood framework is more general than squared error, and why the choice of likelihood encodes assumptions about the data-generating process.

A likelihood curve shows how the plausibility of the observed data changes as we move the parameter. The peak is the maximum likelihood estimate. The width of the curve tells you how much information the data provide about the parameter. To compare curves with different sample sizes, we often shift each log-likelihood so that its maximum is zero. This removes an arbitrary additive constant and makes the shapes easier to compare.

For Bernoulli samples with the same interior observed proportion, the log-likelihood is concave and becomes narrower as sample size grows. The peak stays at the sample proportion, but the data become more informative, so small departures from the peak cost more log-likelihood. That is the visual counterpart of the standard error formula $\sqrt{\hat p(1-\hat p)/n}$.

The plotted Gaussian log-likelihood is a downward parabola in $\mu$, because the variance is held fixed. Its horizontal peak is the estimate and its curvature controls the local precision. Exponentiating the unshifted log-likelihood would give a bell-shaped likelihood, but neither plot by itself assigns posterior probabilities to $\mu$. Bayesian parameter probabilities also require a prior and normalization.

We now work through the likelihood calculations for two small examples. The first example is a Bernoulli problem with a count. The second is a Gaussian problem with five observations. In both cases, the same pattern appears. Build the likelihood, take its logarithm, differentiate, solve the score equation, then inspect curvature.

Suppose twenty visitors arrive, and twelve click. We model each visit as an independent Bernoulli trial with unknown click probability $p$. The observed count is $k=12$ and the number of trials is $n=20$. The likelihood is

$$L(p)=p^{12}(1-p)^{8}.$$

where:

• $L(p)$: the likelihood of the click probability $p$ given the observed counts
• $p$: the unknown click probability
• the exponents $12$ and $8$: the numbers of clicks and non-clicks, so the terms carry the success and failure counts

The log-likelihood is

$$\ell(p)=12\log p+8\log(1-p).$$

The score is

$$\ell'(p)=\frac{12}{p}-\frac{8}{1-p}.$$

Setting the score to zero gives

$$\frac{12}{p}=\frac{8}{1-p}.$$

Cross-multiplying gives $12(1-p)=8p$, so $12-12p=8p$, and therefore $12=20p$. Hence

$$\hat p=\frac{12}{20}=0.6.$$

The second derivative is

$$\ell''(p)=-\frac{12}{p^2}-\frac{8}{(1-p)^2}.$$

At $\hat p=0.6$, this is

$$\begin{aligned} \ell''(0.6) &= -\frac{12}{0.6^2} - \frac{8}{(1-0.6)^2} && \text{(substitute } \hat p = 0.6 \text{)} \\ &= -\frac{12}{0.36} - \frac{8}{0.16} && \text{(square the denominators)} \\ &= -33.\overline{3} - 50 && \text{(divide each term)} \\ &= -83.\overline{3} && \text{(add)}. \end{aligned}$$

where:

• $\ell''(0.6)$: the second derivative of the Bernoulli log-likelihood evaluated at $\hat p = 0.6$
• $\frac{12}{0.6^2}$: the magnitude of the negative second-derivative contribution from the 12 successes
• $\frac{8}{(1-0.6)^2}$: the magnitude of the negative second-derivative contribution from the 8 failures

The observed information is therefore $J(0.6) = 83.\overline{3}$, so the approximate standard error is

$$\widehat{\mathrm{se}}(\hat p) =\sqrt{\frac{1}{83.333}} \approx 0.1095.$$

where:

• $\widehat{\mathrm{se}}(\hat p)$: the approximate standard error of the Bernoulli estimate
• $83.333$: the observed information $J(0.6)$ rounded

This matches the formula $\sqrt{\hat p(1-\hat p)/n}=\sqrt{0.6\cdot0.4/20}=\sqrt{0.012}\approx0.1095$.

We can compare the log-likelihood at a few candidate values. At $p=0.5$,

$$\begin{aligned} \ell(0.5) &= 12\log(0.5) + 8\log(0.5) && \text{(substitute } p = 0.5 \text{)} \\ &= 20\log(0.5) && \text{(combine like terms)} \\ &\approx -13.863 && \text{(evaluate)}. \end{aligned}$$

At $p=0.6$,

$$\begin{aligned} \ell(0.6) &= 12\log(0.6) + 8\log(0.4) && \text{(substitute } p = 0.6 \text{)} \\ &\approx -13.460 && \text{(evaluate)}. \end{aligned}$$

At $p=0.7$,

$$\begin{aligned} \ell(0.7) &= 12\log(0.7) + 8\log(0.3) && \text{(substitute } p = 0.7 \text{)} \\ &\approx -13.912 && \text{(evaluate)}. \end{aligned}$$

The largest value is at $p=0.6$, as expected. The differences in log-likelihood are modest, which is another way of seeing that $p=0.5$ is not absurdly implausible with only twenty observations. A rough 95% Wald interval is $0.6\pm1.96(0.1095)$, or approximately $(0.385,0.815)$. The interval is symmetric around $0.6$ and can be unreliable near the boundaries, but it gives a useful first summary of uncertainty.

Suppose five readings from a sensor are

$$y=(9.8,\ 10.3,\ 9.9,\ 10.2,\ 9.8).$$

We model them as independent draws from $N(\mu,\sigma^2)$. The sample mean is

$$\hat\mu=\bar y=\frac{9.8+10.3+9.9+10.2+9.8}{5}=10.0.$$

where:

• $\hat\mu$: the maximum likelihood estimate of the mean
• $\bar y$: the sample mean
• the numerator: the sum of the five observed readings
• the denominator $5$: the number of observations $n$

The deviations from the mean are $-0.2,0.3,-0.1,0.2,-0.2$. Squaring and summing gives

$$\begin{aligned} \mathrm{SSE} &= (y_1-\bar y)^2 + (y_2-\bar y)^2 + (y_3-\bar y)^2 + (y_4-\bar y)^2 + (y_5-\bar y)^2 && \text{(definition)} \\ &= (-0.2)^2 + (0.3)^2 + (-0.1)^2 + (0.2)^2 + (-0.2)^2 && \text{(substitute deviations)} \\ &= 0.04 + 0.09 + 0.01 + 0.04 + 0.04 && \text{(square each term)} \\ &= 0.22. && \text{(sum)} \end{aligned}$$

where:

• $\mathrm{SSE}$: the sum of squared deviations of the readings from the sample mean
• each squared term: one observation's deviation from the mean, squared

The maximum likelihood variance is therefore

$$\hat\sigma^2=\frac{\mathrm{SSE}}{n}=\frac{0.22}{5}=0.044.$$

The maximum likelihood standard deviation is

$$\hat\sigma=\sqrt{0.044}\approx0.2098.$$

where:

• $\hat\sigma$: the maximum likelihood standard deviation
• $0.044$: the maximum likelihood variance estimate

For comparison, the unbiased sample variance divides by $n-1=4$:

$$\begin{aligned} s^2 &= \frac{\mathrm{SSE}}{n-1} = \frac{0.22}{4} = 0.055, \\ s &= \sqrt{s^2} = \sqrt{0.055} \approx 0.2345. \end{aligned}$$

where:

• $s^2$: the unbiased sample variance, which divides by $n-1$ rather than $n$
• $s$: the square root of the unbiased sample variance; $s$ itself is generally biased for $\sigma$
• $n-1 = 4$: the degrees of freedom, one less than the number of observations

The maximum likelihood variance is smaller because it divides by $n$. On average, it underestimates the true variance by a factor of $(n-1)/n=4/5$. With only five observations, the bias is not negligible.

The bias is easiest to see as a ratio. The maximum likelihood variance has expected value $(n-1)\sigma^2/n$, while the unbiased sample variance has expected value $\sigma^2$. The gap closes as $n$ increases.

The score for $\mu$ at the estimate is

$$\begin{aligned} \left.\frac{\partial\ell}{\partial\mu}\right|_{\mu=10} &= \frac{1}{\hat\sigma^2}\sum_{i=1}^{n}(y_i-10) && \text{(score formula at } \mu=10\text{)} \\ &= \frac{1}{0.044}\left(-0.2+0.3-0.1+0.2-0.2\right) && \text{(substitute deviations)} \\ &= \frac{1}{0.044}(0) && \text{(sum the deviations)} \\ &= 0 && \text{(simplify)}. \end{aligned}$$

where:

• $\left.\frac{\partial\ell}{\partial\mu}\right|_{\mu=10}$: the score evaluated at $\mu = 10$
• $\hat\sigma^2$: the maximum likelihood variance estimate, $0.044$
• $y_i-10$: each observation's deviation from the candidate mean

The score is zero at $\hat\mu = \bar y = 10$, which is exactly the first-order condition for the maximum.

The approximate standard error of the mean, using $\hat\sigma^2$ as the variance, is

$$\begin{aligned} \widehat{\mathrm{se}}(\hat\mu) &= \sqrt{\frac{\hat\sigma^2}{n}} && \text{(information } n/\hat\sigma^2 \text{)} \\ &= \sqrt{\frac{0.044}{5}} && \text{(substitute values)} \\ &\approx \sqrt{0.0088} && \text{(divide)} \\ &\approx 0.0938. && \text{(square root)} \end{aligned}$$

where:

• $\widehat{\mathrm{se}}(\hat\mu)$: the approximate standard error of the estimated mean
• $\hat\sigma^2$: the maximum likelihood variance estimate
• $n$: the number of observations

A rough 95% interval using the normal approximation is $10.0\pm1.96(0.0938)$, or approximately $(9.816,10.184)$. With only five observations, a $t$ interval based on four degrees of freedom would be wider, which is a reminder that the curvature approximation is asymptotic. The likelihood still gives us the point estimate and a sense of precision, but the finite-sample interval should be chosen with care.

You will now implement the same calculations from primitives. First you compute the closed-form estimates directly. Then you compute a likelihood grid and run numerical optimization with SciPy. Finally, we compute curvature-based standard errors. The data are small and deterministic so that every printed value can be checked by hand.

The next cell imports SciPy's optimizer. The earlier hidden setup already imported NumPy, Matplotlib and the normal distribution and seeded the shared plot environment.