Computing Maximum-Likelihood Estimates for 

Parameters of the National Descriptive Model 

of Mercury in Fish

Open-File Report 2012–1181

U.S. Department of the Interior

U.S. Geological Survey
Computing Maximum-Likelihood Estimates for 

Parameters of the National Descriptive Model 

of Mercury in Fish

By David I. Donato

Open-File Report 2012–1181

U.S. Department of the Interior

U.S. Geological Survey
ii

U.S. Department of the Interior

KEN SALAZAR, Secretary

U.S. Geological Survey

Marcia K. McNutt, Director

U.S. Geological Survey, Reston, Virginia: 2012

For more information on the USGS—the Federal source for science about the Earth,

its natural and living resources, natural hazards, and the environment—visit 

http://www.usgs.gov or call 1–888–ASK–USGS

For an overview of USGS information products, including maps, imagery, and publications,

visit http://www.usgs.gov/pubprod

Suggested citation:

Donato, D.I., 2012, Computing maximum-likelihood estimates for parameters of the National Descriptive Model of 

Mercury in Fish: U.S. Geological Survey Open-File Report 2012–1181, 16 p., at http://pubs.usgs.gov/of/2012/1181.

Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the 

U.S. Government.

Although this information product, for the most part, is in the public domain, it also may contain copyrignted materials 

as noted in the text. Permission to reproduce copyrighted items must be secured from the copyright owner.
iii

Contents

1 Introduction ............................................................................................................................................................ 1

2 Maximum-Likelihood Estimates for a Linear Model................................................................................................ 3

2.1 A Linear Model................................................................................................................................................ 3

2.2 The Likelihood Function for the Linear Model ................................................................................................. 3

2.3 The Log-Likelihood Function for the Linear Model .......................................................................................... 6

2.4 Maximizing the Likelihood and Log-Likelihood Functions for the Linear Model............................................... 6

3 Maximum-Likelihood Estimates for the National Descriptive Model of Mercury in Fish.......................................... 7

3.1 Contributions to the Likelihood from Censored Observations ......................................................................... 8

3.2 The Likelihood Function for the National Descriptive Model of Mercury in Fish .............................................10

3.3 The Log-Likelihood Function for the National Descriptive Model of Mercury in Fish ......................................11

3.4 Maximizing the Likelihood and Log-Likelihood Functions for the National Descriptive Model of Mercury in

Fish ................................................................................................................................................................12

4 Numerical Computation of Parameter Estimates for the National Descriptive Model of Mercury in Fish...............14

4.1 Newton-Raphson Iteration .............................................................................................................................14

5 Summary...............................................................................................................................................................15

References Cited..........................................................................................................................................................16

Appendix.......................................................................................................................................................................16
1

Computing Maximum-Likelihood Estimates for Parameters 

of the National Descriptive Model of Mercury in Fish

By David I. Donato

1 Introduction

The National Descriptive Model of Mercury in Fish Tissue (NDMMF) is a statistical model 

(Wente, 2004) used to predict the concentration in fish tissue of methylmercury, a potent neurotoxin and 

a known health hazard for humans and wildlife. The utility of the NDMMF derives from its power to 

explain much of the observed variation in fish-tissue methylmercury concentrations as (1) variation by 

geographic location, (2) variation over time, and (3) variation by species and fish length due to 

bioaccumulation. Because it is based on a national database of methylmercury observations, this model 

has the potential to be used locally throughout the conterminous United States and in parts of Canada 

for planning ongoing sampling of methylmercury in fish, and for preparing public health warnings (fish 

consumption advisories) identifying specific species and sizes of fish from particular locations known or 

likely to be high in methylmercury content (Hearn and others, 2006).

The NDMMF may be expressed formally as follows:

1 1

n m

ijk j j k k ijk

j k

y vg  

= =

=++ ∑ ∑  (1)

where

1: if the event is 

0 : if the event is not ,

ln( 1),

 the parameter for

 the parameter for t

 the th sampling e

he th species/cut combination, and

ln(

ve

1) : if t

t

h

n ,

j

k

ijk

k

ijk ijk

j

j

j

y C

v

k

length

g

j







=





=

=

+

+

=

=

e species/cut is 

 0 : if the species/cut is not 

k

k







and where

• Cijk is the observed methylmercury concentration for the ith sample from sampling event j with 

species/cut combination (or type) k ,

• n is the number of sampling events and m is the number of species/cut combinations,

• ijk length is the length in inches of the fish for the ith sample from sampling event j with 

species/cut combination (or type) k ,

• some of the { } Cijk may be left censored, and
2

• ijk  is a random variable specifying the error for the ith sample from sampling event j with 

species/cut combination (or type) k .

The term sampling event as it applies to the NDMMF means the collection of samples from a 

particular location (such as a lake or part of a stream) at one particular time. The term species/cut in this 

context refers to the combination of a particular species of fish and a type of laboratory treatment 

(Wente, 2004). In this publication, the notations “ ln( ) a ” or “ ln[a] ” are used for the natural logarithm

of a and are equivalent to “ log ( ) e a ”.

Before the NDMMF can be used to predict methylmercury concentrations in fish tissue, its 

parameters (the sets of values { }  j and { } k ) must be estimated on the basis of a dataset of 

