PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
J Comput Graph Stat. Author manuscript; available in PMC 2010 September 25.
Published in final edited form as:
J Comput Graph Stat. 2010 September 1; 19(3): 645–665.
doi:  10.1198/jcgs.2010.09014
PMCID: PMC2945396
NIHMSID: NIHMS205488

MM Algorithms for Some Discrete Multivariate Distributions

Abstract

The MM (minorization–maximization) principle is a versatile tool for constructing optimization algorithms. Every EM algorithm is an MM algorithm but not vice versa. This article derives MM algorithms for maximum likelihood estimation with discrete multivariate distributions such as the Dirichlet-multinomial and Connor–Mosimann distributions, the Neerchal–Morel distribution, the negative-multinomial distribution, certain distributions on partitions, and zero-truncated and zero-inflated distributions. These MM algorithms increase the likelihood at each iteration and reliably converge to the maximum from well-chosen initial values. Because they involve no matrix inversion, the algorithms are especially pertinent to high-dimensional problems. To illustrate the performance of the MM algorithms, we compare them to Newton’s method on data used to classify handwritten digits.

Keywords: Dirichlet and multinomial distributions, Inequalities, Maximum likelihood, Minorization

1. INTRODUCTION

The MM algorithm generalizes the celebrated EM algorithm (Dempster, Laird, and Rubin 1977). In this article we apply the MM (minorization–maximization) principle to devise new algorithms for maximum likelihood estimation with several discrete multivariate distributions. A series of research papers and review articles (Groenen 1993; de Leeuw 1994; Heiser 1995; Hunter and Lange 2004; Lange 2004; Wu and Lange 2010) have argued that the MM principle can lead to simpler derivations of known EM algorithms. More importantly, the MM principle also generates many new algorithms of considerable utility. Some statisticians encountering the MM principle for the first time react against its abstraction, unfamiliarity, and dependence on the mathematical theory of inequalities. This is unfortunate because real progress can be made applying a few basic ideas in a unified framework. The current article relies on just three well-known inequalities. For most of our examples, the derivation of a corresponding EM algorithm appears much harder, the main hindrance being the difficulty of choosing an appropriate missing data structure.

Discrete multivariate distributions are seeing wider use throughout statistics. Modern data mining employs such distributions in image reconstruction, pattern recognition, document clustering, movie rating, network analysis, and random graphs. High-dimension data demand high-dimensional models with ten to hundreds of thousands of parameters. Newton’s method and Fisher scoring are capable of finding the maximum likelihood estimates of these distributions via the parameter updates

θ(n+1)=θ(n)+M(θ(n))1L(θ(n)),

where [nabla]L(θ) is the score function and M(θ) is the observed or the expected information matrix, respectively. Several complications can compromise the performance of these traditional algorithms: (a) the information matrix M(θ) may be expensive to compute, (b) it may fail to be positive definite in Newton’s method, (c) in high dimensions it is expensive to solve the linear system M(θ)x = [nabla]L(n)), and (d) if parameter constraints and parameter bounds intrude, then the update itself requires modification. Although mathematical scientists have devised numerous remedies and safeguards, these all come at a cost of greater implementation complexity. The MM principle offers a versatile weapon for attacking optimization problems of this sort. Although MM algorithms have at best a linear rate of convergence, their updates are often very simple. This can tip the computational balance in their favor. In addition, MM algorithms are typically easy to code, numerically stable, and amenable to acceleration. For the discrete distributions considered here, there is one further simplification often missed in the literature. These distributions involve gamma functions. To avoid the complications of evaluating the gamma function and its derivatives, we fall back on a device suggested by Haldane (1941) that replaces ratios of gamma functions by rising polynomials.

Rather than tire the skeptical reader with more preliminaries, it is perhaps best to move on to our examples without delay. The next section defines the MM principle, discusses our three driving inequalities, and reviews two simple acceleration methods. Section 3 derives MM algorithms for some standard multivariate discrete distributions, namely the Dirichlet-multinomial and Connor–Mosimann distributions, the Neerchal–Morel distribution, the negative-multinomial distribution, certain distributions on partitions, and zero-truncated and zero-inflated distributions. Section 4 describes a numerical experiment comparing the performance of the MM algorithms, accelerated MM algorithms, and Newton’s method on model fitting of handwritten digit data. Our discussion concludes by mentioning directions for further research and by frankly acknowledging the limitations of the MM principle.

2. OVERVIEW OF THE MM ALGORITHM

As we have already emphasized, the MM algorithm is a principle for creating algorithms rather than a single algorithm. There are two versions of the MM principle, one for iterative minimization and another for iterative maximization. Here we deal only with the maximization version. Let f (θ) be the objective function we seek to maximize. An MM algorithm involves minorizing f (θ) by a surrogate function g(θ|θn) anchored at the current iterate θn of a search. Minorization is defined by the two properties

f(θn)=g(θn|θn),
(2.1)

f(θ)g(θ|θn),   θθn.
(2.2)

In other words, the surface θ [mapsto] g(θ|θn) lies below the surface θ [mapsto] f (θ) and is tangent to it at the point θ = θn. Construction of the surrogate function g(θ|θn) constitutes the first M of the MM algorithm.

In the second M of the algorithm, we maximize the surrogate g(θ|θn) rather than f (θ). If θn+1 denotes the maximum point of g(θ|θn), then this action forces the ascent property fn+1) ≥ fn). The straightforward proof

f(θn+1)g(θn+1|θn)g(θn|θn)=f(θn)

reflects definitions (2.1) and (2.2) and the choice of θn+1. The ascent property is the source of the MM algorithm’s numerical stability. Strictly speaking, it depends only on increasing g(θ|θn), not on maximizing g(θ|θn).

