The EM algorithm with gradient function update for discrete mixtures with known (fixed) number of components

The paper is focussing on some recent developments in nonparametric mixture distributions. It discusses nonparametric maximum likelihood estimation of the mixing distribution and will emphasize gradient type results, especially in terms of global results and global convergence of algorithms such as vertex direction or vertex exchange method. However, the NPMLE (or the algorithms constructing it) provides also an estimate of the number of components of the mixing distribution which might be not desirable for theoretical reasons or might be not allowed from the physical interpretation of the mixture model. When the number of components is fixed in advance, the before mentioned algorithms can not be used and globally convergent algorithms do not exist up to now. Instead, the EM algorithm is often used to find maximum likelihood estimates. However, in this case multiple maxima are often occuring. An example from a meta-analyis of vitamin A and childhood mortality is used to illustrate the considerable, inferential importance of identifying the correct global likelihood. To improve the behavior of the EM algorithm we suggest a combination of gradient function steps and EM steps to achieve global convergence leading to the EM algorithm with gradient function update (EMGFU). This algorithms retains the number of components to be exactly k and typically converges to the global maximum. The behavior of the algorithm is highlighted at hand of several examples.     

Parameter estimation of a composite production function model based on improved artificial fish swarm algorithm and model application

The traditional Cobb–Douglas production function uses the compact mathematical form to describe the relationship between the production results and production factors in the production technology process. However, in macro-economic production, multi-structured production exists universally. In order to better demonstrate such input–output relation, a composite production function model is proposed in this article. In aspect of model parameter estimation, artificial fish swarm algorithm is applied. The algorithm has satisfactory performance in overcoming local extreme value and acquiring global extreme value. Moreover, realization of the algorithm does not need the gradient value of the objective function. For this reason, it is adaptive to searching space. Through the improved artificial fish swarm algorithm, convergence rate and precision are both considerably improved. In aspect of model application, the composite production function model is mainly used to calculate economic growth factor contribution rate. In this article, a relatively more accurate calculating method is proposed. In the end, empirical analysis on economic growth contribution rate of China is implemented.     

Likelihood estimation of missing cell means in the fixed model analysis of variance

This paper examines the formation of maximum likelihood estimates of cell means in analysis of variance problems for cells with missing observations. Methods of estimating the means for missing cells has a long history which includes iterative maximum likelihood techniques, approximation techniques and ad hoc techniques. The use of the EM algorithm to form maximum likelihood estimates has resolved most of the issues associated with this problem. Implementation of the EM algorithm entails specification of a reduced model. As demonstrated in this paper, when there are several missing cells, it is possible to specify a reduced model that results in an unidentifiable likelihood. The EM algorithm in this case does not converge, although the slow divergence may often be mistaken by the unwary as convergence. This paper presents a simple matrix method of determining whether or not the reduced model results in an identifiable likelihood, and consequently in an EM algorithm that converges. We also show the EM algorithm in this case to be equivalent to a method which yields a closed form solution.     

PARAMETER ESTIMATION FOR LOW-ORDER AUTO-REGRESSIVE MODELS WITH MISSING VALUES

Longitudinal data analysis in epidemiological settings is complicated by large multiplicities of short time series and the occurrence of missing observations. To handle such difficulties Rosner & Muñoz (1988) developed a weighted non-linear least squares algorithm for estimating parameters for first-order autoregressive (AR1) processes with time-varying covariates. This method proved efficient when compared to complete case procedures. Here that work is extended by (1) introducing a different estimation procedure based on the EM algorithm, and (2) formulating estimation techniques for second-order autoregressive models. The second development is important because some of the intended areas of application (adult pulmonary function decline, childhood blood pressure) have autocorrelation functions which decay more slowly than the geometric rate imposed by an AR1 model. Simulation studies are used to compare the three methodologies (non-linear, EM based and complete case) with respect to bias, efficiency and coverage both in the presence and in the absence of time-varying covariates. Differing degrees and mechanisms of missingness are examined. Preliminary results indicate the non-linear approach to be the method of choice: it has high efficiency and is easily implemented. An illustrative example concerning pulmonary function decline in the Netherlands is analyzed using this method.     