observations, ordinarily for thousands of sampling events and hundreds of species/cut combinations. 

Among many statistical approaches that might be used to compute a set of parameters for this model in 

order to fit the parameters to the observed data, two widely used methods are those of least-squares and 

maximum-likelihood estimation (Cramér, 1946; Miller and Miller, 2004). Because some observations 

are left censored (that is, some are at or below some detection limit), the method of maximumlikelihood estimation is more suitable for estimating parameters for the NDMMF because it makes 

better use of the information provided by censored observations than the method of least-squares 

estimation (Helsel, 2004).

Computing the maximum-likelihood estimates of the parameters for the NDMMF requires

substantial computer-processor time and memory because the NDMMF defines an unusually large 

number of parameters and variables (tens of thousands). Consequently, the task of fitting the parameters 

of the NDMMF to a growing database of fish-tissue mercury observations now strains or exceeds the 

computational capabilities of readibly available statistical software packages and computer 

workstations. This report addresses this problem by deriving and presenting the mathematical and 

statistical results required to enable the development of custom computer software capable of fitting the 

parameters of the NDMMF more quickly and using less computer memory than general-purpose 

statistical software.

The detailed mathematical explanation of the application of the method of maximum likelihood 

to the NDMMF in this publication includes the derivation of the mathematical expressions used in 

computing maximum-likelihood estimates of the parameters of the NDMMF. These details are included 

in order to do the following:

1. Provide mathematical and computational expressions suitable for use in 

implementing custom software for fitting parameters to the NDMMF in the 

maximum-likelihood sense; and

2. Present the derivation of maximum-likelihood-estimation (MLE) expressions for the

current NDMMF in sufficient detail to facilitate future derivation of revised MLE 

expressions in case the NDMMF is revised.

To aid understanding of the derivation of the MLE expressions for the NDMMF, this publication 

first presents in section 2 the derivation of MLE expressions for a much simpler but structurally 

similar model.1 Then section 3 derives and presents the mathematical MLE expressions for the 

 

 1This report uses the symbols α and β differently than the seminal report on the NDMMF (Wente, 2004) in order to clarify 

and emphasize the structural similarity between a simple linear statistical model and the NDMMF. The linear model 

described in section 2 illustrates concepts in a simple setting; these concepts are then applied to the more complex NDMMF 

in section 3.
3

NDMMF. Finally, section 4 discusses the actual computation of maximum-likelihood estimates for 

NDMMF parameters using the results of section 3.

2 Maximum-Likelihood Estimates for a Linear Model

This section applies the method of maximum-likelihood estimation to a particular linear 

statistical model involving left-censored data. The derivation of maximum-likelihood estimates for this 

simple model illustrates the basic techniques and principles to be applied in the following section for the 

more complicated model—the NDMMF.

2.1 A Linear Model

The particular linear statistical model used for illustration in this section is:

y a bx =+ +  , (2)

where

•  ~ (0, ) N  ,

• some observed values of y may be left censored and, thus, known only to be at or

below a detection limit that may differ from observation to observation,

2

• y is the response variable,

• x is the independent variable, and

• a, b, and σ are the unknown parameters to be estimated.

2.2 The Likelihood Function for the Linear Model

In general, a likelihood function is the joint probability or probability density of a set of 

observations (that is, of a sample of observed values of the independent and response variables) 

construed as a function of the unknown parameters of a statistical model; in general, finding the 

maximum-likelihood estimates of the unknown parameters means finding a set of values for the 

unknown parameters that maximizes this joint probability (or probability density) of the observations 

under the model. Like least-squares estimation, maximum-likelihood estimation fits the values of the 

parameters of a statistical model to a set of observed values of independent variables and response 

variables (Cramér, 1946; Miller and Miller, 2004).

For this particular linear model, the likelihood function (,, ) a b  is the joint probability density 

of the uncensored observations multiplied by the joint probability of the censored observations, with this 

entire product regarded as a function of the unknown parameters a, b, and σ. Thus, the likelihood 

function for this linear model is:

( ) 1

2 2 2

1

2

1

1 ( ) (,, ) exp { |,, , }

2

n m i i

jj j

i j

y a bx a b   P y D ab x = =  

    − +

= − ⋅≤    

   

 ∏ ∏ , (3)

where

 

 2That is, if j y is a left-censored observation and Dj is its associated detection limit, then j j y D≤ .
4

• { |,, , } P y D ab x jj j ≤  is the probability that j y is less than or equal to Dj given 

particular values for a, b, σ, and j x ,

• Dj is the detection limit for censored observation j,

• n is the number of uncensored observations, and

• m is the number of left-censored observations.

The likelihood function is an application of the general multiplicative formula of probability

theory. The general multiplicative formula states that when each member of a set of observations is 

independently drawn from a population with a known probability distribution, then the joint probability 

(or joint probability density) of the set of observations is the product of the probabilities (or probability 

densities) of the individual observations (Chung, 1974; Wilks, 1962). The likelihood function

“ ( , , )”  a b  defined here follows this principle in that it is the product of probability densities for all n 

uncensored observations and of probabilities for all m censored observations.

To understand the first repetitive product on the right-hand side of equation 3, first recall that the 

probability density function (PDF) of a normal distribution3 is

2