The art in devising an MM algorithm revolves around intelligent choices of minorizing functions. This brings us to the first of our three basic minorizations

ln(i=1mαi)i=1mαinj=1mαjnln(j=1mαjnαinαi),
(2.3)

invoking the chord below the graph property of the concave function ln x. Note here that all parameter values are positive and that equality obtains whenever αi=αin for all i. Our second basic minorization

ln(c+α)ln(c+αn)1c+αn(ααn)
(2.4)

restates the supporting hyperplane property of the convex function −ln(c + x). Our final basic minorization

ln(1α)ln(1αn)+αn1αnln(ααn)
(2.5)

is just a rearrangement of the two-point information inequality

αnlnαn+(1αn)ln(1αn)αnlnα+(1αn)ln(1α).

Here α and αn must lie in (0, 1). Any standard text on inequalities, for example, the book by Steele (2004), proves these three inequalities. Because piecemeal minorization works well, our derivations apply the basic minorizations only to strategic parts of the overall objective function, leaving other parts untouched.

The convergence theory of MM algorithms is well known (Lange 2004). Convergence to a stationary point is guaranteed provided five properties of the objective function f (θ) and the MM algorithm map M(θ) hold: (a) f (θ) is coercive on its open domain; (b) f (θ) has only isolated stationary points; (c) M(θ) is continuous; (d) θ* is a fixed point of M(θ) if and only if it is a stationary point of f (θ); (e) f [M(θ*)] ≤ f (θ*), with equality if and only if θ* is a fixed point of M(θ). Most of these conditions are easy to verify for our examples, so the details will be omitted.

A common criticism of EM and MM algorithms is their slow convergence. Fortunately, MM algorithms can be easily accelerated (Jamshidian and Jennrich 1995; Lange 1995a; Jamshidian and Jennrich 1997; Varadhan and Rolland 2008). We will employ two versions of the recent square iterative method (SQUAREM) developed by Varadhan and Roland (2008). These simple vector extrapolation techniques require computation of two MM updates at each iteration. Denote the two updates by Mn) and M * Mn), where M(θ) is the MM algorithm map. These updates in turn define two vectors

u=M(θn)θn,   υ=MM(θn)M(θn)u.

The versions diverge in how they compute the steplength constant s. SqMPE1 (minimal polynomial extrapolation) takes s=utuutυ, while SqRRE1 (reduced rank extrapolation) takes s=utυutυ. Once s is specified, we define the next accelerated iterate by θn+1 = θn − 2su + s2υ. Readers should consult the original article for motivation of SQUAREM. Whenever θn+1 decreases the log-likelihood L(θ), we revert to the MM update θn+1 =M * Mn). Finally, we declare convergence when

|L(θn)L(θn1)||L(θn1)|+1<ε.
(2.6)

In the numerical examples that follow, we use the stringent criterion ε = 10−9. More sophisticated stopping criteria based on the gradient of the objective function and the norm of the parameter increment lead to similar results.

3. APPLICATIONS

3.1 Dirichlet-Multinomial and Connor–Mosimann Distributions

When count data exhibit overdispersion, the Dirichlet-multinomial distribution is often substituted for the multinomial distribution. The multinomial distribution is characterized by a vector p = (p1,…, pd) of cell probabilities and a total number of trials m. In the Dirichlet-multinomial sampling, p is first drawn from a Dirichlet distribution with parameter vector α = (α1,…,αd). Once the cell probabilities are determined, multinomial sampling commences. This leads to the admixture density

h(x|α)=Γ(|α|)Γ(α1)Γ(αd)Δd(mx)i=1dpixi+αi1dp1dpd=(mx)Γ(α1+x1)Γ(αd+xd)Γ(|α|+m)Γ(|α|)Γ(α1)Γ(αd),
(3.1)

where |α|=i=1dαi, Δd is the unit simplex in Rd, and x = (x1,…, xd) is the vector of cell counts. Note that the count total |x|=i=1dxi is fixed at m. Standard calculations show that a random vector X drawn from h(x|α) has the means, variances, and covariances

E(Xi)=mαi|α|,Var(Xi)=mαi|α|(1αi|α|)|α|+m|α|+1,Cov(Xi,Xj)=mαi|α|αj|α||α|+m|α|+1,    ij.

If the fractions αi|α| tend to constants pi as |α| tends to ∞, then these moments collapse to the corresponding moments of the multinomial distribution with proportions p1,…, pd.

One of the most unappealing features of the density function h(x|α) is the occurrence of the gamma function. Fortunately, very early on Haldane (1941) noted the alternative representation

h(x|α)=(mx)j=1dαj(αj+1)(αj+xj1)|α|(|α|+1)(|α|+m1).
(3.2)

The replacement of gamma functions by rising polynomials is a considerable gain in simplicity. Bailey (1957) later suggested the reparameterization

πj=αj|α|,    j=1,,d,    θ=1|α|

in terms of the proportion vector π = (π1,…, πd) and the overdispersion parameter θ. In this setting, the discrete density function becomes

h(x|π,θ)=(mx)j=1dπj(πj+θ)[πj+(xj1)θ](1+θ)[1+(m1)θ].
(3.3)

This version of the density function is used to good effect by Griffiths (1973) in implementing Newton’s method for maximum likelihood estimation with the beta-binomial distribution.

In maximum likelihood estimation, we pass to log-likelihoods. This introduces logarithms and turns factors into sums. To construct an MM algorithm under the parameterization (3.2), we need to minorize terms such as ln(αj + k) and −ln(|α| + k). The basic inequalities (2.3) and (2.4) are directly relevant. Suppose we draw t independent samples x1,…, xt from the Dirichlet-multinomial distribution with mi trials for sample i. The term −ln(|α| + k) occurs in the log-likelihood for xi if and only if mik + 1. Likewise the term ln(αj + k) occurs in the log-likelihood for xi if and only if xijk + 1. It follows that the log-likelihood for the entire sample can be written as