Computing the log concave NPMLE for interval censored data

In analyzing interval censored data, a non-parametric estimator is often desired due to difficulties in assessing model fits. Because of this, the non-parametric maximum likelihood estimator (NPMLE) is often the default estimator. However, the estimates for values of interest of the survival function, such as the quantiles, have very large standard errors due to the jagged form of the estimator. By forcing the estimator to be constrained to the class of log concave functions, the estimator is ensured to have a smooth survival estimate which has much better operating characteristics than the unconstrained NPMLE, without needing to specify a parametric family or smoothing parameter. In this paper, we first prove that the likelihood can be maximized under a finite set of parameters under mild conditions, although the log likelihood function is not strictly concave. We then present an efficient algorithm for computing a local maximum of the likelihood function. Using our fast new algorithm, we present evidence from simulated current status data suggesting that the rate of convergence of the log-concave estimator is faster (between \(n^{2/5}\) and \(n^{1/2}\)) than the unconstrained NPMLE (between \(n^{1/3}\) and \(n^{1/2}\)).     

Use of vector divisions in solving quasi-likelihood equations for a Poisson model

Maximum-likelihood estimation technique is known to provide consistent and most efficient regression estimates but often this technique is tedious to implement, particularly in the modelling of correlated count responses. To overcome this limitation, researchers have developed semi- or quasi-likelihood functions that depend only on the correct specification of the mean and variance of the responses rather than on the distribution function. Moreover, quasi-likelihood estimation provides consistent and equally efficient estimates as the maximum-likelihood approach. Basically, the quasi-likelihood estimating function is a non-linear equation constituting of the gradient, Hessian and basic score matrices. Henceforth, to obtain estimates of the regression parameters, the quasi-likelihood equation is solved iteratively using the Newton–Raphson technique. However, the inverse of the Jacobian matrix involved in the Newton–Raphson method may not be easy to compute since the matrix is very close to singularity. In this paper, we consider the use of vector divisions in solving quasi-likelihood equations. The vector divisions are implemented to form secant method formulas. To assess the performance of the use of vector divisions with the secant method, we generate cross-sectional Poisson counts using different sets of mean parameters. We compute the estimates of the regression parameters using the Newton–Raphson technique and vector divisions and compare the number of non-convergent simulations under both algorithms.     

The one-step-late PXEM algorithm

The EM algorithm is a popular method for computing maximum likelihood estimates or posterior modes in models that can be formulated in terms of missing data or latent structure. Although easy implementation and stable convergence help to explain the popularity of the algorithm, its convergence is sometimes notoriously slow. In recent years, however, various adaptations have significantly improved the speed of EM while maintaining its stability and simplicity. One especially successful method for maximum likelihood is known as the parameter expanded EM or PXEM algorithm. Unfortunately, PXEM does not generally have a closed form M-step when computing posterior modes, even when the corresponding EM algorithm is in closed form. In this paper we confront this problem by adapting the one-step-late EM algorithm to PXEM to establish a fast closed form algorithm that improves on the one-step-late EM algorithm by insuring monotone convergence. We use this algorithm to fit a probit regression model and a variety of dynamic linear models, showing computational savings of as much as 99.9%, with the biggest savings occurring when the EM algorithm is the slowest to converge.     

On the $$L_p$$ norms of kernel regression estimators for incomplete data with applications to classification

We consider kernel methods to construct nonparametric estimators of a regression function based on incomplete data. To tackle the presence of incomplete covariates, we employ Horvitz–Thompson-type inverse weighting techniques, where the weights are the selection probabilities. The unknown selection probabilities are themselves estimated using (1) kernel regression, when the functional form of these probabilities are completely unknown, and (2) the least-squares method, when the selection probabilities belong to a known class of candidate functions. To assess the overall performance of the proposed estimators, we establish exponential upper bounds on the \(L_p\) norms, \(1\le p<\infty \), of our estimators; these bounds immediately yield various strong convergence results. We also apply our results to deal with the important problem of statistical classification with partially observed covariates.     