1

, 2 2 2

1 () ( ) exp

2

w f w  



 

  − = −   

, (4)

where w is the independent variable and μ and σ are the mean and standard deviation of the normal 

distribution, respectively (Miller and Miller, 2004).

Then observe that equation 2 can be rewritten as

y a bx −+ = ( )  (5)

to emphasize that the difference between the observed response y and the predictor ( ) a bx + is 

equivalent to the error term  . Since  is a normally distributed random variable with  = 0 , equation 5 

indicates that y a bx − + ( )is also a stochastic quantity following the normal probability distribution with 

a mean of zero. Therefore,

( y a bx y a bx − + −= − + ( )0 ( ) ) (6)

can be substituted for ( ) w−  in equation 4 to express the probability density for each uncensored 

observation (, ) i i x y as a function of the unknown parameters(,, ) a b  . The probability density thus 

derived,

( ) 1

2 2 2

2

1 ( )

exp

2

i i y a bx

 

  − +   −

 

, (7)

has the properties required for the likelihood function. This probability density tends to be relatively 

large when the error term ( )

2

( ) i i y a bx − + is near zero, and it tends to be relatively small when the error 

term is further from zero. This means that choices of (,, ) a b  that produce small error terms (errors of 

estimation near zero) also produce higher values of the likelihood function. Intuitively, higher likelihood 

corresponds to a "better" fit of the model to the data.

 

 3Although it is conventional to use the variable x in expressing the PDF and cumulative distribution function (CDF) of the 

normal distribution, w is used here to avoid confusion with the specific censored observations{ }i x .
5

To understand the second repetitive product on the right-hand side of equation 3, consider what 

it means for an observation to be left censored. An observation is left censored if we know that its 

value lies somewhere on the real number line (or perhaps on the nonnegative side of the real number 

line) at or below a particular value Dj , generally referred to as the detection limit4 (Helsel, 2004).

Thus, for the m left-censored observations, the multiplicative contribution of observation j to the 

likelihood is the probability that observation j is at or below the detection limit Dj . Since  ~ (0, ) N  , 

we can express the contribution of left-censored observations to the likelihood in terms of the 

cumulative distribution function (CDF) of the normal distribution. Let F0, denote the CDF of a normal 

distribution with mean μ = 0 and standard deviation σ. 

2

1

0, 2 2 2

1 ( 0) ( ) exp

2

A

w F A  dw

  −∞

  − = −   ⌠



⌡

(8)

denotes the probability of observing a value less than or equal to some specific value A.

Since

( ) j j jj j j j j y D a bx D ≤ → + +≤ → ≤ −+   D a bx , (9)

then

( ) 2

1

2 2 2

1 { | , , , } exp .

2

D a bx j j

jj j

w P y D a b x  dw

 

− +

−∞

  ≤= −   

⌠



⌡ (10)

Thus, equation 3 can be rewritten as follows:

( ) ( ) 2

1 1

2 2 2 2 2 2

1

2

1

1 1 ( ) ( , , ) exp exp . 2 2

D a bx j j n m i i

i j

y a bx w a b  dw

   

− +

= = −∞

    − +     =− ⋅ −                

⌠



⌡

 ∏ ∏ (11)

For many physical measurements, values are constrained to be nonnegative. Laboratory 

measurements of concentrations of chemicals in samples of organic material are an example of 

nonnegative values that may be subject to detection limits. A measurement that is both below a 

detection limit and also constrained to be nonnegative can be described as interval censored. An 

observation is interval censored if we know that its value lies in an interval on the real number line 

between an upper and a lower limit. The term left censored is, however, widely applied to censored 

observations of concentrations; therefore, when the term left censored is applied to concentrations, the 

restriction to nonnegative values is implicit. When the response variable j y for the linear model 

represents a concentration value, its multiplicative contribution to the likelihood differs from that of an

unconstrained left-censored observation. To determine the contribution to likelihood when j y is 

 

 4A detection limit is both a lower bound and an upper bound. It is a lower bound for meaningful point measures. When a 

measurement procedure (such as a laboratory analysis of a chemical concentration) cannot provide a measurement below a 

certain value, that value is the detection limit. In such a case, a measurement below the detection limit has no meaning as a 

point value, but it does still provide quantitative information that can be used in statistical analysis: the measurement conveys 

the information that a value has been observed that is less than or equal to the detection limit. The detection limit is also the 

upper bound for the unknown value that could not be measured as a point value. A measurement at or below the detection 

limit is referred to as a left-censored observation.
6

constrained to be nonnegative, this constraint on j y must be used to find the constraint on the error term 

 . As

0 0 () j j j j j y ≥ → + + ≥ → ≥− + a bx   a bx , (12)

then in the case of censored concentrations,

( ) 2

1

2 2 2

( )

1 {0 | , , , } exp

2

j j

j

D a bx

jj j

a bx

w P y D a b x  dw

 

− +

− +

  ≤≤ = −   ⌠



⌡ (13)

is the contribution to the likelihood for each observation.

For simplicity, the following two subsections on the linear model treat censored observations as

general left-censored observations that may potentially have negative values.

2.3 The Log-Likelihood Function for the Linear Model

Since the logarithm function is defined for all positive real numbers and is monotone increasing 