L(α)=krkln(|α|+k)+jksjkln(αj+k),rk=i1{mik+1},     sjk=i1{xijk+1}.
(3.4)

The index k in these formulas ranges from 0 to maxi mi −1.

Applying our two basic minorizations to L(α) yields the surrogate function

g(α|αn)=krk1|αn|+k|α|+jksjkαjnαjn+klnαj

up to an irrelevant additive constant. Equating the partial derivative of the surrogate with respect to αj to 0 produces the simple MM update

αjn+1=(ksjkαjnαjn+k)/(krk|αn|+k).
(3.5)

Minka (2003) derived these updates from a different perspective.

Under the parameterization (3.3), matters are slightly more complicated. Now we minorize the terms −ln(1 + kθ) and ln(πj + kθ) via

log(1+kθ)log(1+kθn)11+kθn(kθkθn)

and

log(πj+kθ)πjnπjn+kθnlog(πjn+kθnπjnπj)+kθnπjn+kθnlog(πjn+kθnkθnkθ).

These minorizations lead to the surrogate function

krkk1+kθnθ+jksjk{πjnπjn+kθnlogπj+kθnπjn+kθnlogθ}

up to an irrelevant constant. Setting the partial derivative with respect to θ equal to 0 yields the MM update

θn+1=(jksjkkθnπjn+kθn)/(krkk1+kθn).
(3.6)

The update of the proportion vector π must be treated as a Lagrange multiplier problem owing to the constraint Σj πj = 1. Familiar arguments produce the MM update

πjn+1=(ksjkπjnπjn+kθn)/(lkslkπlnπln+kθn).
(3.7)

The two updates summarized by (3.5), (3.6), and (3.7) enjoy several desirable properties. First, parameter constraints are built in. Second, stationary points of the log-likelihood are fixed points of the updates. Virtually all MM algorithms share these properties. The update (3.7) also reduces to the maximum likelihood estimate

π^j=ksjklkslk=ixijimi
(3.8)

of the corresponding multinomial proportion when θn = 0.

The estimate (3.8) furnishes a natural initial value πj0. To derive an initial value for the overdispersion parameter θ, consider the first two moments

E(Pj)=αj|α|,    E(Pj2)=αj(αj+1)|α|(|α|+1)

of a Dirichlet distribution with parameter vector α. These identities imply that

j=1dE(Pj2)E(Pj)=|α|+d|α|+1=ρ,

which can be solved for θ = 1/|α| in terms of ρ as θ = (ρ − 1)/(d − ρ). Substituting the estimate

ρ^=ji(xij/mi)2i(xij/mi)

for ρ gives a sensible initial value θ0.

To test our two MM algorithms, we now turn to the beta-binomial data of Haseman and Soares (1976) on male mice exposed to various mutagens. The two outcome categories are (a) dead implants and (b) survived implants. In their first dataset, there are t = 524 observations with between m = 1 and m = 20 trials per observation. Table 1 presents the final log-likelihood, number of iterations, and running time (in seconds) of the two MM algorithms and their SQUAREM accelerations on these data. All MM algorithms converge to the maximum point previously found by the scoring method (Paul, Balasooriya, and Banerjee 2005). For the choice ε = 10−9 in stopping criterion (2.6), the MM algorithm (3.5) takes 700 iterations and 0.1580 sec to converge on a laptop computer. The alternative MM algorithm given in the updates (3.6) and (3.7) takes 339 iterations and 0.1626 sec. Figure 1 depicts the progress of the MM iterates on a contour plot of the log-likelihood. The conventional MM algorithm crawls slowly along the ridge in the contour plot; the accelerated versions SqMPE1 and SqRRE1 significantly reduce both the number of iterations and the running time until convergence.

Figure 1
MM Ascent of the Dirichlet-multinomial log-likelihood surface. A color version of this figure is available in the electronic version of this article.
Table 1
MM algorithms for the Haseman and Soares beta-binomial data.

The Dirichlet-multinomial distribution suffers from two restrictions that limit its applicability, namely the negative correlation of coordinates and the determination of variances by means. It is possible to overcome these restrictions by choosing a more flexible mixing distribution as a prior for the multinomial. Connor and Mosimann (1969) suggested a generalization of the Dirichlet distribution that meets this challenge. The resulting admixed distribution, called the generalized Dirichlet-multinomial distribution, has proved its worth in machine learning problems such as the modeling and clustering of images, handwritten digits, and text documents (Bouguila 2008). It is therefore helpful to derive an MM algorithm for maximum likelihood estimation with this distribution that avoids the complications of gamma/digamma/trigamma functions arising with Newton’s method (Bouguila 2008). The Connor–Mosimann distribution is constructed inductively by the mechanism of stick breaking. Imagine breaking the interval [0, 1] into d subintervals of lengths P1,…, Pd by choosing d − 1 independent beta variates Zi with parameters αi and βi. The length of subinterval 1 is P1 = Z1. Given P1 through Pi, the length of subinterval i + 1 is Pi+1 = Zi+1(1 − P1(...)Pi). The last length Pd = 1 − (P1 + (...) + Pd−1) takes up the slack. Standard calculations show that the Pi have the joint density

g(p|α,β)=j=1d1Γ(αj+βj)Γ(αj)Γ(βj)pjαj1(1i=1jpi)γj,    pΔd,

