In this note we show that the Monte Carlo EM algorithm, appropriately constructed with importance re-weighting, monotonically increases a corresponding simulated likelihood. This is result is formally proved but also intuitively explained by a formulation of the problem using auxiliary variables

Papaspiliopoulos, Omiros

Sermaidis, Giorgos

English

Warwick Research Archives Portal Repository

 University of Warwick institutional repository: http://go.warwick.ac.uk/wrap  This paper is made available online in accordance with publisher policies. Please scroll down to view the document itself. Please refer to the repository record for this item and our policy information available from the repository home page for further information.  To see the final version of this paper please visit the publisher’s website. Access to the published version may require a subscription. Author(s):  O Papaspiliopoulos and G Sermaidis Article Title: Monotonicity Properties of the Monte Carlo EM Algorithm and Connections with Simulated Likelihood Year of publication: 2007 Link to published article: http://www2.warwick.ac.uk/fac/sci/statistics/crism/research/2007/paper07-24 Publisher statement:  None   Monotonicity properties of the Monte Carlo EM algorithmand connections with simulated likelihoodBy OMIROS PAPASPILIOPOULOS AND GIORGOS SERMAIDISDepartment of Statistics, Warwick University, Coventry, CV4 7AL, U.K.O.Papaspiliopoulos@warwick.ac.uk, G.Sermaidis@warwick.ac.ukSUMMARYIn this note we show that the Monte Carlo EM algorithm, appropriately constructedwith importance re-weighting, monotonically increases a corresponding simulated likeli-hood. This is result is formally proved but also intuitively explained by a formulation ofthe problem using auxiliary variables.Some keywords: Importance sampling, fixed random seeds, incomplete data, latentstochastic processes1. IntroductionA fairly general description of the incomplete data framework is as follows. Let Qθ be thecommon probability distribution of a pair of processes (Y,X), assumed to be absolutelycontinuous with respect to a dominating measure W with densitydQθ(Y,X) = pi(Y,X | θ) dW(Y,X) .Without loss of generality we will assume that W is a probability measure. Only Y isobserved, with corresponding marginal densitypi(Y | θ) = E[pi(Y,X | θ) | Y ] , (1)where the expectation is taken with respect to dW(X | Y ) and (1) is a density withrespect to the marginal law dW(Y ). Y , X and (Y,X) are referred to as the observed,missing and complete data respectively, and L(θ) := pi(Y | θ) and pi(Y,X | θ) as theobserved and complete likelihood respectively. The aim is to infer θ on the basis of the1CRiSM Paper No. 07-24, www.warwick.ac.uk/go/crismobserved likelihood, the complication being that the latter is typically intractable due tothe expectation involved in (1).A popular deterministic algorithm for finding the maximum likelihood estimatorand the observed information matrix is the Expectation-Maximization (EM) algorithm(Dempster et al., 1977). For a current estimate θi−1, the E-step of the ith iteration ofthe EM algorithm requires the computation of the following function of θ,Q(θ, θi−1) = E[log pi(Y,X | θ) | Y ] , (2)where the expectation is taken with respect to dQθi−1(X | Y ). The M-step replacesθi−1 with some θi, such that Q(θi, θi−1) > Q(θi−1, θi−1). Each iteration of the algorithmincreases the observed likelihood, L(θi) ≥ L(θi−1), thus under regularity conditions (Wu,1983) the θis are guaranteed to converge to local maximizer of the likelihood. The EMalgorithm is a very useful tool when both the E- and M-step can be performed analyt-ically, e.g. when the complete likelihood is in the regular exponential family. However,in many applications, although the complete likelihood is explicit, it is not possible toperform the E-step analytically. A variety of stochastic algorithms which involve simu-lation from dQθ(X | Y ) have been developed for this case (see for example Jank, 2006,for a recent review). We focus on the simulated likelihood (Geyer, 1994) and the MonteCarlo EM algorithm (Wei & Tanner, 1990).Let X1, . . . , XN be a stationary sequence with marginal law dW(X | Y ) and defineL(N)(θ) :=N∑j=1pi(Y,Xj | θ) , (3)Q(N)(θ, θi−1) :=N∑j=1log pi(Y,Xj | θ)pi(Y,Xj | θi−1)/L(N)(θi−1) . (4)Notice that L(N)(θ)/N is an unbiased estimator of (1), thus it forms a simulated like-lihood in the sense of Geyer (1994). Additionally, (4) is an unbiased estimator of (2)thus it defines a Monte Carlo EM algorithm, where samples from dQθi−1(X | Y ) are ob-tained by importance sampling with proposals from dW(X | Y ). The maximizer of (3)under certain conditions (Geyer, 1994; Beskos et al., 2007) converges to the maximum2CRiSM Paper No. 07-24, www.warwick.ac.uk/go/crismlikelihood estimator as N → ∞. The Monte Carlo EM algorithm does not anymoremonotonically increase the likelihood due to the extra randomness introduced by thesimulations in the E-step. As a result it is necessary to increase appropriately N withthe iterations to achieve convergence (Chan & Ledolter, 1995; Fort & Moulines, 2003).Additionally, it is hard to devise stopping rules for the algorithm, because of its randomoscillations once it reaches a high likelihood region. McCulloch (1997) compares the twomethods empirically in the context of generalised linear mixed models and advocates acombination of both as a good strategy.Notice that we have used importance re-weighting to construct the Monte Carlo E-step. This has been considered before, see in particular Quintana et al. (1999); Levine& Casella (2001) for computational reasons when Markov chain Monte Carlo is used toproduce the Xjs, since the same draws can be used for multiple iterations. Our mainresult contained in Theorem 1 states that if the Xjs are fixed throughout the iterationsof the Monte Carlo EM algorithm, then L(N)(θi) ≥ L(N)(θi−1). Thus, this Monte CarloEM monotonically increases the corresponding simulated likelihood. We obtain a similarresult for monotonic increase of the likelihood ratio L(N)(θ)/L(N)(θ0) when we chooseW = Qθ0 for some fixed θ0.2. The main resultTheorem 1. Let X1, . . . , XN be a stationary sequence with marginal law dW(X | Y ),and L(N) and Q(N) defined as in (3) and (4) respectively. We assume that the Xjs arekept fixed throughout the iterations of the Monte Carlo EM algorithm. Then, for any pair(θi−1, θi) such that Q(N)(θi, θi−1) > Q(N)(θi−1, θi−1) it holds that L(N)(θi) > L(N)(θi−1).Proof. Letwj(θ) := pi(Y,Xj | θ) , pij(θ) = wj(θ)/N∑j=1wj(θ) .Then, a direct calculation givesQ(N)(θ, θi−1) =N∑j=1log pij(θ)pij(θi−1) + logL(N)(θ) .3CRiSM Paper No. 07-24, www.warwick.ac.uk/go/crismThen, any pair (θi−1, θi) such that Q(N)(θi, θi−1) > Q(N)(θi−1, θi−1) implies thatN∑j=1logpij(θi)pij(θi−1)pij(θi−1) > logL(N)(θi−1)L(N)(θi),which by Jensen’s inequality proves the result.We have two observations on this main result. Firstly, a similar argument can beemployed when W = Qθ0 for some fixed θ0, to yield that the Monte Carlo EM increasesmonotonically the Monte Carlo likelihood ratio L(N)(θ)/L(N)(θ0). Secondly, there is analternative way to prove and motivate this monotonicity property. One can view thesimulated likelihood L(N)(θ)/N as a marginal of the following “complete” likelihoodpi(Y,X1, . . . , XN , J = j | θ) = 1Npi(Y,Xj | θ) , (5)where this density is with respect to the product of the counting measure on {1, . . . , N},dW(Y ) and (dW(X | Y ))N . A direct calculation verifies that (5) indeed defines aprobability measure (this is particularly easy to check if one assumes that both Qθ andWhave densities with respect to another dominating measure, say the Lebesgue measure).In this formulation, J can be treated as missing data, whereas Y and X1, . . . , XN asobserved data. Then, one can check that up to a constant (4) is the corresponding Qfunction of an ordinary EM algorithm for finding the maximizer of (5).AcknowledgementsThe second author would like to thank the Greek State Scholarship Foundation (IKY)for financial support.ReferencesBeskos, A., Papaspiliopoulos, O. & Roberts, G. O. (2007). Monte Carlo max-imum likelihood estimation for discretely observed diffusion processes. Ann. Statist.To appear.Chan, K. S. & Ledolter, J. (1995). Monte Carlo EM estimation for time series modelsinvolving counts. J. Amer. Statist. Assoc. 90 242–252.4CRiSM Paper No. 07-24, www.warwick.ac.uk/go/crismDempster, A. P., Laird, N. M. & Rubin, D. B. (1977). Maximum likelihood fromincomplete data via the EM algorithm. J. Roy. Statist. Soc. Ser. B 39 1–38. Withdiscussion.Fort, G. & Moulines, E. (2003). Convergence of the Monte Carlo expectation maxi-mization for curved exponential families. Ann. Statist. 31 1220–1259.Geyer, C. J. (1994). On the convergence of Monte Carlo maximum likelihood calcula-tions. J. Roy. Statist. Soc. Ser. B 56 261–274.Jank, W. (2006). The EM algorithm, its stochastic implementation and globaloptimization: Some challenges and opportunities for OR. Available fromhttp://www.smith.umd.edu/faculty/wjank/GA-EM-SaulGass.pdf.Levine, R. A. & Casella, G. (2001). Implementations of the Monte Carlo EM algo-rithm. J. Comput. Graph. Statist. 10 422–439.McCulloch, C. E. (1997). Maximum likelihood algorithms for generalized linear mixedmodels. J. Amer. Statist. Assoc. 92 162–170.Quintana, F., Liu, J. & del Pino, G. (1999). Monte Carlo EM with importancereweighting and its applications in random effects models. Computational Statistics &Data Analysis 29 429–444.Wei, G. C. G. & Tanner, M. A. (1990). A Monte Carlo implementation of the EMalgorithm and the poor man’s data augmentation algorithms. J. Amer. Statist. Assoc.85 699–704.Wu, C. F. J. (1983). On the convergence properties of the EM algorithm. Ann. Statist.11 95–103.5CRiSM Paper No. 07-24, www.warwick.ac.uk/go/crism

A Monte Carlo implementation of the EM algorithm and the poor man's data augmentation algorithms.

Convergence of the Monte Carlo expectation maximization for curved exponential families.

Implementations of the Monte Carlo EM algorithm.

Maximum likelihood algorithms for generalized linear mixed models.

Maximum likelihood from incomplete data via the EM algorithm.

Monte Carlo EM estimation for time series models involving counts.

Monte Carlo EM with importance reweighting and its applications in random e®ects models.

Monte Carlo maximum likelihood estimation for discretely observed di®usion processes.

On the convergence of Monte Carlo maximum likelihood calculations.

On the convergence properties of the EM algorithm.

The EM algorithm, its stochastic implementation and global optimization: Some challenges and opportunities for OR. Available from http://www.smith.umd.edu/faculty/wjank/GA-EM-SaulGass.pdf.

Monotonicity properties of the Monte Carlo EM algorithm and connections with simulated likelihood

http://wrap.warwick.ac.uk/35557/1/WRAP_Papaspilioupoulos_07-24w.pdf

Monotonicity properties of the Monte Carlo EM algorithm and connections with simulated likelihood

Abstract

Similar works

Full text

Available Versions

Warwick Research Archives Portal Repository