on the positive side of the real line, the logarithm of the likelihood function, or the log-likelihood 

function, will be maximized by the same parameter set as the likelihood function. For some models, a 

third option of minimizing the negative log-likelihood function will find the same parameter set and 

may prove to be either analytically or computationally easier than either of the other two equivalent 

methods of maximum-likelihood estimation for the model. Because of this equivalence of optimization 

problems, we can choose to solve whichever is easiest. In general, numerical solutions for maximumlikelihood problems tend to use either the positive or negative log-likelihood function—especially for 

models (like the NDMMF) that involve large numbers of parameters—because machine computations 

involving logarithms are less prone to the problems of overflow, underflow, and loss of precision.

For the linear model with some left-censored observations, the log-likelihood function is the 

following:

( ) 2 ( ) 0, ( )

2 1 1

1 1 2

ln ( , , ) ln ( ) ln ( ) . 2 (2 )

n m

n i i j j

n i j

a b y a bx F D a bx   

  = =

 

= − −+ + −+          

 ∑ ∑ (14)

Equation 14 is derived from equations 8 and 11. Notice that the natural logarithm of an indexed 

product is an indexed sum. The log-likelihood function is easier to use than the likelihood function itself 

in finding the optimal parameter set for the linear model.

2.4 Maximizing the Likelihood and Log-Likelihood Functions for the Linear Model

In order to maximize the log-likelihood function (and the likelihood function as well), we must 

find the ordered triple ˆ (,, ) a b ˆ ˆ , if it exists, defined so that for any ordered triple (,, ) a b′′ ′  it will be 

true that ( ) ( ) ˆ ln ( , , ) ln ( , , )   ab a b ˆ   ˆ ≥ ′′ ′ , or equivalently, ˆ   (,, ) ( , , ) ab a b ˆ   ˆ ≥ ′′ ′ .

If a maximizing ordered triple ˆ (,, ) a b ˆ ˆ exists, it will satisfy the three simultaneous conditions 

determined by equating to zero (0) each of the three partial derivatives of the log-likelihood function 

with respect to a, b, and 2  (Miller and Miller, 2004). The three simultaneous equations that must be 

satisfied by ˆ (,, ) a b ˆ ˆ are:
7

( ) ( )

( ) 2

1 1 0,

1 ( ) (ln[ ( , , )]) ( ) 0

( )

n m j j

i i

i j j j

f D a bx a b y a bx

a F D a bx







 = =

∂ − +

= −+ − = ∂ − +

 ∑ ∑ , (15)

( ) ( )

( ) 2

1 1 0,

1 ( ) (ln[ ( , , )]) ( ) ( ) 0 ( )

n m j j

i ii j

i j j j

f D a bx a b y a bx x x

b F D a bx







 = =

∂ − +

= −+ − = ∂ − +

 ∑ ∑ , and (16)

( ) ( ( ) )

( )

2 0,

2 2 4

1 1 0,

2 ( ) 1 (ln[ ( , , )]) ( ) 0 .

2 2 ( )

n m j j

i i

i j j j

d F D a bx

n d a b y a bx F D a bx





 

   = =

− + ∂ −

= + −+ + = ∂ − +

 ∑ ∑ (17)

In each of equations 15 and 16, the contribution to the partial derivative of the log-likelihood 

function from uncensored observations equals the contribution from the left-censored observations. 

Equation 17 showing the partial derivative of the log-likelihood function with respect to 2  could be 

replaced by an equation involving the partial derivative with respect to σ; either choice is acceptable 

analytically, but the equation involving the partial derivative with respect to 2  is preferable for 

computation.

In equation 17, the expression

2 ( 0, ( j j ( )) ) d F D a bx

d  

− +

denotes the derivative of the normal CDF with respect to 2  evaluated at the point ( ) D a bx j j − + ;

in this expression, the normal CDF is construed as a function of 2  . Thus,

( )

( ) 2

1

0, 2 2 2

1 ( ) exp

2

D a bx j j

j j

w F D a bx  dw

 

− +

−∞

 

−+ = −   

⌠



⌡ (18)

where the variable of integration is denoted as w in order to avoid notational confusion with the 

observed values of the independent variable j x . (Please see equation 4.) In solving the system of 

equations 15, 16, and 17 for a,b, and σ, this expression for the derivative of the normal CDF with 

respect to 2  must be evaluated numerically. The Newton-Raphson Method (also known as Newton’s 

Method) may be used to solve this system of three simultaneous equations numerically. The NewtonRaphson Method for finding maximum-likelihood parameters is described in section 4.1 (Weisstein, 

1999).

3 Maximum-Likelihood Estimates for the National Descriptive Model of Mercury in Fish

The methods used in section 2 to estimate the three unknown parameters of a simple linear 

model in the maximum-likelihood sense may now be specialized and applied to estimate the n m+ +1

parameters of the NDMMF in the maximum-likelihood sense. Comparison of the linear model (as 

follows)

y a bx =+ + 
8

with the NDMMF (as follows)

1 1

n m

ijk j j k k ijk

j k

y vg  

= =

=++ ∑ ∑ 

shows some structural similarity between the the linear model and the NDMMF but significant

dissimilarity as well. In particular, although the NDMMF is structurally linear, it is effectively nonlinear 