where γj = βj − αj+1 − βj+1 for j = 1,…, d − 2 and γd−1 = βd−1 − 1. The univariate case (d = 2) corresponds to the beta distribution. The Dirichlet distribution is recovered by taking βj = αj+1 +(...)d. With d − 2 more parameters than the Dirichlet distribution, the Connor–Mosimann distribution is naturally more versatile.

The Connor–Mosimann distribution is again conjugate to the multinomial distribution, and the marginal density of a count vector X over m trials is easily shown to be

Pr(X=x)=Δd(mx)j=1d1pjxjg(p|α,β)dp=(mx)j=1d1Γ(αj+xj)Γ(αj)Γ(βj+yj+1)Γ(βj)Γ(αj+βj)Γ(αj+βj+yj),

where yj=k=jdxk. If we adopt the reparameterization

θj=1αj+βj,    πj=αjαj+βj,    j=1,,d1,

and use the fact that xj + yj+1 = yj, then the density can be re-expressed as

(mx)j=1d1πj[πj+(xj1)θj]×(1πj)[1πj+(yj+11)θj]1[1+(yj1)θj].
(3.9)

Thus, maximum likelihood estimation of the parameter vectors π = (π1,…, πd−1) and θ = (θ1,…, θd−1) by the MM algorithm reduces to the case of d − 1 independent beta-binomial problems.

Let x1,…, xt be a random sample from the generalized Dirichlet-multinomial distribution (3.9) with mi trials for observation xi. Following our reasoning for estimation with the Dirichlet-multinomial, we define the associated counts

rjk=i=1t1{xijk+1},    sjk=i=1t1{yijk+1}

for 1 ≤ jd − 1. In this notation, the reader can readily check that the MM updates become

πjn+1=(krjkπjnπjn+kθjn)/(k[rjkπjnπjn+kθjn+sj+1,k(1πjn)1πjn+kθjn]),θjn+1=(k[rjkkθjnπjn+kθjn+sj+1,kkθjn1πjn+kθjn])/(ksjk1+kθjn).

3.2 Neerchal–Morel Distribution

Neerchal and Morel (1998, 2005) proposed an alternative to the Dirichlet-multinomial distribution that accounts for overdispersion by finite admixture. If x represents count data over m trials and d categories, then their discrete density is

h(x|π,ρ)=j=1dπj(mx)[(1ρ)π1]x1[(1ρ)πj+ρ]xj[(1ρ)πd]xd,
(3.10)

where π = (π1,…, πd) is a vector of proportions and ρ [set membership] [0, 1] is an overdispersion parameter. The Neerchal–Morel distribution collapses to the multinomial distribution when ρ = 0. Straightforward calculations show that the Neerchal–Morel distribution has means, variances, and covariances

E(Xi)=mπi,Var(Xi)=mπi(1πi)[1ρ2+mρ2],Cov(Xi,Xj)=mπiπj[1ρ2+mρ2],    ij.

These are precisely the same as the first- and second-order moments of the Dirichlet-multinomial distribution provided we identify πi = αi/|α| and ρ2 = 1/(|α| + 1).

If we draw t independent samples x1,…, xt from the Neerchal–Morel distribution with mi trials for sample i, then the log-likelihood is

iln{jπj(mixi)[(1ρ)π1]xi1[(1ρ)πj+ρ]xij[(1ρ)πd]xid}.
(3.11)

It is worth bearing in mind that every mixture model yields to the minorization (2.3). This is one of the secrets to the success of the EM algorithm. As a practical matter, explicit minorization via inequality (2.3) is more mechanical and often simpler to implement than performing the E step of the EM algorithm. This is particularly true when several minorizations intervene before we reach the ideal surrogate. Here two successive minorizations are needed.

To state the first minorization, let us abbreviate

Πij=πj[(1ρ)π1]xi1[(1ρ)πj+ρ]xij[(1ρ)πd]xid

and denote by ijn the same quantity evaluated at the nth iterate. In this notation it follows that

ln(jΠij)jwijnln(Πijwijn)=jwijnlnΠijjwijnlnwijn

with weights wijn=ijnlijn. The logarithm splits ln [product]ij into the sum

lnΠij=miln(1ρ)+lnπj+xi1lnπ1++xijln(πj+θ)++xidlnπd

for θ = ρ/(1 − ρ). To separate the parameters πj and θ in the troublesome term ln(πj + θ), we apply the minorization (2.3) again. This produces

ln(πj+θ)πjnπjn+θnln(πjn+θnπjnπj)+θnπjn+θnln(πjn+θnθnθ),

and up to a constant the surrogate function takes the form

ijwijn[kxiklnπk+(1xijθnπjn+θn)lnπj]+ijwijn[(mixijθnπjn+θn)ln(1ρ)+xijθnπjn+θnlnρ].

Standard arguments now yield the updates

πkn+1=(ijwijnxik+jwikn(1xikθnπkn+θn))/(lijwijnxil+liwiln(1xilθnπln+θn)),ρn+1=(ijwijnxijθnπjn+θn)/(imi),θn+1=ρn+11ρn+1.

Table 2 lists convergence results for this MM algorithm and its SQUAREM accelerations on the previously discussed Haseman and Soares data.

Table 2
Performance of the Neerchal–Morel MM algorithms.

3.3 Negative-Multinomial

The motivation for the negative-multinomial distribution comes from multinomial sampling with d + 1 categories assigned probabilities π1,…, πd+1. Sampling continues until category d + 1 accumulates β outcomes. At that moment we count the number of outcomes xi falling in category i for 1 ≤ id. For a given vector x = (x1,…, xd ), elementary combinatorics gives the probability

h(x|β,π)=(β+|x|1|x|)(|x|x)i=1dπixiπd+1β=β(β+1)(β+|x|1)x1!xd!i=1dπixiπd+1β.
(3.12)

