Score Equations and Numerical Maximum Likelihood

From a Likelihood to a Score Equation

A likelihood turns observed data into a function of unknown model parameters. Maximizing that function is often easier after taking its logarithm: the logarithm preserves the location of the maximum and converts products into sums. The score is the derivative of the log-likelihood with respect to a parameter. Its zeros are useful candidates for a maximum, but a zero alone does not establish that the candidate is an admissible maximum.

Definition: For observed data \(x\) with likelihood \(L(\theta;x)\), the log-likelihood is \(\ell(\theta;x)=\log L(\theta;x)\). For a scalar parameter \(\theta\), the score is \(U(\theta)=\partial\ell(\theta;x)/\partial\theta\). For a parameter vector \(\boldsymbol{\theta}\), the score is the vector of partial derivatives, \(\nabla\ell(\boldsymbol{\theta};x)\).

For independent observations with densities or probability mass functions \(f(x_j;\theta)\), the likelihood is \(L(\theta;x)=\prod_{j=1}^n f(x_j;\theta)\), and the log-likelihood is \(\ell(\theta;x)=\sum_{j=1}^n\log f(x_j;\theta)\). The same principle applies to a likelihood built for a different observation scheme, but the score must be derived from that scheme’s likelihood; a complete-data formula should not be used unchanged for censored or otherwise incomplete observations.

Interior maximum condition: If the parameter is scalar, the log-likelihood is differentiable, and a maximum occurs at an interior point of the allowed parameter space, then its score is zero there. For a vector parameter, the gradient is zero at an interior maximum. These conditions are necessary, not sufficient: a root could be a minimum, a saddle point, or one of several candidates.

Why Differentiating the Log-Likelihood Works

Suppose \(L(\theta)>0\) and is differentiable. Since the logarithm is strictly increasing, maximizing \(L\) and maximizing \(\ell=\log L\) gives the same parameter value. Also,

\frac{d\ell}{d\theta}=\frac{1}{L(\theta)}\frac{dL}{d\theta}.

At an interior stationary point, \(dL/d\theta=0\) therefore implies \(d\ell/d\theta=0\), and conversely. The second derivative helps classify a scalar stationary point: a negative value of \(d^2\ell/d\theta^2\) indicates a local maximum, while a positive value indicates a local minimum. Even a local maximum still needs to be compared with other candidates and any allowed boundary points.

Write the likelihood for the observed data.
Include the sampling assumptions and any known exposure or observation restrictions that define the data.

Take logarithms and differentiate.
Keep track of which terms depend on the unknown parameter; constants may be omitted only after the likelihood is correctly specified.

Find and check candidates.
Solve the score equations when possible, check curvature or compare log-likelihood values, and evaluate admissible boundaries.

Enforce the parameter space.
A numerical root outside the permitted parameter range is not an admissible estimate. Optimize over the allowed set and compare feasible candidates.

Worked Examples with Analytic Score Equations

Worked Example: A Constrained Bernoulli Probability

A hypothetical portfolio records whether each of 12 independent trial policies has a claim during a fixed observation period. Eight policies have a claim. Suppose each claim indicator is Bernoulli with probability \(p\), and the model restricts \(0.20\leq p\leq0.75\). Find the maximum likelihood estimate of \(p\).

Method. The data contain 8 successes and 4 failures, so the likelihood, apart from no omitted parameter-dependent terms, is \(L(p)=p^8(1-p)^4\). Its log-likelihood and score are

\ell(p)=8\log p+4\log(1-p),\qquad U(p)=\frac{8}{p}-\frac{4}{1-p}.

Set the score equal to zero and solve:

\frac{8}{p}=\frac{4}{1-p} \quad\Longrightarrow\quad 8(1-p)=4p \quad\Longrightarrow\quad p=\frac{8}{12}=\frac{2}{3}.

The candidate \(2/3\) lies in the allowed interval. The second derivative is \(\ell''(p)=-8/p^2-4/(1-p)^2\), which is negative throughout \(0<p<1\); thus the log-likelihood is strictly concave and this feasible stationary point is the unique maximum. The estimate is \(0.667\) to three decimals.