because the covariate term ln( 1) g length k = + ijk is a nonlinear function of the physically measured 

covariate (fish length) and because the response variable ln( 1) ijk ijk y C = + is a nonlinear function of the 

measured responseCijk .

Before discussing the fitting of the parameters of the NDMMF, it is instructive to consider how

the NDMMF—once its parameters have been determined—is used to estimate or predict a fish-tissue 

methylmercury concentration. The simple linear model can be used to compute a response y from a 

value for the single independent variable x along with values for the two parameters(,) a b . In 

comparison, the NDMMF also requires three items of data in order to compute an estimate of a fishtissue methylmercury concentration: the α parameter associated with a particular sampling event, the β

parameter associated with a particular species/cut combination, and a fish length (in inches). The linear 

model and the NDMMF are similar in that any single estimate makes use of exactly two parameters. An 

estimate made with the linear model makes use of the parameters a and b, while an estimate made with 

the NDMMF makes use of exactly one  j and exactly one k because, except for the selected event j

and the selected species/cut k, the indicator variables “ j v ” and “ k g ” are zero for all other sampling 

events and for all other species/cut combinations, respectively. Although NDMMF model parameters 

are used only two at a time for prediction, all parameters are involved in fitting the NDMMF to 

observed data in the maximum-likelihood sense. This is shown by the equations and systems of 

equations for maximizing the likelihood and log-likelihood functions of the NDMMF presented below.

3.1 Contributions to the Likelihood from Censored Observations

As previously explained in section 2.2, the term left censored is applied to observations known 

to be at or below some detection limit Dk . In some cases, an observation described as left censored may 

be subject to additional constraints, such as the restriction of chemical concentrations to be nonnegative. 

The related term interval censored refers to an observation known to fall in an interval on the real 

number line between a lower limit and an upper limit. Left censoring is a special case of interval 

censoring with an implicit (rather than explicit) lower limit for the censoring interval. The implicit lower 

limit may depend on either, or both, of the form of the statistical model and the physical characteristics 

of the observed data. The response variable y for the linear model considered in section 2 is an example 

of a response constrained neither by the form of the model nor by the physical characteristics of the 

observations; this response variable y can assume any value on the real number line, so this model’s 

implied lower limit5 for a left-censored observation is −∞ . In contrast to the simple linear model, the 

response variable for the NDMMF is implicitly constrained both by the form of the model and by the 

physical characteristics of chemical concentrations. The functional form of the NDMMF response 

variable ln( 1) ijk ijk y C = + constrains the measured value Cijk to be strictly greater than −1 in order that 

the argument of the ln function also be strictly greater than zero (as it must be in order for the logarithm 

function to be defined for the argument). Ignoring for the moment that the concentration Cijk must be 

 

 5As this example illustrates, the interval applied in interval censoring may be an improper interval.
9

nonnegative, left censoring for the NDMMF is equivalent to interval censoring with −< ≤ 1 C D ijk k . 

(The detection limit is designated here by Dk because it can vary from observation to observation.)

For the sake of generality, and also because a number of studies have estimated NDMMF 

parameters assuming interval censoring with −< ≤ 1 C D ijk k , the likelihood and log-likelihood functions 

for the NDMMF are expressed in the following subsections in terms of this least restrictive form of 

interval censoring. These generalized expressions can be readily specialized by constraining 

concentrations to be nonnegative, as explained in the final paragraph of section 3.2 below.

The likelihood function for the NDMMF presented in the following section 3.2 depends in part 

on the probabilities of the censored observations. Determining the probability of a censored observation 

requires that the upper and lower limits for measured concentrations be used in finding the 

corresponding upper and lower limits on the stochastic error term ijk  . The upper and lower limits for 

ijk  determine the probability of a censored observation as the integral of the PDF of the error term ijk 

from the minimum error value to the maximum error value; or, equivalently, as the CDF of the error 

term evaluated at its upper limit minus the CDF of the error term evaluated at its lower limit, as shown 

below in equation 21 in section 3.2.

Throughout the remainder of section 3, let U Dc denote the upper limit for a censored 

concentration and let DL denote the lower limit.

6 The upper and lower limits for the error term ijk  can 

be determined from the upper and lower bounds U Dc and DL for censored concentrations as follows:

Finding the upper bound for ijk  :

ln( 1) 1 1 ijk ijk y y

ijk ijk ijk ijk y C e C Ce = + → =+→ =− .

Therefore

1 1 ln 1 . ( ) ijk ijk U UU U y y CD e D e D y D ijk c ≤ → −≤ → ≤ + → ≤ + c c ijk c

Substituting c c c c ijk   v g + +  for ijk y in the inequality above results in

ln 1 ( ) ln 1 ( ) . ( ) U U

c c c c ijk c ijk c c c c c   vg D + + ≤ + → ≤ +− +   D vg  

Taking the natural logarithm of both sides of the inequality is valid (preserves the inequality) because 

the natural logarithm function is monotone increasing over its entire domain.

Finding the lower bound for ijk  :

 1 1 ln 1 ( ) . ijk ijk y y CD e D e D y D ijk L ≥ → −≥ → ≥ + → ≥ + L L ijk L

So, therefore, substituting c c c c ijk   v g + +  for ijk y as before results in

