EM Iteration Research Articles

Statistical inference of Poisson censored δ-shock model

In this article, the sufficient statistics, complete statistics, and a uniformly minimum variance unbiased estimate (UMVUE) of the failure parameters of the Poisson censored δ -shock model are obtained, which are based on the sample data of shock arrival time. The maximum likelihood estimates of the shock parameters are also obtained by using the EM algorithm, which is based on the total number of observed shock samples. Finally, the estimation errors and convergence rate of the EM iteration are simulated numerically using Matlab R2016b.

Randomized extrapolation for accelerating EM-type fixed-point algorithms

Several extrapolation strategies have been proposed in the literature to accelerate the EM algorithm, with varying degrees of success. One advantage of extrapolation methods is their ease of implementation, as they only require working with the EM iterations and do not need auxiliary quantities, such as gradient and Hessian. In this paper, we introduce a new family of iterative schemes based on vector extrapolation methods. We also construct and numerically test a randomly relaxed version of the scheme. Our results demonstrate that these new strategies can significantly and stably accelerate the convergence of the EM algorithm compared to existing methods. Moreover, these strategies are highly versatile as they can accelerate any linearly convergent fixed point iteration, including EM-type algorithms. Finally, we provide statistical modeling experiments at the end of the paper to demonstrate the applicability and interest of these convergence acceleration schemes, whether applied to the EM algorithm or one of its variants.

Randomly initialized EM algorithm for two-component Gaussian mixture achieves near optimality in $O(\sqrt{n})$ iterations

We analyze the classical EM algorithm for parameter estimation in the symmetric two-component Gaussian mixtures in d dimensions. We show that, even in the absence of any separation between components, provided that the sample size satisfies n=\Omega(d \log^4 d) , the randomly initialized EM algorithm converges to an estimate in at most O(\sqrt{n}) iterations with high probability, which is at most O((d / n)^{1/4} \log n) in Euclidean distance from the true parameter and within logarithmic factors of the minimax rate of (d / n)^{1/4} . Both the nonparametric statistical rate and the sublinear convergence rate are direct consequences of the zero Fisher information in the worst case. Refined pointwise guarantees beyond worst-case analysis and convergence to the MLE are also shown under mild conditions. This improves the previous result of Balakrishnan, Wainwright, and Yu (2017), which requires strong conditions on both the separation of the components and the quality of the initialization, and that of Daskalakis, Tzamos, and Zampetakis (2017), which requires sample splitting and restarting the EM iteration.

Quadratic extrapolation for accelerating convergence of the EM fixed point problem

Various extrapolation methods have been applied to accelerate convergence of the EM algorithm. These methods are easy to implement, since they work only with EM basic iterations. In other words, auxiliary quantities, such as gradient and hessian, are not needed. In this paper, we define a new family of iterative schemes based on nonlinear extrapolation methods. It is shown that these strategies can accelerate convergence of the EM algorithm much more stably than competing methods. They are extremely general in the sense that they can accelerate any linearly convergent fixed point iterative method, and hence, any EM-type algorithm. A randomly relaxed version is also deduced and numerically tested.

Scattering Noise Model Enhanced EM-TV Algorithm for Benchtop X-ray Fluorescence Computed Tomography Image Reconstruction

Current benchtop x-ray fluorescence computed tomography (XFCT) devices, which use x-ray tubes to stimulate x-ray fluorescence (XRF) photons, suffer from the contamination of Compton scatter background produced by the polychromatic incident beam. The conventional maximum-likelihood expectation-maximization (ML-EM) algorithm only considers the noise model of the XRF signal, which results in high statistical noise in reconstructed images caused by scattered photons. In this study, we proposed a scattering noise model enhanced EM-TV algorithm for benchtop XFCT image reconstruction in order to reduce the noise of scatter background and improve the sensitivity of XFCT images. The statistical noise of scattered photons was considered in the likelihood function and the EM iteration step was modified correspondingly to suppress the statistical noise caused by Compton scattered photons. The robustness of the EM iteration was improved by applying the reweighted total variation (TV) norm as the penalty function. Numerical simulations and imaging experiments of a PMMA phantom consisting of gadolinium (Gd) solutions were performed to validate the proposed algorithm. The phantom was irradiated by a cone-beam polychromatic source and the projection was recorded by a linear-array photon counting detector. For comparison, the XFCT images of Gd were reconstructed using different algorithms. Results indicate that compared with the conventional ML-EM algorithm, the proposed algorithm can obtain XFCT images with lower background noise and higher contrast, which may further improve the sensitivity and image performance of current benchtop XFCT systems.

