Solving Likelihood Score Equations Numerically

Why Score Equations Need Numerical Methods

For some severity models, the likelihood equations give parameter estimates directly. For a gamma model with both shape and rate unknown, the shape estimate generally does not have a closed-form expression. It must be found by solving a nonlinear equation. The calculation is manageable if the likelihood is first simplified and the numerical root is checked rather than accepted on appearance alone.

This tutorial uses complete, independent, positive observations and the gamma parameterization below. The task is specifically to solve the shape score equation; numerical optimization is not a substitute for identifying the correct likelihood and score.

Definition: Suppose $X_1,\ldots,X_n$ are independent and identically distributed gamma random variables with shape $\alpha>0$ and rate $\beta>0$. Their density is $$ f_X(x;\alpha,\beta)=\frac{\beta^\alpha}{\Gamma(\alpha)}x^{\alpha-1}e^{-\beta x},\qquad x>0. $$ Here $E[X]=\alpha/\beta$, and $\operatorname{Var}(X)=\alpha/\beta^2$. The function $\psi(\alpha)=\frac{d}{d\alpha}\log\Gamma(\alpha)$ is the digamma function.

Deriving the Shape Equation

Write $\bar{x}=\frac{1}{n}\sum_{i=1}^n x_i$ and $\overline{\log x}=\frac{1}{n}\sum_{i=1}^n\log x_i$. Taking the logarithm of the joint density gives the complete-data log likelihood

\ell(\alpha,\beta) =n\alpha\log\beta-n\log\Gamma(\alpha) +(\alpha-1)\sum_{i=1}^n\log x_i-\beta\sum_{i=1}^n x_i.

For fixed $\alpha$, differentiating with respect to $\beta$ and setting the derivative equal to zero gives $n\alpha/\beta-\sum_i x_i=0$, so $\widehat{\beta}(\alpha)=\alpha/\bar{x}$. Substituting this profile estimate into the shape score yields the equation below.

Result: Define $h(\alpha)=\log\alpha-\psi(\alpha)$ and $c=\log\bar{x}-\overline{\log x}$. The gamma shape score equation is $$ h(\alpha)=c. $$ After solving for $\widehat{\alpha}$, the rate estimate is $\widehat{\beta}=\widehat{\alpha}/\bar{x}$. These results assume positive, complete observations from the stated gamma model; censored or truncated observations require different likelihood contributions.

The equation follows directly by differentiating the profile log likelihood. In detail, substituting $\beta=\alpha/\bar{x}$ into the log likelihood and differentiating with respect to $\alpha$ gives

\frac{d\ell_p}{d\alpha} =n\left(\log\alpha-\psi(\alpha)+\overline{\log x}-\log\bar{x}\right) =n\left(h(\alpha)-c\right).

Thus a root of $h(\alpha)-c$ is exactly a zero of the profiled shape score. The quantity $c$ is nonnegative by the arithmetic-geometric mean inequality: $\bar{x}\geq \exp(\overline{\log x})$. When observations are not all equal, $c>0$.

Choosing and Checking a Root-Finding Method

A robust approach is to first find a bracket: two positive values at which $h(\alpha)-c$ has opposite signs. Bisection then repeatedly halves that interval and cannot leave it. Newton's method can converge faster, but its update must remain in the positive parameter space:

\alpha_{\text{new}} =\alpha-\frac{h(\alpha)-c}{h'(\alpha)}, \qquad h'(\alpha)=\frac{1}{\alpha}-\psi_1(\alpha),

where $\psi_1$ is the trigamma function. A useful safeguard is to reject a Newton step that is nonpositive or outside a known bracket and use a bisection step instead. Stop only when the root estimate is sufficiently stable and the score residual $U(\alpha)=n(h(\alpha)-c)$ is suitably small. A rounded parameter value can have a larger residual than the unrounded value used internally.

There is at most one finite root when $c>0$. Indeed, the standard series for the trigamma function gives $\psi_1(\alpha)=\sum_{k=0}^{\infty}(\alpha+k)^{-2}$. Since this sum is greater than $\int_0^\infty(\alpha+t)^{-2}\,dt=1/\alpha$, it follows that $h'(\alpha)<0$. Also, $h(\alpha)$ decreases from infinity toward zero as $\alpha$ increases. Therefore every positive $c$ has a unique positive solution.

Worked Example: Newton's Method for Four Hypothetical Losses

A hypothetical portfolio has four complete gamma-severity observations, measured in thousands of dollars: $1,2,3,4$. Estimate the shape and rate. Use the gamma shape equation and check the score residual.

Method and summaries. The observations are positive and complete, so the stated gamma likelihood applies. Their mean is $\bar{x}=(1+2+3+4)/4=2.5$. Their mean log is $\overline{\log x}=\log(24)/4\approx0.794513458$, while $\log\bar{x}=\log(2.5)\approx0.916290732$. Hence $c\approx0.121777274$.

Bracket. Numerical evaluation gives $h(4)\approx0.13017669$ and $h(5)\approx0.10332024$. In particular, $\psi(5)=1+1/2+1/3+1/4-\gamma\approx1.506117668$, so $h(5)=\log 5-\psi(5)\approx0.103320244$. Since $h(4)>c>h(5)$, the unique root is between 4 and 5.

Newton refinements. Starting at $\alpha_0=4$, the derivative is $h'(4)=1/4-\psi_1(4)\approx-0.03382296$. The first update is

\alpha_1 =4-\frac{0.13017669-0.121777274}{-0.03382296} \approx4.24834.

Evaluating the same update again gives approximately $4.2654$; one further refinement gives $\widehat{\alpha}\approx4.26543$. The unrounded solution has $h(\widehat{\alpha})$ approximately $0.1217773$, matching $c$. Using the displayed five-decimal estimate gives a total shape-score residual with absolute value less than $0.0001$.

Rate and check. The estimated rate, in inverse thousands of dollars, is $\widehat{\beta}=4.26543/2.5\approx1.70617$. The fitted mean is $\widehat{\alpha}/\widehat{\beta}=2.5$, agreeing with the sample mean as required by the profiled rate equation. The common mistake avoided here is using an inaccurate endpoint value or claiming convergence merely because the parameter updates look small; the score residual provides a direct check.

Bracketing and a Second Data Set

Bisection is valuable when a Newton step is unstable, when a derivative is difficult to evaluate accurately, or when a transparent numerical certificate is preferred. Each midpoint replaces the bracket endpoint having the same sign of $h(\alpha)-c$. The root remains enclosed at every iteration.

Worked Example: Bracketing a Shape Estimate for Two Losses

For a second hypothetical portfolio, two complete losses are $1$ and $10$, in thousands of dollars. Find a bracket for the shape, obtain a numerical estimate, and calculate the rate.

Form the equation. Here $\bar{x}=(1+10)/2=5.5$, and $\overline{\log x}=(\log1+\log10)/2\approx1.151292546$. Thus $c=\log5.5-\overline{\log x}\approx1.704748092-1.151292546=0.553455546$.

Refine the root. Since $h(1)=\gamma\approx0.57721566>c$ and $h(1.1)\approx0.51907<c$, the root lies in $(1,1.1)$. For a faster safeguarded calculation, Newton's method from $\alpha_0=1$ uses $h'(1)=1-\pi^2/6\approx-0.644934$, giving

\alpha_1 =1-\frac{0.57721566-0.553455546}{-0.644934} \approx1.03684.

At $1.03684$, numerical evaluation gives $h(\alpha_1)\approx0.55437$, still above $c$. The next Newton refinement is approximately $1.03835$, and a further refinement gives $\widehat{\alpha}\approx1.03838$. This remains inside the original bracket. At the displayed estimate, the total score residual has absolute value less than $0.00001$.

Calculate the rate. The estimate is $\widehat{\beta}=1.03838/5.5\approx0.188796$ per thousand dollars. The rate is positive, and $\widehat{\alpha}/\widehat{\beta}=5.5$, which checks the profiled mean. Newton's method is useful here, but the bracket remains the safeguard if an iterate fails to behave well.

Boundary Behavior and a Third Example

When every observation has the same value, $c=0$. Since $h(\alpha)>0$ for every finite positive $\alpha$, the score equation then has no finite root. A root-finding program cannot fix this by extending its search: the data have led to a boundary behavior, not a missed finite solution. This distinction is important when diagnosing a numerical failure.

Worked Example: A Nearly Boundary-Adjacent Data Pattern

A hypothetical set of four complete losses, in thousands of dollars, is $1,1,4,4$. Determine the shape equation and solve it numerically. This example also illustrates that repeated values do not by themselves imply the boundary case; all observations must be equal for $c=0$.

Calculate the summaries. The arithmetic mean is $\bar{x}=(1+1+4+4)/4=2.5$. The mean log is $(0+0+\log4+\log4)/4=\log2\approx0.693147181$. Therefore $c=\log2.5-\log2=\log1.25\approx0.223143551>0$, so there is a unique finite shape root.

Numerical solution. At $\alpha=2.4$, $h(2.4)\approx0.222567$, just below the target. The derivative there is approximately $-0.09878$. One Newton update gives

\alpha_{\text{new}} =2.4-\frac{0.222567-0.223143551}{-0.09878} \approx2.3942.

Refining the calculation gives $\widehat{\alpha}\approx2.3942$. At the four-decimal estimate, the absolute total score residual is less than $0.0001$. The rate is $\widehat{\beta}=2.3942/2.5\approx0.95768$ per thousand dollars, and the fitted mean is again 2.5. This method avoids treating a small score discrepancy at an intermediate trial value as evidence that no root exists.

Common Numerical Traps

Reversing the target. The equation is $h(\alpha)=\log\bar{x}-\overline{\log x}$, not the reverse. Reversing the difference can produce an impossible negative target for unequal positive observations.
Mixing gamma rate and scale. This tutorial uses the rate $\beta$, so $\widehat{\beta}=\widehat{\alpha}/\bar{x}$. The reciprocal is the scale, not the rate.
Trusting a Newton step without checking it. Newton's method is local. A step can leave the positive domain or a known bracket; use a safeguarded step instead.
Using only the update size as a convergence check. A small update does not necessarily mean the score is small. Evaluate $U(\widehat{\alpha})=n[h(\widehat{\alpha})-c]$ at the reported estimate.
Searching for a finite root when all observations are equal. In that case $c=0$, while $h(\alpha)>0$ for finite $\alpha$. Report the boundary behavior rather than an artificial numerical root.

Key takeaway: Profile out the gamma rate, solve the one-dimensional equation $h(\alpha)=c$, and verify the result with a score residual. A sign-changing bracket makes the numerical search dependable; a positive target ensures a unique finite shape root.

Check Your Understanding

For each problem, use the stated gamma shape equation and distinguish a numerical root from a boundary case.

A hypothetical complete sample is $2,2,8,8$, measured in thousands of dollars. Its scale is twice that of the $1,1,4,4$ sample above. Give the shape estimate and rate estimate. (Answer: $\widehat{\alpha}\approx2.3942$, $\widehat{\beta}\approx0.4788$ per thousand dollars.)
For four observations with $\bar{x}=2.5$ and $\overline{\log x}=\log(24)/4$, state the target $c$, a bracket for the root, and the approximate shape estimate. (Answer: $c\approx0.121777274$, bracket $(4,5)$, $\widehat{\alpha}\approx4.26543$.)
A numerical calculation gives $c=0.110000$. Using $h(4)\approx0.13017669$ and $h(5)\approx0.10332024$, identify a valid bracket and state why it contains the unique root. (Answer: $(4,5)$; $h(4)>c>h(5)$, and $h$ is strictly decreasing.)
All six complete observations equal $7$, in the chosen loss units. What is $c$, and does the finite shape score equation have a root? (Answer: $c=0$; no finite positive root, because $h(\alpha)>0$ for every finite $\alpha>0$.)
Every loss in the sample $1,2,3,4$ is multiplied by 3 to express the data in a larger unit. State how the gamma shape estimate and rate estimate change. (Answer: the shape remains approximately $4.26543$; the rate is divided by 3, becoming approximately $0.568724$ per new unit.)

What You'll Learn