c c c c ijk ln 1 ( L ) ijk ln 1 ( ) . ( L ) c c c c   vg D + + ≥ + → ≥ +− +   D vg  

 

 6 This notation indicates that the upper censoring limit may vary from observation to observation; however, the lower limit 

is the same for all censored observations.
10

3.2 The Likelihood Function for the National Descriptive Model of Mercury in Fish

In its general form, the NDMMF does not necessarily restrict the form of the various random 

error terms{ } ijk  . For purposes of determining the likelihood function for the NDMMF, however, the 

model will be restricted by the simplifying assumption that all members of the set of random error terms 

{ } ijk  are independently and identically distributed as normal variates with mean  = 0 and standard 

deviation σ (where σ is an unknown parameter of the NDMMF), as follows:

{ } are i.i.d. with ~ (0, ) for all { }   ijk ijk N  ijk .

With this simplifying assumption about the error terms, the likelihood function for the NDMMF 

may be expressed as:

({ },{ }, ) ·    j k uc = P P , (19)

where

2

1

2 2

1

y = y ,

,

,

,

1 ( ) exp , 2

ln( 1), and

the number of uncensored observations. 

u n

u uu u u

u

u

u j

u k

u ijk

u

u ijk

u j

yvg P

g length

n

v v

 

  

 

 

=

  − − = −   

=

=

= +

=

=

∏

(20)

and

( ) ( [ ] )

( )

1

0, 0,

1

t

{ |{ },{ }, },

ln 1 ln 1 ,

,

,

,

ln 1 ,

 he number of censored observations,

 the upper detection limit 

c

c

n

U

c L ijk c j k

c

n

U

c cc c c L cc c c

c

c j

c k

c j

c ijk

c

U

c

P PD C D

F D v gF D v g

v v

g length

n

D

 

 

   

 

 

=

=

= ≤≤

= +− −   +− −   

=

=

=

= +

−

=





=



∏

∏

for the censored observation , and

 the lower limit for the censoring interval for the censored observation . L

c

D c =

(21)

In equation groups 20 and 21, each index u or c identifies a particular uncensored or censored 

observation, respectively, which is associated with exactly one event j and exactly one species/cut k.
11

In the derivation of the log-likelihood function, it is convenient to express Pu in the following 

alternative but equivalent form:

[ ] 2

2

2

1

1 1

exp .

(2 ) 2

u

nu u

n

u u uu u u n u

P yvg  

   =

    = − −−         ∑ (22)

The likelihood function for the NDMMF is the product of two composite factors contributing to 

the likelihood: a contribution from uncensored observations and another contribution from censored 

observations. The contribution Pu from uncensored observations is a product of probability densities, 

and the contribution Pc for censored observations is a product of probabilities. Each factor in the 

indexed product Pu is the probability density for a particular observation of an uncensored mercury 

concentration and fish length in association with a specific event j and a specific species/cut k for any 

given values of the parameters{ }  j , { } k , and σ. Each factor in the indexed product Pc is the 

probability that a censored mercury concentration was in the range of possible values below the 

detection limit U Dc and simultaneously above the lower censoring limit DL for any given values of the 

parameters.

Equation 21 can be readily specialized for two different interpretations of censored observations 

of concentrations. Case 1: As stated previously in the first paragraph of section 3.1, the functional form 

of the NDMMF requires thatCijk > −1; this is equivalent to requiring that the lower limit of the 

censoring interval, DL , approach (but not reach) −1 from above. As 1 DL → − , the term 

0, (ln[ 1] ) F D vg  L cc c c +− −   in equation 21 approaches 0, F  ( )0 −∞ = because as 1 0 DL + → , then

ln[ 1] DL + → −∞ . Therefore, when 1 DL → − , the term 0, (ln[ 1] ) F D vg  L cc c c +− −   can be dropped 

from equation 21. Case 2: If Cijk is constrained to be nonnegative, then 0 DL = and

ln[ 1] ln[1] 0. DL += = In this latter case, the second term in equation 21 becomes 0, ( ). F vg  − −   cc c c

3.3 The Log-Likelihood Function for the National Descriptive Model of Mercury in Fish

As explained previously for the linear model presented in section 2.1, the likelihood function 

will be maximized by the same parameter set that maximizes the log-likelihood function.

7 In this sense, 

the problem of maximizing the likelihood is equivalent to the problem of maximizing the log-likelihood. 

For the NDMMF, as for many maximum-likelihood estimation problems, it is advantageous both 

analytically and computationally to work with the log-likelihood function rather that the likelihood 

function.

For the NDMMF, the log-likelihood function is:

( )

[ ]

2

2

1 2

0, 0,

1

ln ({ },{ }, ) ln( ) ln( )

1 1 ln

2 (2 )

ln (ln[ 1] ) (ln[ 1] ) .

u

u

u

c

j k u c

n

n u uu u u

n u

n

U

c cc c c L cc c c

c

P P

yvg

F D v gF D v g  

 

 



 

   

=

=

= +

    = − −−          

+   +− − − +− −  

∑

∑



(23)

 

 7Although a maximizing parameter set can be found for reasonably consistent datasets, there are datasets for which a 

maximizing parameter set may not exist or may be difficult to compute using conventional fixed-precision machine numbers. 