This formula continues to make sense even if the positive parameter β is not an integer. For arbitrary β > 0, the most straightforward way to construct the negative-multinomial distribution is to run d independent Poisson processes with intensities π1,…, πd. Wait a gamma distributed length of time with shape parameter β and intensity parameter πd+1. At the expiration of this waiting time, count the number of random events Xi of each type i among the first d categories. The random vector X has precisely the discrete density (3.12).

The Poisson process perspective readily yields the moments

E(Xi)=βπiπd+1,Var(Xi)=βπiπd+1(1+πiπd+1),Cov(Xi,Xj)=βπiπd+1πjπd+1,  ij.
(3.13)

Compared to a Poisson distributed random variable with the same mean, the component Xi is overdispersed. Also in contrast to the multinomial and Dirichlet-multinomial distributions, the counts from a negative-multinomial are positively correlated. Negative-multinomial sampling is therefore appealing in many applications.

Let x1,…, xt be a random sample from the negative-multinomial distribution with mi = |xi|. To maximize the log-likelihood

L(β,π)=krkln(β+k)+j=1dx·jlnπj+tβlnπd+1ijlnxij!,rk=i1{mik+1},    x·j=ixij,

we must deal with the terms ln(β + k). Fortunately, the minorization (2.4) implies

ln(β+k)βnβn+kln(βn+kβnβ)+kβn+kln(βn+kkk),

leading to the surrogate function

g(β,π|βn,πn)=krkβnβn+klnβ+j=1dx·jlnπj+tβlnπd+1

up to an irrelevant constant. In view of the constraint πd+1=1j=1dπj the stationarity conditions for a maximum of the surrogate reduce to

0=1βkrkβnβn+k+tlnπd+1,    0=x·jπjtβπd+1,    1jd.
(3.14)

Unfortunately, it is impossible to solve this system of equations analytically. There are two resolutions to the dilemma. One is block relaxation (de Leeuw 1994) alternating the updates

βn+1=(krkβnβn+k)/(tlnπd+1n)

and

πd+1n+1=tβn+1k=1dx·k+tβn+1,    πjn+1=x·jk=1dx·k+tβn+1,    1jd.

This strategy enjoys the ascent property of all MM algorithms.

The other possibility is to solve the stationarity equations numerically. It is clear that the system of equations (3.14) reduces to the single equation

0=1βkrkβnβn+k+tln(tβk=1dx·k+tβ)

for β. Equivalently, if we let

α=β1,    m¯=1tijxij,    cn=krkβnβn+k,

then we must find a root of the equation f (α) = αcnt ln(αm + 1) = 0. It is clear that f (α) is a strictly convex function with f (0) = 0 and limα→∞ f (α) = ∞. Furthermore, a little reflection shows that f′ (0) = cntm < 0. Thus, there is a single root of f (α) on the interval (0,∞). Owing to the convexity of f (α), Newton’s method will reliably find the root if started to the right of the minimum of f (α) at α = t/cn − 1/m.

To find initial values, we again resort to the method of moments. Based on the moments (3.13), the mean and variance of |X| = Σk Xj are

E(|X|)=β(1πd+1)πd+1,    Var(|X|)=β(1πd+1)πd+12.

These suggest that we take

β0=x¯2s2x¯,    πd+10=x¯s2,    πj0=πd+10β0x·jt,    1jd,

where

x¯=1ti=1tmi,    s2=1t1i=1t(mix¯)2.

When the data are underdispersed (s2 < [x with macron]), our proposed initial values are not meaningful, but a negative-multinomial model is a poor choice anyway.

3.4 Distributions on Partitions

A partition of a positive integer m into k parts is a vector a = (a1,…, am) of non-negative integers such that Σi ai = k and |a| = Σi iai = m. In population genetics, the partition distributions of Ewens (2004) and Pitman (Pitman 1995; Johnson, Kotz, and Balakrishnan 1997) find wide application. We now develop an MM algorithm for Pitman’s distribution, which generalizes Ewens’s distribution. Pitman’s distribution

Pr(A=a|m,α,θ)=m!i=1|a|1(θ+iα)(θ+1)(θ+m1)×j=1m[(1α)(1α+j2)]i!]aj1aj!

involves two parameters 0 ≤ α < 1 and θ > −α. Ewens’s distribution corresponds to the choice α = 0. We will restrict θ to be positive.

To estimate parameters given u independent partitions a1,…, au from Pitman’s distribution, we use the minorizations (2.3) and (2.4) to derive the minorizations

ln(θ+iα)θnθn+iαnlnθ+iαnθn+iαnlnα+c,ln(1α+i)1αn1αn+iln(1α)+c,ln(θ+i)1θn+iθ+c,

where c is a different irrelevant constant in each case. Assuming aj is a partition of the integer mj, it follows that the log-likelihood is minorized by

iriθnθn+iαnlnθ+iriiαnθn+iαnlnα+isi(1αn)1αn+iln(1α)itiθn+iθ+c,

where

ri=j1{|aj|i+1},    si=jki+2ajk,    ti=j1{mii+1}.

Standard arguments now yield the simple updates

αn+1=(iriiαnθn+iαn)/(iriiαnθn+iαn+isi(1αn)1αn+i),θn+1=(iriθnθn+iαn)/(itiθn+i).

If we set α0 = 0, then in all subsequent iterates αn = 0, and we get the MM updates for Ewens’s distribution. Despite the availability of the moments of the parts Ai (Charalambides 2007), it is not clear how to initialize α and θ. Unfortunately, the alternative suggestion of Nobuaki (2001) does not guarantee that the initial values satisfy the constraints α [set membership] [0, 1) and θ > 0.

3.5 Zero-Truncated and Zero-Inflated Data