Check and trap. The estimate is the observed success fraction, and \(8/12=2/3\) is within the stated bounds. A common mistake is to use a boundary estimate simply because the parameter is constrained. Here neither boundary is the maximum: the interior score root is feasible. If the constraint had excluded \(2/3\), the feasible boundary would need to be checked instead.

Worked Example: Joint Maximum Likelihood Estimates for a Normal Model

Five hypothetical repair costs, measured in thousands of dollars, are \(1,2,3,4,5\). Assume independent observations from \(N(\mu,\sigma^2)\), with both \(\mu\) and \(\sigma^2>0\) unknown. Find the joint maximum likelihood estimates.

Method. The likelihood is the product of five normal densities. Writing \(v=\sigma^2\) to keep the variance parameter explicit, the log-likelihood is

\ell(\mu,v)=-\frac{5}{2}\log(2\pi)-\frac{5}{2}\log v -\frac{1}{2v}\sum_{j=1}^5(x_j-\mu)^2,\qquad v>0.

Differentiate with respect to the mean:

\frac{\partial\ell}{\partial\mu} =\frac{1}{v}\sum_{j=1}^5(x_j-\mu)=0 \quad\Longrightarrow\quad \hat{\mu}=\frac{1+2+3+4+5}{5}=3.

For the variance, first use the mean equation and then differentiate:

\frac{\partial\ell}{\partial v} =-\frac{5}{2v}+\frac{1}{2v^2}\sum_{j=1}^5(x_j-\mu)^2=0 \quad\Longrightarrow\quad \hat v=\frac{1}{5}\sum_{j=1}^5(x_j-\hat\mu)^2.

The deviations from 3 are \(-2,-1,0,1,2\), and their squared values are \(4,1,0,1,4\). Their sum is \(10\), so \(\hat v=10/5=2\) in squared thousands of dollars. The fitted standard deviation is \(\sqrt{2}\) thousand dollars. The variance estimate is positive, as required.

Check and trap. The residuals sum to zero, as required by the mean score equation, and their squared sum is 10. The maximum likelihood variance divides by \(n=5\), not by \(n-1=4\); dividing by \(n-1\) gives a different estimator and does not solve the likelihood score equation. The fitted variance is \(2\), in squared units, not a standard deviation of \(2\).

Numerical Maximum Likelihood and Constraints

Some score equations cannot be rearranged to give a closed-form estimate. A numerical method can then find a root or maximize the log-likelihood directly. One common root-finding method is Newton iteration. For a scalar score \(U(\theta)\), it updates a current value \(\theta_k\) by

\theta_{k+1}=\theta_k-\frac{U(\theta_k)}{U'(\theta_k)} =\theta_k-\frac{\ell'(\theta_k)}{\ell''(\theta_k)}.

Newton iteration can converge quickly near a suitable root, but it can also diverge, move outside the parameter space, or converge to a root that is not a maximum. A safeguarded method such as bisection can be more reliable when a score changes sign over a known interval. For constrained optimization, a numerical optimizer should respect the parameter bounds; then compare its result with any other feasible stationary points and boundaries. In multiple dimensions, the analogous curvature information is the Hessian matrix of second partial derivatives.

Worked Example: A Numerical Score Root with an Upper Constraint

For a hypothetical normalized claim severity \(x\) between 0 and 1, suppose the density is \(f_\theta(x)=\theta e^{\theta x}/(e^\theta-1)\) for \(0<x<1\) and \(\theta>0\), with the uniform density as the limiting model at \(\theta=0\). The model restricts \(0\leq\theta\leq1.5\). Three observed severities are \(0.40000, 0.70000, 0.86955\). Find the constrained maximum likelihood estimate. (The values are normalized units, not dollars.)

Method. The sample size is \(n=3\), and the sum is \(0.40000+0.70000+0.86955=1.96955\). For \(\theta>0\), the log-likelihood and score are

\ell(\theta)=3\log\theta+1.96955\theta-3\log(e^\theta-1), \qquad U(\theta)=\frac{3}{\theta}+1.96955-\frac{3e^\theta}{e^\theta-1}.

The score equation is equivalent to matching the sample mean to the model mean:

\bar{x}=m(\theta),\qquad m(\theta)=\frac{1}{1-e^{-\theta}}-\frac{1}{\theta}.

Here \(\bar{x}=1.96955/3=0.6565167\). There is no elementary rearrangement that isolates \(\theta\) in this equation. To locate the unconstrained score root, use Newton iteration on \(m(\theta)-\bar{x}=0\), where \(m'(\theta)=1/\theta^2-e^{-\theta}/(1-e^{-\theta})^2\). Starting at \(\theta_0=1\), \(m(1)=0.58198\) and \(m'(1)=0.07933\), giving

\theta_1=1-\frac{0.58198-0.6565167}{0.07933}\approx1.9396.

A second iteration gives approximately \(2.000\), the unconstrained score root. This value is outside the allowed interval, so it cannot be the constrained estimate. To check the allowed interval, \(m(\theta)\) is increasing because \(m'(\theta)>0\) for \(\theta>0\); also \(m(1.5)\approx0.62055\). Thus for \(0<\theta\leq1.5\), \(U(\theta)=3\{\bar{x}-m(\theta)\}>0\). The log-likelihood is increasing throughout the permitted interval, so its maximum occurs at the upper boundary:

\hat{\theta}=1.5.

Check and trap. The unconstrained score root is not the answer under the specified model constraint. Reporting approximately \(2.000\) would violate \(0\leq\theta\leq1.5\). The positive score throughout the feasible interval verifies that moving toward its upper endpoint increases the log-likelihood; the constrained estimate is therefore \(1.5\).

Practical Checks Before Reporting an Estimate

Score equations organize the calculation, but a complete maximum likelihood solution also checks that the reported estimate is valid. In particular, a numerical routine’s output should not be accepted merely because it stopped or returned a number. Check that the estimate satisfies the model’s parameter restrictions and that the score is close to zero when the estimate is an interior solution. For a boundary solution, a nonzero score may be appropriate; the direction of change within the feasible parameter space must support the maximum.

Situation	What to check
Interior scalar estimate	The score is zero (or numerically close to zero), and the candidate is a maximum rather than a minimum.
Several score roots	Evaluate the log-likelihood at feasible candidates and compare them.
Constrained parameter	Check admissible boundaries as well as feasible interior roots.
Numerical estimate	Check convergence, parameter bounds, and whether small changes in the estimate materially change the log-likelihood.

A useful numerical alternative to enforcing a positive parameter bound directly is to reparameterize, for example by writing a positive parameter as \(e^\eta\) for unrestricted \(\eta\). This prevents invalid negative values, but it does not remove the need to check convergence or boundary behavior. Some boundaries cannot be represented by a finite \(\eta\), and a constrained maximum may occur exactly there.

Key takeaway: A score equation supplies candidate estimates, not an automatic answer. Derive it from the correct log-likelihood, solve it analytically or numerically, then verify the candidate’s curvature, feasibility, and likelihood relative to other candidates and boundaries.

Check Your Understanding

For each problem, derive or use the score as appropriate, and pay attention to the stated parameter space.

Ten independent hypothetical indicators contain 6 successes. Under a Bernoulli model with \(0.1\leq p\leq0.9\), what is the maximum likelihood estimate of \(p\)? (Answer: \(0.6000\))
For independent observations from \(N(\mu,4)\), where the variance \(4\) is known, three observed values are \(1,3,5\). Find the maximum likelihood estimate of \(\mu\). (Answer: \(3\))
For independent observations from \(N(\mu,\sigma^2)\) with both parameters unknown, the observed values are \(2,2,6,6\). Find the joint maximum likelihood estimates of \(\mu\) and \(\sigma^2\). (Answer: \(\hat{\mu}=4,\ \hat{\sigma}^2=4\))
A differentiable log-likelihood on \(1\leq\theta\leq4\) has score \(U(\theta)=6-2\theta\). Find its constrained maximum likelihood estimate. (Answer: \(\hat{\theta}=3\))
A differentiable log-likelihood on \(0\leq\theta\leq2\) has score \(U(\theta)=5-\theta\), which is positive throughout the permitted interval. Where is the constrained maximum? (Answer: \(\hat{\theta}=2\), the upper boundary)

What You'll Learn