This potential problem is beyond the scope of this publication.
12

3.4 Maximizing the Likelihood and Log-Likelihood Functions for the National Descriptive Model 

of Mercury in Fish

To maximize the likelihood and log-likelihood functions for the NDMMF, it is necessary to 

determine a collection of n m+ +1 parameters ˆ {{ },{ }, } ˆ ˆ   j k such that

ˆ ({ },{ }, ) ({ },{ }, ) ˆ ˆ     j k j k

′ ′′  ≥ (24)

for every ordered (n m+ +1)-tuple of parameters({ },{ }, )   j k

′ ′′ or equivalently such that

ˆ ln ({ },{ }, ) ln ({ },{ }, ) ˆ ˆ     j k j k

′ ′′   ≥         (25)

for every ordered(n m+ +1)-tuple of parameters({ } , ,{ } )   j k

′ ′′ .

If a maximum for the likelihood and log-likelihood functions exists, it will satisfy the n m+ +1

simultaneous conditions determined by equating to zero each of the n m+ +1partial derivatives of the 

log-likelihood function with respect to each of the n  j’s, the m k ’s, and 2  (Miller and Miller, 2004).
13

The partial derivatives of the log-likelihood function have three forms, as shown below:

1. For the  j’s:

( )

( ) 2

1

0, 0,

1 0, 0,

ln ( { },{ }, )

1

(ln[ 1] ) (ln[ 1] )

(ln[ 1] ) (ln[ 1] )

0 .

u

c

j k

j

n

u uu u u

u

n U

c cc c c L cc c c

U

c c cc c c L cc c c

yvg

f D v gf D v g

F D v gF D v g

 

 

 



 



   

   

=

=

∂   =   ∂

− −

  − +− − + +− −   +  +− − − + −   − 

=

∑

∑



(26)

The summations in the jth equation implied by equation 26 only pick up values for the 

observations for event j .

2. For the k ’s:

( )

( ) 2

1

0, 0,

1 0, 0,

ln ( { },{ }, )

1

(ln[ 1] ) (ln[ 1] )

(ln[ 1] ) (ln[ 1

0 .

] )

u

c

j k

k

n

u uu u u u

u

n U

c cc c c L cc c c

U c

c c cc c c L cc c c

y v gg

f D v gf D v g

g F D v gF D v g

 

 

  

 



   

   

=

=

∂   =   ∂

− −

− +− − + +− −

+

+− − − +− −

   

   

=

∑

∑



(27)

The summations in the kth equation implied by equation 27 only pick up values for the observations 

involving species/cut k .

3. For 2  :

2 ( )

2

2 4

1

2 0, 0,

1 0, 0,

ln ( { },{ }, )

1 ( ) 2 2

(ln[ 1] ) (ln[ 1] )

0 .

(ln[ 1] ) (ln[ 1] )

u

c

j k

n

u

u uu u

u

U n c cc c c L cc c c

U

c c cc c c L cc c c

n

yvg

d F D v gF D v g d

F D v gF D v g

 

 

 



 

 

   



   

=

=

∂   =   ∂

−+ − −



  

+− − − +− −   

+ =

+− − − +− − 

∑

∑



(28)

The summations in equation 28 pick up values from all observations for all sampling events and for all 

species/cut combinations.
14

4 Numerical Computation of Parameter Estimates for the National 

Descriptive Model of Mercury in Fish

The simple linear model of section 2 illustrates in principle what is involved in computing 

maximum-likelihood estimates of the parameters for a statistical model. Equations 15, 16, and 17

constitute a system of three simultaneous equations in the unknown parameters (,, ) a b  . The three 

parameters can be determined by solving the system of three equations. In principle, model parameters 

for the NDMMF are computed as they are for the linear model, by solving a system of simultaneous 

equations. The computational problem and task for the NDMMF is, however, on a very different scale 

from that of the linear model. Equations 26, 27, and 28 imply a system of n m+ +1 equations in as 

many unknown parameters with n equaling the number of sampling events and m equaling the number 

of species/cut combinations. In practical cases, this means that finding the maximum-likelihood 

estimates of the parameters of the NDMMF requires solving a system of tens of thousands of equations 

in as many unknowns.

Since there is no analytical solution for the system of n m+ +1 equations, the system of 

equations must in practice be solved numerically. A numerical solution will be computationally feasible 

only if the system of equations converges from an initial set of parameter estimates to a solution within 

a reasonably small number of iterations. The Newton-Raphson Method is an obvious choice for this 

problem because it generally offers relatively rapid convergence in finding the zeroes of functions.

4.1 Newton-Raphson Iteration

The Newton-Raphson Method (also called Newton’s Method) is a well-known numerical 

algorithm for finding the zero or root of a function or equation (Weisstein, 1999). The zero or root of a 

function of a real variable is the value of the independent variable at which the value of the function is 

equal to zero. The Newton-Raphson Method requires that the function be continuous and differentiable. 

In essence, the Newton-Raphson Method makes use of a few low-order terms of the Taylor Series 

representation of the function within a neighborhood of a root in order to find that root. If the function 

has few or no points of inflection in the neighborhood and is otherwise well behaved within the 

neighborhood of the root, the Newton-Raphson Method tends to converge rapidly to an acceptably 