In this section we briefly indicate how the MM perspective sheds fresh light on EM algorithms for zero-truncated and zero-inflated data. Once again mastery of a handful of inequalities rather than computation of conditional expectations drives the derivations.

In many discrete probability models, only data with positive counts are observed. Counts that are 0 are missing. If f (x|θ) represents the density of the complete data, then the density of a random sample x1,…, xt of zero-truncated data amounts to

h(x|θ)=i=1tf(xi|θ)1f(0|θ).

Inequality (2.5) immediately implies the minorization

lnh(x|θ)i=1t[lnf(xi|θ)+f(0|θn)1f(0|θn)lnf(0|θ)]+c,

where c is an irrelevant constant. In many models, maximization of this surrogate function is straightforward.

For instance, with zero-truncated data from the binomial, Poisson, and negative-binomial distributions, the MM updates reduce to

pn+1=(ixi)/(imi1(1pn)mi),    λn+1=(ixi)/(i11eλn),pn+1=(imi1(pn)mi)/(i[xi+mi1(pn)mi]).

For observation i of the binomial model, there are xi successes out of mi trials with success probability p per trial. λ is the mean in the Poisson model. For observation i of the negative-binomial model, there are xi failures before mi required successes.

More complicated models can be handled in similar fashion. The key insight in each case is to augment every ordinary observation xi > 0 by a total of f (0|θn)/[1 − f (0|θn)] pseudo-observations of 0 at iteration n. With this amendment, the two MM algorithms for the beta-binomial distribution implemented in (3.5), (3.6), and (3.7) remain valid except that the count variables rk and sjk defining the updated parameters at iteration n become

s1k=i1{xi1k+1},    s2k=i[1{xi2k+1}+f(0|πn,θn)1f(0|πn,θn)],    rk=i[1+f(0|πn,θn)1f(0|πn,θn)]1{mik+1},

where

f(0|πn,θn)=π2n(π2n+θn)[π2n+(mi1)θn](1+θn)[1+(mi1)θn].

Here category 1 represents success and category 2 failure. If we start with θ0 = 0, then we recover the updates for the zero-truncated binomial distribution.

Zero-inflated data are equally easy to handle. The density function is now

h(x|θ,π)i=1t[(1π)+πf(0|θ)]1{xi=0}[πf(xi|θ)]1{xi>0}.

Inequality (2.3) entails the minorization

lnh(x|θ,π)i=1t1{xi=0}{znln(1π)+(1zn)[lnπ+lnf(0|θn)]}+i=1t1{xi>0}[lnπ+lnf(xi|θ)],zn=1πn1πn+πnf(0|θn).

The MM update of the inflation-admixture parameter clearly is

πn+1=1ti=1t[1{xi>0}+1{xi=0}(1zn)].

As a typical example, consider estimation with the zero-inflated Poisson (Patil 2007). The mean λ of the Poisson component is updated by

λn+1=ixii[(1zn)1{xi=0}+1{xi>0}].

In other words, every 0 observation is discounted by the amount zn at iteration n. This makes intuitive sense.

4. A NUMERICAL EXPERIMENT

As a numerical experiment, we fit the Dirichlet-multinomial (two parameterizations) and the Neerchal–Morel distributions to the 3823 training digits in the handwritten digit data from the UCI machine learning repository (Asuncion and Newman 2007). Each normalized 32 × 32 bitmap is divided into 64 blocks of size 4 × 4, and the black pixels are counted in each block. This generates a 64-dimensional count vector for each bitmap. Bouguila (2008) successfully fit mixtures of Connor–Mosimann to the training data and used the estimated models to cluster the test data. For illustrative purposes we now fit the Dirichlet-multinomial (two parameterizations) and Neerchal–Morel models. Based on the majorization (2.3), it is straightforward to extend our MM algorithms to fit finite mixture models using any of the previously encountered multivariate discrete distributions.

Table 3 lists the final log-likelihoods, number of iterations, and running times of the different algorithms tested. The MM and accelerated MM algorithms were coded in plain Matlab script language. Newton’s method was implemented using the fmincon function in the Marla Optimization Toolbox under the interior-point option with user-supplied analytical gradient and Hessian. All iterations started from the initial points θ0 = 1 and π0=α0=(164,,164). The stopping criterion for Newton’s method was tuned to achieve precision comparable to the stopping criterion (2.6) for the MM algorithms. Running times in seconds were recorded from a laptop computer.

Table 3
Numerical experiment. Row 1: MM; Row 2: SqMPE1 MM; Row 3: SqRRE1 MM; Row 4: Newton’s method using the (fmincon) function available in the Matlab Optimization Toolbox.

Inspection of Table 3 demonstrates that the MM algorithms outperform Newton’s method and that acceleration is often very beneficial. The cost of evaluating and inverting the observed information matrices of the Neerchal–Morel model significantly slows Newton’s method even in these problems with only 64 parameters. The observed information matrix of the Dirichlet-multinomial distribution possesses a special structure (diagonal plus rank-1 perturbation) that makes matrix inversion far easier. Table 3 does not show the human effort in devising, programming, and debugging the various algorithms. For Newton’s method, derivation and programming took in excess of one day. Formulas for the score and observed information of the Dirichlet-multinomial and Neerchal–Morel distributions are omitted for the sake of brevity. Fisher’s scoring algorithm was not implemented because it is even more cumbersome than Newton’s method (Neerchal and Morel 2005).

This numerical comparison is merely for illustrative purpose. Numerical analysts have developed quasi-Newton algorithms to mend the defects of Newton’s method. The limited-memory BFGS (LBFGS) algorithm (Nocedal and Wright 2006) is especially pertinent to high-dimensional problems. A systematic comparison of the two methods is worth pursuing.