AR(1) latent class models for longitudinal count data.

In a variety of applications involving longitudinal or repeated-measurements data, it is desired to uncover natural groupings or clusters that exist among study subjects. Motivated by the need to recover clusters of longitudinal trajectories of conduct problems in the field of developmental psychopathology, we propose a method to address this goal when the response data in question are counts. We assume the subject-specific observations are generated from a first-order autoregressive process that is appropriate for count data. A key advantage of our approach is that the class-specific likelihood function arising from each subject's data can be expressed in closed form, circumventing common computational issues associated with random effects models. To further improve computational efficiency, we propose an approximate EM procedure for estimating the model parameters where, within each EM iteration, the maximization step is approximated by solving an appropriately chosen set of estimating equations. We explore the effectiveness of our procedures through simulations based on a four-class model, placing a special emphasis on recovery of the latent trajectories. Finally, we analyze data and recover trajectories of conduct problems in an important nationally representative sample. The methods discussed here are implemented in the R package inarmix, which is available from the Comprehensive R Archive Network (http://cran.r-project.org).

Approximation of Unknown Multivariate Probability Distributions by Using Mixtures of Product Components: A Tutorial

In literature the references to EM estimation of product mixtures are not very frequent. The simplifying assumption of product components, e.g. diagonal covariance matrices in case of Gaussian mixtures, is usually considered only as a compromise because of some computational constraints or limited dataset. We have found that the product mixtures are rarely used intentionally as a preferable approximating tool. Probably, most practitioners do not “trust” the product components because of their formal similarity to “naive Bayes models.” Another reason could be an unrecognized numerical instability of EM algorithm in multidimensional spaces. In this paper we recall that the product mixture model does not imply the assumption of independence of variables. It is even not restrictive if the number of components is large enough. In addition, the product components increase numerical stability of the standard EM algorithm, simplify the EM iterations and have some other important advantages. We discuss and explain the implementation details of EM algorithm and summarize our experience in estimating product mixtures. Finally we illustrate the wide applicability of product mixtures in pattern recognition and in other fields.

Graph theoretic analysis on large DCM models

Graph theoretic analysis on large DCM models

Estimating Parameters of a Probabilistic Heterogeneous Block Model via the EM Algorithm

We introduce a semiparametric block model for graphs, where the within- and between-cluster edge probabilities are not constants within the blocks but are described by logistic type models, reminiscent of the 50-year-old Rasch model and the newly introducedα-βmodels. Our purpose is to give a partition of the vertices of an observed graph so that the induced subgraphs and bipartite graphs obey these models, where their strongly interlaced parameters give multiscale evaluation of the vertices at the same time. In this way, a profoundly heterogeneous version of the stochastic block model is built via mixtures of the above submodels, while the parameters are estimated with a special EM iteration.

A Kernel Smooth Approach for Joint Modeling of Accelerated Failure Time and Longitudinal Data

Joint likelihood approaches have been widely used to handle survival data with time-dependent covariates. In construction of the joint likelihood function for the accelerated failure time (AFT) model, the unspecified baseline hazard function is assumed to be a piecewise constant function in the literature. However, there are usually no close form formulas for the regression parameters, which require numerical methods in the EM iterations. The nonsmooth step function assumption leads to very spiky likelihood function which is very hard to find the globe maximum. Besides, due to nonsmoothness of the likelihood function, direct search methods are conducted for the maximization which are very inefficient and time consuming. To overcome the two disadvantages, we propose a kernel smooth pseudo-likelihood function to replace the nonsmooth step function assumption. The performance of the proposed method is evaluated by simulation studies. A case study of reproductive egg-laying data is provided to demonstrate the usefulness of the new approach.

A hybrid electromagnetism-like algorithm for two-stage assembly flow shop scheduling problem

This paper presents a study on the two-stage assembly flow shop scheduling problem for minimising the weighed sum of maximum makespan, earliness and lateness. There are m machines at the first stage, each of which produces a component of a job. A single machine at the second stage assembles the m components together to complete the job. A novel model for solving the scheduling problem is built to optimise the maximum makespan, earliness and lateness simultaneously. Two optimal operation sequences of jobs are determined and verified. As the problem is known to be NP-hard, a hybrid variable neighbourhood search – electromagnetism-like mechanism (VNS-EM) algorithm is proposed for its handling. To search beyond local optima for a global one, VNS algorithm is embedded in each iteration of EM, whereby the fine neighbourhood search of optimum individuals can be realised and the solution is thus optimised. Simulation results show that the proposed hybrid VNS-EM algorithm outperforms the EM and VNS algorithms in both average value and standard deviation.