precise solution within a few iterations.

Computing maximum-likelihood estimates of the numerous parameters of the NDMMF requires 

that the Newton-Raphson Method be applied repeatedly to each equation in the system of n m+ +1

equations implied by equations 26, 27, and 28. A single pass through all n m+ +1equations, in order to 

find a parameter for each equation through a separate Newton-Raphson iteration for each equation, does 

not suffice to determine the values of all parameters because of the dependencies among the equations 

in the system of n m+ +1equations. In the n equations implied by equation 26, the u and c

parameters depend on the current values of the u and c parameters. Similarly, in the m equations 

implied by equation 27, the u and c parameters depend on the current values of the u and c

parameters. In equation 28, the value of 2  depends on the values of all of the  and  parameters. 

Thus, as each successive Newton-Raphson iteration determines a value for one parameter, the values of 

many other parameters are affected and must be corrected through another pass. Multiple passes through 

the n m+ +1equations (applying the Newton-Raphson Method to each equation in sequence during each 

pass) are required in order to achieve convergence to the overall collection of parameter values that 

solves the system of n m+ +1equations and, thereby, computes the set of parameters that maximizes the 

likelihood and log-likelihood functions.
15

Despite the abundance of opportunity during this computation for the dependencies among the 

parameters to prevent convergence towards an overall solution, experience shows that convergence does 

in fact take place in practice when the starting set of initial parameter estimates is sufficiently close to 

the maximum-likelihood estimates and the data meet minimal requirements for consistency and 

completeness. Experience also shows that when the starting estimates of parameters used in a NewtonRaphson iteration are those determined by ordinary least-squares estimation, these estimates are usually 

sufficiently close to the maximum-likelihood estimates to enable convergence of the multiple-pass, 

sequential Newton-Raphson iteration described in this section.

5 Summary

This report explains how the method of maximum-likelihood estimation (MLE) can be applied 

to find the parameters of the National Descriptive Model of Mercury in Fish (NDMMF), a statistical 

model used to predict the concentration of methylmercury in fish tissue. Computing the parameters for 

the NDMMF in the maximum-likelihood sense is a complicated procedure because the NDMMF often 

involves a large number of parameters (tens of thousands) and because the NDMMF makes use of the 

information provided by left-censored observations of methylmercury concentrations. Because of the 

large number of parameters and because of the additional computation required to make use of leftcensored observations, the task of computing parameters for the NDMMF strains or exceeds the 

computational capabilities of readily available statistical software packages and computer workstations. 

Consequently, researchers and practitioners need custom software to fit NDMMF parameters in the 

MLE sense more quickly using less computer memory.

The two main purposes of this publication are the following:

1. To provide mathematical and computational expressions suitable for use in 

implementing custom software for fitting parameters to the NDMMF in the 

maximum-likelihood sense; and

2. To present the derivation of maximum-likelihood estimation (MLE) expressions for 

the current NDMMF in sufficient detail to facilitate derivation of revised MLE 

expressions for any future revisions of the NDMMF.

As an aid to understanding, this publication illustrates the mathematics of maximum-likelihood 

parameter estimation for a simple linear model before deriving the mathematical results for the more 

complex NDMMF. Finally, this publication explains in general terms how the Newton-Raphson Method 

can be used in computations to find the set of best-fit parameters for the NDMMF in the maximumlikelihood sense.
16

References Cited

Chung, Kai Lai, 1974, Elementary probability theory with stochastic processes: New York, SpringerVerlag, 325 p.

Cramér, Harald, 1946, Mathematical methods of statistics: Princeton University Press, 575 p.

Hearn, Paul, and others, 2006, EMMMA; a Web-based system for environmental mercury mapping, 

modeling, and analysis: U.S. Geological Survey Open-File Report 2006–1086, 13 p.

Helsel, D.R., 2004, Nondetects and data analysis; statistics for censored environmental data: WileyInterscience, 268 p.

Miller, Irwin, and Miller, Marylees, 2004, John E. Freund’s mathematical statistics with applications 

(7th ed.): Upper Saddle River, N.J., Pearson Prentice Hall, 614 p.

Weisstein, E.W., 1999, "Newton-Raphson Method" and "Newton’s Method" in CRC concise 

encyclopedia of mathematics: Chapman & Hall/CRC, 1969 p.

Wente, S.P., 2004, A statistical model and national data set for partitioning fish-tissue mercury 

concentration variation between spatiotemporal and sample characteristic effects: U.S. Geological 

Survey Scientific Investigations Report 2004–5199, 15 p.

Wilks, S.S., 1962, Mathematical statistics: John Wiley & Sons, Inc., 644 p.

Appendix

This publication includes the source code for custom software to compute the best-fit parameters 

for the National Descriptive Model of Mercury in Fish (NDMMF) in the maximum-likelihood sense. 

The source code is provided as a plain-text (ASCII) file of C code. This C code is based directly on the 

mathematical expressions presented in this publication, with obvious correspondences between key 

variable names in the C code and the terms and symbols used in this report. The software is provided 

subject to the disclaimer contained within the source-code file.
Donato—Computing Maximum-Likelihood Estimates for Parameters of the National Descriptive Model of Mercury in Fish—Open-File Report 2012–1181