Partitioned algorithms for maximum likelihood and other non-linear estimation

There are a variety of methods in the literature which seek to make iterative estimation algorithms more manageable by breaking the iterations into a greater number of simpler or faster steps. Those algorithms which deal at each step with a proper subset of the parameters are called in this paper partitioned algorithms. Partitioned algorithms in effect replace the original estimation problem with a series of problems of lower dimension. The purpose of the paper is to characterize some of the circumstances under which this process of dimension reduction leads to significant benefits.Four types of partitioned algorithms are distinguished: reduced objective function methods, nested (partial Gauss-Seidel) iterations, zigzag (full Gauss-Seidel) iterations, and leapfrog (non-simultaneous) iterations. Emphasis is given to Newton-type methods using analytic derivatives, but a nested EM algorithm is also given. Nested Newton methods are shown to be equivalent to applying to same Newton method to the reduced objective function, and are applied to separable regression and generalized linear models. Nesting is shown generally to improve the convergence of Newton-type methods, both by improving the quadratic approximation to the log-likelihood and by improving the accuracy with which the observed information matrix can be approximated. Nesting is recommended whenever a subset of parameters is relatively easily estimated. The zigzag method is shown to produce a stable but generally slow iteration; it is fast and recommended when the parameter subsets have approximately uncorrelated estimates. The leapfrog iteration has less guaranteed properties in general, but is similar to nesting and zigzagging when the parameter subsets are orthogonal.     

Likelihood inference for lognormal data with left truncation and right censoring with an illustration

The lognormal distribution is quite commonly used as a lifetime distribution. Data arising from life-testing and reliability studies are often left truncated and right censored. Here, the EM algorithm is used to estimate the parameters of the lognormal model based on left truncated and right censored data. The maximization step of the algorithm is carried out by two alternative methods, with one involving approximation using Taylor series expansion (leading to approximate maximum likelihood estimate) and the other based on the EM gradient algorithm (Lange, 1995). These two methods are compared based on Monte Carlo simulations. The Fisher scoring method for obtaining the maximum likelihood estimates shows a problem of convergence under this setup, except when the truncation percentage is small. The asymptotic variance-covariance matrix of the MLEs is derived by using the missing information principle (Louis, 1982), and then the asymptotic confidence intervals for scale and shape parameters are obtained and compared with corresponding bootstrap confidence intervals. Finally, some numerical examples are given to illustrate all the methods of inference developed here.     

Partially projected gradient algorithms for computing nonparametric maximum likelihood estimates of mixing distributions

There exist primarily three different types of algorithms for computing nonparametric maximum likelihood estimates (NPMLEs) of mixing distributions in the literature, which are the EM-type algorithms, the vertex direction algorithms such as VDM and VEM, and the algorithms based on general constrained optimization techniques such as the projected gradient method. It is known that the projected gradient algorithm may run into stagnation during iterations. When a stagnation occurs, VDM steps need to be added. We argue that the abrupt switch to VDM steps can significantly reduce the efficiency of the projected gradient algorithm, and is usually unnecessary. In this paper, we define a group of partially projected directions, which can be regarded as hybrids of ordinary projected gradient directions and VDM directions. Based on these directions, four new algorithms are proposed for computing NPMLEs of mixing distributions. The properties of the algorithms are discussed and their convergence is proved. Extensive numerical simulations show that the new algorithms outperform the existing methods, especially when a NPMLE has a large number of support points or when high accuracy is required.     

An ADMM with continuation algorithm for non-convex SICA-penalized regression in high dimensions