The matrix ridge approximation: algorithms and applications

We are concerned with an approximation problem for a symmetric positive semidefinite matrix due to motivation from a class of nonlinear machine learning methods. We discuss an approximation approach that we call matrix ridge approximation. In particular, we define the matrix ridge approximation as an incomplete matrix factorization plus a ridge term. Moreover, we present probabilistic interpretations using a normal latent variable model and a Wishart model for this approximation approach. The idea behind the latent variable model in turn leads us to an efficient EM iterative method for handling the matrix ridge approximation problem. Finally, we illustrate the applications of the approximation approach in multivariate data analysis. Empirical studies in spectral clustering and Gaussian process regression show that the matrix ridge approximation with the EM iteration is potentially useful.

Optimally Sequencing Mixed-model Assembly Lines Optimally with Skip Utility Work Strategy

To improve the work efficiency of the mixed model assembly line, the products sequencing problem with the skip utility work strategy is addressed, where the idle cost and the utility cost are to be optimized simultaneously. Then, a necessary condition of skip utility work and a lower bound of utility work cost are given. The strong NP-hardness of the problem is proved. Since the problem is strongly NP-hard, a hybrid algorithm based on embeded VNS-EM (variable neighborhood search–electromagnetism-like mechanism) algorithm is developed. To escape from the local optima, the enhanced VNS algorithm is embedded in each iteration of EM. With the aid of the good local search ability of VNS algorithm, the fine neighhood search of the optimum individual can be made and the solution is improved. Simulation results confirm the feasibility and validity of this proposed method.

A Fast EM Algorithm for BayesA-Like Prediction of Genomic Breeding Values

Prediction accuracies of estimated breeding values for economically important traits are expected to benefit from genomic information. Single nucleotide polymorphism (SNP) panels used in genomic prediction are increasing in density, but the Markov Chain Monte Carlo (MCMC) estimation of SNP effects can be quite time consuming or slow to converge when a large number of SNPs are fitted simultaneously in a linear mixed model. Here we present an EM algorithm (termed “fastBayesA”) without MCMC. This fastBayesA approach treats the variances of SNP effects as missing data and uses a joint posterior mode of effects compared to the commonly used BayesA which bases predictions on posterior means of effects. In each EM iteration, SNP effects are predicted as a linear combination of best linear unbiased predictions of breeding values from a mixed linear animal model that incorporates a weighted marker-based realized relationship matrix. Method fastBayesA converges after a few iterations to a joint posterior mode of SNP effects under the BayesA model. When applied to simulated quantitative traits with a range of genetic architectures, fastBayesA is shown to predict GEBV as accurately as BayesA but with less computing effort per SNP than BayesA. Method fastBayesA can be used as a computationally efficient substitute for BayesA, especially when an increasing number of markers bring unreasonable computational burden or slow convergence to MCMC approaches.

EM algorithms for nonlinear mixed effects models

Implementing the Monte Carlo EM algorithm (MCEM) algorithm for finding maximum likelihood estimates (MLEs) in the nonlinear mixed effects model (NLMM) has encountered a great deal of difficulty in obtaining samples used for estimating the E step due to the intractability of the target distribution. Sampling methods such as Markov chain techniques and importance sampling have been used to alleviate such difficulty. The advantage of Markov chains is that they are applicable to a wider range of distributions than the approaches based on independent samples. However, in many cases the computational cost of Markov chains is significantly greater than that of independent samplers. The MCEM algorithms based on independent samples allow for straightforward assessment of Monte Carlo error and can be considerably more efficient than those based on Markov chains when an efficient candidate distribution is chosen, which forms the motivation of this paper. The proposed MCEM algorithm in this paper uses samples obtained from an easy-to-simulate and efficient importance distribution so that the computational intensity and complexity is much reduced. Moreover, the proposed MCEM algorithm preserves the flexibility introduced by independent samples in gauging Monte Carlo error and thus allows the Monte Carlo sample size to increase with the number of EM iterations. We also introduce an EM algorithm using Gaussian quadrature approximations (GQEM) for the E step. In low-dimensional cases, the GQEM algorithm is more efficient than the proposed MCEM algorithm and thus can be used as an alternative. The performances of the proposed EM methods are compared to the existing ML estimators using real data examples and simulations.

