Newton Raphson Algorithm

Let f(x) be a function (possibly multivariate) and suppose we are interested in determining
the maximum of f and, often more importantly, the value of x which maximizes f. The most
common statistical application of this problem is finding a Maximum Likelihood Estimate (MLE).
This document discusses the Newton Raphson method

1 Motivation

Newton Raphson maximization is based on a Taylor series expansion of the function f(x). Specifically,
if we expand f (x) about a point a,

where f'(·) is the gradient vector and f''(·) is the hessian matrix of second derivatives. This creates
a quadratic approximation for f. We know how to maximize a quadratic function (take derivatives,
set equal to zero , and solve)

The Newton Raphson process iterates this equation. Specifically, let be a starting point for
the algorithm and define successive estimates recursively through the equation

If the function f(x) is quadratic , then of course the quadratic “approximation” is exact and the
Newton Raphson method converges to the maximum in one iteration. If the function is concave,
then the Newton Raphson method is gauranteed to converge to the correct answer. If the function
is convex for some values of x , then the algorithm may or may not converge. The NR algorithm may
converge to a local maximum and not the global maximum, it might converge to a local minimum,
or it might cycle between two points . Starting the algorithm near the global maximum is the best
practical method for helping convergence to the global maximum.

Fortunately, loglikelihoods are typically approximately quadratic (the reason asymptotically
normality occurs for many random variables ). Thus, the NR algorithm is an obvious choice for
finding MLEs. The starting value for the algorithm is often a “ simpler ” estimate (in terms of ease
of computation) of the parameter, such as a method of moments estimator.


Let and suppose we want the joint maximum likelihood estimate
of . The loglikelihood is

Solving for the MLE analytically requires solving the equations

These two equations cannot be solved analytically (the Gamma function is difficult to work with).
The two equations do provide us with the gradient

The hessian matrix is

The starting values of the algorithm may be found using the method of moments. Since
and , the method of moment estimators are and .
Suppose we have data (in truth actually generated with α = 2 and β = 3) such that n = 1000,
, and s2 = 0.2235679. The algorithm begins at
2.073976 and = 3.04577. After 3 iterations, the NR algorithm stabilizes at and

Note that there is no need to treat this situation as a multivariate problem. The partial
derivative for β may be solved in terms of α.

Thus, for each value of α, the maximum of β is attained at . Thus, we can reduce the
to the one-dimensional problem of maximizing

This is called the profile loglikelihood. The first and second derivatives are

Again starting the algorithm at 2.073976, the algorithm converges to after
???? iterations.

Prev Next