5. DISCUSSION

In designing algorithms for maximum likelihood estimation, Newton’s method and Fisher scoring come immediately to mind. In the last generation, statisticians have added the EM principle. These are good mental reflexes, but the broader MM principle also deserves serious consideration. In many problems, the EM and MM perspectives lead to the same algorithm. In other situations such as image reconstruction in transmission tomography, it is possible to construct different EM and MM algorithms for the same purpose (Lange 2004). One of the most appealing features of the EM perspective is that it provides a statistical interpretation of algorithm intermediates. Although it is a matter of taste and experience whether inequalities or missing data offer an easier path to algorithm development, the fact that there are two routes adds to the possibilities for new algorithms.

One can argue that applications of minorizations (2.3) and (2.5) are just disguised EM algorithms. This objection misses the point in three ways. First, it does not suggest missing data structures explaining the minorization (2.4) and other less well-know minorizations. Second, it fails to weigh the difficulties of invoking simple inequalities versus calculating conditional expectations. When the creation of an appropriate surrogate function requires several minorizations, the corresponding conditional expectations become harder to execute. For example, although the EM principle dictates adding pseudo-observations for zero-truncated data, it is easy to lose sight of this simple interpretation in complicated examples such as the beta-binomial distribution. The genetic segregation analysis example appearing in chapter 2 of the book by Lange (2002) falls into the same category. Third, it fails to acknowledge the conceptual clarity of the MM principle, which shifts focus away from the probability spaces connected with missing data to the simple act of minorization. For instance, when one undertakes maximum a posteriori estimation, should the E step of the EM algorithm take into account the prior?

Some EM and MM algorithms are notoriously slow to converge. As we noted earlier, slow convergence is partially offset by the simplicity of each iteration. There is a growing body of techniques for accelerating MM algorithms (Jamshidian and Jennrich 1995; Lange 1995a; Jamshidian and Jennrich 1997; Varadhan and Rolland 2008). These techniques often lead to a ten-fold or even a hundred-fold reduction in the number of iterations. The various examples appearing in this article are typical in this regard. On problems with boundaries or nondifferentiable objective functions, acceleration may be less helpful.

Our negative-multinomial example highlights two useful tactics for overcoming complications in solving the maximization step of the EM and MM algorithms. It is a mistake to think of the various optimization algorithms in isolation. Often block relaxation (de Leeuw 1994) and Newton’s method can be combined creatively with the MM principle. Systematic application of Newton’s method in solving the maximization step of the MM algorithm is formalized in the MM gradient algorithm (Lange 1995b).

Parameter asymptotic standard errors are a natural byproduct of Newton’s method and scoring. With a modicum of additional effort, the EM and MM algorithms also deliver asymptotic standard errors (Meng and Rubin 1991; Hunter and Lange 2004). Virtually all optimization algorithms are prone to converge to inferior modes. For this reason, we have emphasized finding reasonable initial values. The overlooked article of Ueda and Nakano (1998) suggested an annealing approach to maximization with mixture models. Here the idea is to flatten the likelihood surface and eliminate all but the dominant mode. As the iterations proceed, the flat surface gradually warps into the true bumpy surface. Our recent work (Zhou and Lange 2010) extends this idea to many other EM and MM algorithms. A similar idea, called graduated non-convexity (GNC), appears in computer vision and signal processing literature (Blake and Zisserman 1987). In the absence of a good annealing procedure, one can fall back on starting an optimization algorithm from multiple random points, but this inevitably increases the computational load. The reassurance that a log-likelihood is concave is always welcome.

Readers may want to try their hands at devising their own MM algorithms. For instance, the Dirichlet-negative-multinomial distribution, the bivariate Poisson (Johnson, Kotz, and Balakrishnan 1997), and truncated multivariate discrete distributions yield readily to the techniques described. The performance of the MM algorithm on these problems is similar to that in our fully developed examples. Of course, many objective functions are very complicated, and devising a good MM algorithm is a challenge. The greatest payoffs are apt to be on high-dimensional problems. For simplicity of exposition, we have not tackled any extremely high-dimensional problems, but these certainly exist (Sabatti and Lange 2002; Ayers and Lange 2008; Lange and Wu 2008). In any event, most mathematicians and statisticians keep a few tricks up their sleeves. The MM principle belongs there, waiting for the right problems to come along.

Supplementary Material

Matlab Code

ACKNOWLEDGMENTS

The authors thank the editors and referees for their many valuable comments.

Footnotes

SUPPLEMENTAL MATERIALS

Datasets and Matlab codes: The supplementary material (a single zip package) contains all datasets appearing here and the Matlab codes generating our numerical results and graphs. The readme.txt file describes the contents of each file in the package. (supp_material.zip)

Contributor Information

Hua Zhou, Post-Doctoral Fellow, Department of Human Genetics, University of California, Los Angeles, CA 90095-7088 (ude.alcu@uohzauh).

Kenneth Lange, Professor, Departments of Biomathematics, Human Genetics, and Statistics, University of California, Los Angeles, CA 90095-7088.