Incremental Model-Based Clustering for Large Datasets With Small Clusters

Clustering is often useful for analyzing and summarizing information within large datasets. Model-based clustering methods have been found to be effective for determining the number of clusters, dealing with outliers, and selecting the best clustering method in datasets that are small to moderate in size. For large datasets, current model-based clustering methods tend to be limited by memory and time requirements and the increasing difficulty of maximum likelihood estimation. They may fit too many clusters in some portions of the data and/or miss clusters containing relatively few observations. We propose an incremental approach for data that can be processed as a whole in memory, which is relatively efficient computationally and has the ability to find small clusters in large datasets. The method starts by drawing a random sample of the data, selecting and fitting a clustering model to the sample, and extending the model to the full dataset by additional EM iterations. New clusters are then added incrementally, initialized with the observations that are poorly fit by the current model. We demonstrate the effectiveness of this method by applying it to simulated data, and to image data where its performance can be assessed visually.

Initializing EM using the properties of its trajectories in Gaussian mixtures

A strategy is proposed to initialize the EM algorithm in the multivariate Gaussian mixture context. It consists in randomly drawing, with a low computational cost in many situations, initial mixture parameters in an appropriate space including all possible EM trajectories. This space is simply defined by two relations between the two first empirical moments and the mixture parameters satisfied by any EM iteration. An experimental study on simulated and real data sets clearly shows that this strategy outperforms classical methods, since it has the nice property to widely explore local maxima of the likelihood function.

Linear mixed models and penalized least squares

Linear mixed-effects models are an important class of statistical models that are used directly in many fields of applications and also are used as iterative steps in fitting other types of mixed-effects models, such as generalized linear mixed models. The parameters in these models are typically estimated by maximum likelihood or restricted maximum likelihood. In general, there is no closed-form solution for these estimates and they must be determined by iterative algorithms such as EM iterations or general nonlinear optimization. Many of the intermediate calculations for such iterations have been expressed as generalized least squares problems. We show that an alternative representation as a penalized least squares problem has many advantageous computational properties including the ability to evaluate explicitly a profiled log-likelihood or log-restricted likelihood, the gradient and Hessian of this profiled objective, and an ECME update to refine this objective.

On the simulation size and the convergence of the Monte Carlo EM algorithm via likelihood-based distances

When the conditional expectation of a complete-data likelihood in an EM algorithm is analytically intractable, Monte Carlo integration is often used to approximate the E-step. While the resulting Monte Carlo EM algorithm (MCEM) is flexible, assessing convergence of the algorithm is a more difficult task than the original EM algorithm, because of the uncertainty involved in the Monte Carlo approximation. In this note, we propose a convergence criterion using a likelihood-based distance. Because the likelihood is approximated by Monte Carlo integration, we make the distance small with a large probability by selecting the Monte Carlo sample size adaptively at each step of the MCEM algorithm. We implement the proposed convergence criterion along with the simulation size selection in a one-way random effects model. The result shows that our MCEM iterations match the exact EM iterations closely.

Implementations of the Monte Carlo EM Algorithm

The Monte Carlo EM (MCEM) algorithm is a modification of the EM algorithm where the expectation in the E-step is computed numerically through Monte Carlo simulations. The most exible and generally applicable approach to obtaining a Monte Carlo sample in each iteration of an MCEM algorithm is through Markov chain Monte Carlo (MCMC) routines such as the Gibbs and Metropolis–Hastings samplers. Although MCMC estimation presents a tractable solution to problems where the E-step is not available in closed form, two issues arise when implementing this MCEM routine: (1) how do we minimize the computational cost in obtaining an MCMC sample? and (2) how do we choose the Monte Carlo sample size? We address the first question through an application of importance sampling whereby samples drawn during previous EM iterations are recycled rather than running an MCMC sampler each MCEM iteration. The second question is addressed through an application of regenerative simulation. We obtain approximate independent and identical samples by subsampling the generated MCMC sample during different renewal periods. Standard central limit theorems may thus be used to gauge Monte Carlo error. In particular, we apply an automated rule for increasing the Monte Carlo sample size when the Monte Carlo error overwhelms the EM estimate at any given iteration. We illustrate our MCEM algorithm through analyses of two datasets fit by generalized linear mixed models. As a part of these applications, we demonstrate the improvement in computational cost and efficiency of our routine over alternative MCEM strategies.