The smooth integration of counting and absolute deviation (SICA) penalty has been demonstrated theoretically and practically to be effective in non-convex penalization for variable selection. However, solving the non-convex optimization problem associated with the SICA penalty when the number of variables exceeds the sample size remains to be enriched due to the singularity at the origin and the non-convexity of the SICA penalty function. In this paper, we develop an efficient and accurate alternating direction method of multipliers with continuation algorithm for solving the SICA-penalized least squares problem in high dimensions. We establish the convergence property of the proposed algorithm under some mild regularity conditions and study the corresponding Karush–Kuhn–Tucker optimality condition. A high-dimensional Bayesian information criterion is developed to select the optimal tuning parameters. We conduct extensive simulations studies to evaluate the efficiency and accuracy of the proposed algorithm, while its practical usefulness is further illustrated with a high-dimensional microarray study.     

Reparameterizing the generalized linear model to accelerate gibbs sampler convergence

Albert and Chib introduced a complete Bayesian method to analyze data arising from the generalized linear model in which they used the Gibbs sampling algorithm facilitated by latent variables. Recently, Cowles proposed an alternative algorithm to accelerate the convergence of the Albert-Chib algorithm. The novelty in this latter algorithm is achieved by using a Hastings algorithm to generate latent variables and bin boundary parameters jointly instead of individually from their respective full conditionals. In the same spirit, we reparameterize the cumulative-link generalized linear model to accelerate the convergence of Cowles’ algorithm even further. One important advantage of our method is that for the three-bin problem it does not require the Hastings algorithm. In addition, for problems with more than three bins, while the Hastings algorithm is required, we provide a proposal density based on the Dirichlet distribution which is more natural than the truncated normal density used in the competing algorithm. Also, using diagnostic procedures recommended in the literature for the Markov chain Monte Carlo algorithm (both single and multiple runs) we show that our algorithm is substantially better than the one recently obtained. Precisely, our algorithm provides faster convergence and smaller autocorrelations between the iterates. Using the probit link function, extensive results are obtained for the three-bin and the five-bin multinomial ordinal data problems.     

Stochastic proximal-gradient algorithms for penalized mixed models

Motivated by penalized likelihood maximization in complex models, we study optimization problems where neither the function to optimize nor its gradient has an explicit expression, but its gradient can be approximated by a Monte Carlo technique. We propose a new algorithm based on a stochastic approximation of the proximal-gradient (PG) algorithm. This new algorithm, named stochastic approximation PG (SAPG) is the combination of a stochastic gradient descent step which—roughly speaking—computes a smoothed approximation of the gradient along the iterations, and a proximal step. The choice of the step size and of the Monte Carlo batch size for the stochastic gradient descent step in SAPG is discussed. Our convergence results cover the cases of biased and unbiased Monte Carlo approximations. While the convergence analysis of some classical Monte Carlo approximation of the gradient is already addressed in the literature (see Atchadé et al. in J Mach Learn Res 18(10):1–33, 2017), the convergence analysis of SAPG is new. Practical implementation is discussed, and guidelines to tune the algorithm are given. The two algorithms are compared on a linear mixed effect model as a toy example. A more challenging application is proposed on nonlinear mixed effect models in high dimension with a pharmacokinetic data set including genomic covariates. To our best knowledge, our work provides the first convergence result of a numerical method designed to solve penalized maximum likelihood in a nonlinear mixed effect model.

     

A class of absolutely continuous bivariate distributions

Block and Basu bivariate exponential distribution is one of the most popular absolutely continuous bivariate distributions. Extensive work has been done on the Block and Basu bivariate exponential model over the past several decades. Interestingly it is observed that the Block and Basu bivariate exponential model can be extended to the Weibull model also. We call this new model as the Block and Basu bivariate Weibull model. We consider different properties of the Block and Basu bivariate Weibull model. The Block and Basu bivariate Weibull model has four unknown parameters and the maximum likelihood estimators cannot be obtained in closed form. To compute the maximum likelihood estimators directly, one needs to solve a four dimensional optimization problem. We propose to use the EM algorithm for computing the maximum likelihood estimators of the unknown parameters. The proposed EM algorithm can be carried out by solving one non-linear equation at each EM step. Our method can be also used to compute the maximum likelihood estimators for the Block and Basu bivariate exponential model. One data analysis has been preformed for illustrative purpose.     

A computational framework for empirical Bayes inference