REFERENCES

  • Asuncion A, Newman DJ. (UCI) Machine Learning Repository. 2007 available at http://www.ics.uci.edu/~mlearn/Repository.html.
  • Ayers KL, Lange K. Penalized Estimation of Haplotype Frequencies. Bioinformatics. 2008;24:1596–1602. [PubMed]
  • Bailey NTJ. The Mathematical Theory of Epidemics. London: Charles Griffin & Company; 1957.
  • Blake A, Zisserman A. Visual Reconstruction. Cambridge, MA: MIT Press; 1987.
  • Bouguila N. Clustering of Count Data Using Generalized Dirichlet Multinomial Distributions. IEEE Transactions on Knowledge and Data Engineering. 2008;20(4):462–474.
  • Charalambides CA. Distributions of Random Partitions and Their Applications. Methodology and Computing in Applied Probability. 2007;9:163–193.
  • Connor RJ, Mosimann JE. Concepts of Independence for Proportions With a Generalization of the Dirichlet Distribution. Journal of the American Statistical Association. 1969;64:194–206.
  • Dempster AP, Laird NM, Rubin DB. Maximum Likelihood From Incomplete Data via the EM Algorithm. Journal of the Royal Statistical Society, Ser. B. 1977;39(1):1–38. (with discussion)
  • de Leeuw J. Block Relaxation Algorithms in Statistics. In: Bock HH, Lenski W, Richter MM, editors. Information Systems and Data Analysis. Berlin: Springer-Verlag; 1994.
  • Ewens WJ. Mathematical Population Genetics. 2nd ed. New York: Springer-Verlag; 2004.
  • Griffiths DA. Maximum Likelihood Estimation for the Beta-Binomial Distribution and an Application to the Household Distribution of the Total Number of Cases of a Disease. Biometrics. 1973;29:637–648. [PubMed]
  • Groenen PJF. The Majorization Approach to Multidimensional Scaling: Some Problems and Extensions. Leiden: The Netherlands: DSWO Press; 1993.
  • Haldane JBS. The Fitting of Binomial Distributions. Annals of Eugenics. 1941;11:179–181.
  • Haseman JK, Soares ER. The Distribution of Fetal Death in Control Mice and Its Implications on Statistical Tests for Dominant Lethal Effects. Mutation Research/Fundamental and Molecular Mechanisms of Mutagenesis. 1976;41:277–288. [PubMed]
  • Heiser WJ. Convergent Computing by Iterative Majorization: Theory and Applications in Multidimensional Data Analysis. In: Krzanowski WJ, editor. Recent Advances in Descriptive Multivariate Analysis. Oxford: Clarendon Press; 1995. pp. 157–189.
  • Hunter DR, Lange K. A Tutorial on MM Algorithms. The American Statistician. 2004;58:30–37.
  • Jamshidian M, Jennrich RI. Acceleration of the EM Algorithm by Using Quasi-Newton Methods. Journal of the Royal Statistical Society, Ser. B. 1995;59:569–587.
  • Jamshidian M, Jennrich RI. Quasi-Newton Acceleration of the EM Algorithm. Journal of the Royal Statistical Society, Ser. B. 1997;59:569–587.
  • Johnson NL, Kotz S, Balakrishnan N. Discrete Multivariate Distributions. New York: Wiley; 1997.
  • Lange K. A Quasi-Newton Acceleration of the EM Algorithm. Statistica Sinica. 1995a;5:1–18.
  • Lange K. A Gradient Algorithm Locally Equivalent to the EM Algorithm. Journal of the Royal Statistical Society, Ser. B. 1995b;57(2):425–437.
  • Lange K. Mathematical and Statistical Methods for Genetic Analysis. 2nd ed. New York: Springer-Verlag; 2002.
  • Lange K. Optimization. New York: Springer-Verlag; 2004.
  • Lange K, Wu TT. An MM Algorithm for Multicategory Vertex Discriminant Analysis. Journal of Computational and Graphical Statistics. 2008;17:527–544.
  • Meng XL, Rubin DB. Using EM to Obtain Asymptotic Variance–Covariance Matrices: The SEM Algorithm. Journal of the American Statistical Association. 1991;86:899–909.
  • Minka TP. Estimating a Dirichlet Distribution. Technical report, Microsoft. 2003.
  • Neerchal NK, Morel JG. Large Cluster Results for Two Parametric Multinomial Extra Variation Models. Journal of the American Statistical Association. 1998;93(443):1078–1087.
  • Neerchal NK, Morel JG. An Improved Method for the Computation of Maximum Likelihood Estimates for Multinomial Overdispersion Models. Computational Statistics & Data Analysis. 2005;49(1):33–43.
  • Nobuaki H. Applying Pitman’s Sampling Formula to Microdata Disclosure Risk Assessment. Journal of Official Statistics. 2001;17:499–520.
  • Nocedal J, Wright S. Numerical Optimization. New York: Springer; 2006.
  • Patil MK, Shirke DT. Testing Parameter of the Power Series Distribution of a Zero Inflated Power Series Model. Statistical Methodology. 2007;4:393–406.
  • Paul SR, Balasooriya U, Banerjee T. Fisher Information Matrix of the Dirichlet-Multinomial Distribution. Biometrical Journal. 2005;47(2):230–236. [PubMed]
  • Pitman J. Exchangeable and Partially Exchangeable Random Partitions. Probability Theory and Related Fields. 1995;102(2):145–158.
  • Sabatti C, Lange K. Genomewide Motif Identification Using a Dictionary Model. Proceedings of the IEEE. 2002;90:1803–1810.
  • Steele JM. The Cauchy–Schwarz Master Class. MAA Problem Books Series. Washington, DC: Mathematical Association of America; 2004.
  • Ueda N, Nakano R. Deterministic Annealing EM Algorithm. Neural Networks. 1998;11:271–282. [PubMed]
  • Varadhan R, Roland C. Simple and Globally Convergent Methods for Accelerating the Convergence of Any EM Algorithm. Scandinavian Journal of Statistics. 2008;35:335–353.
  • Wu TT, Lange K. The MM Alternative to EM. Statistical Science. 2010 to appear.
  • Zhou H, Lange K. On the Bumpy Road to the Dominant Mode. Scandinavian Journal of Statistics. 2010 to appear. [PMC free article] [PubMed]