In empirical Bayes inference one is typically interested in sampling from the posterior distribution of a parameter with a hyper-parameter set to its maximum likelihood estimate. This is often problematic particularly when the likelihood function of the hyper-parameter is not available in closed form and the posterior distribution is intractable. Previous works have dealt with this problem using a multi-step approach based on the EM algorithm and Markov Chain Monte Carlo (MCMC). We propose a framework based on recent developments in adaptive MCMC, where this problem is addressed more efficiently using a single Monte Carlo run. We discuss the convergence of the algorithm and its connection with the EM algorithm. We apply our algorithm to the Bayesian Lasso of Park and Casella (J. Am. Stat. Assoc. 103:681–686, 2008) and on the empirical Bayes variable selection of George and Foster (J. Am. Stat. Assoc. 87:731–747, 2000).     

A new absolute continuous bivariate generalized exponential distribution

The generalized exponential is the most commonly used distribution for analyzing lifetime data. This distribution has several desirable properties and it can be used quite effectively to analyse several skewed life time data. The main aim of this paper is to introduce absolutely continuous bivariate generalized exponential distribution using the method of Block and Basu (1974). In fact, the Block and Basu exponential distribution will be extended to the generalized exponential distribution. We call the new proposed model as the Block and Basu bivariate generalized exponential distribution, then, discuss its different properties. In this case the joint probability distribution function and the joint cumulative distribution function can be expressed in compact forms. The model has four unknown parameters and the maximum likelihood estimators cannot be obtained in explicit form. To compute the maximum likelihood estimators directly, one needs to solve a four dimensional optimization problem. The EM algorithm has been proposed to compute the maximum likelihood estimations of the unknown parameters. One data analysis is provided for illustrative purposes. Finally, we propose some generalizations of the proposed model and compare their models with each other.     

Estimating a Finite Mixed Exponential Distribution under Progressively Type-II Censored Data

The Type-II progressive censoring scheme has become very popular for analyzing lifetime data in reliability and survival analysis. However, no published papers address parameter estimation under progressive Type-II censoring for the mixed exponential distribution (MED), which is an important model for reliability and survival analysis. This is the problem that we address in this paper. It is noted that maximum likelihood estimation of unknown parameters cannot be obtained in closed form due to the complicated log-likelihood function. We solve this problem by using the EM algorithm. Finally, we obtain closed form estimates of the model. The proposed methods are illustrated by both some simulations and a case analysis.     

Simultaneous Wavelet Deconvolution in Periodic Setting

Abstract.  The paper proposes a method of deconvolution in a periodic setting which combines two important ideas, the fast wavelet and Fourier transform-based estimation procedure of Johnstone et al . [ J. Roy. Statist. Soc. Ser. B 66 (2004) 547] and the multichannel system technique proposed by Casey and Walnut [ SIAM Rev . 36 (1994) 537]. An unknown function is estimated by a wavelet series where the empirical wavelet coefficients are filtered in an adapting non-linear fashion. It is shown theoretically that the estimator achieves optimal convergence rate in a wide range of Besov spaces. The procedure allows to reduce the ill-posedness of the problem especially in the case of non-smooth blurring functions such as boxcar functions: it is proved that additions of extra channels improve convergence rate of the estimator. Theoretical study is supplemented by an extensive set of small-sample simulation experiments demonstrating high-quality performance of the proposed method.     

Ml estimation in the poisson binomial distribution with grouped data via the em algorithm

The maximum likelihood estimation of parameters of the Poisson binomial distribution, based on a sample with exact and grouped observations, is considered by applying the EM algorithm (Dempster et al, 1977). The results of Louis (1982) are used in obtaining the observed information matrix and accelerating the convergence of the EM algorithm substantially. The maximum likelihood estimation from samples consisting entirely of complete (Sprott, 1958) or grouped observations are treated as special cases of the estimation problem mentioned above. A brief account is given for the implementation of the EM algorithm when the sampling distribution is the Neyman Type A since the latter is a limiting form of the Poisson binomial. Numerical examples based on real data are included.