Monte Carlo EM Research Articles

Fast maximum likelihood estimation for general hierarchical models

Hierarchical statistical models are important in applied sciences because they capture complex relationships in data, especially when variables are related by space, time, sampling unit, or other shared features. Existing methods for maximum likelihood estimation that rely on Monte Carlo integration over latent variables, such as Monte Carlo Expectation Maximization (MCEM), suffer from drawbacks in efficiency and/or generality. We harness a connection between sampling-stepping iterations for such methods and stochastic gradient descent methods for non-hierarchical models: many noisier steps can do better than few cleaner steps. We call the resulting methods Hierarchical Model Stochastic Gradient Descent (HMSGD) and show that combining efficient, adaptive step-size algorithms with HMSGD yields efficiency gains. We introduce a one-dimensional sampling-based greedy line search for step-size determination. We implement these methods and conduct numerical experiments for a Gamma-Poisson mixture model, a generalized linear mixed models (GLMMs) with single and crossed random effects, and a multi-species ecological occupancy model with over 3000 latent variables. Our experiments show that the accelerated HMSGD methods provide faster convergence than commonly used methods and are robust to reasonable choices of MCMC sample size.

Hyperspectral Image Denoising by Pixel-Wise Noise Modeling and TV-Oriented Deep Image Prior

Model-based hyperspectral image (HSI) denoising methods have attracted continuous attention in the past decades, due to their effectiveness and interpretability. In this work, we aim at advancing model-based HSI denoising, through sophisticated investigation for both the fidelity and regularization terms, or correspondingly noise and prior, by virtue of several recently developed techniques. Specifically, we formulate a novel unified probabilistic model for the HSI denoising task, within which the noise is assumed as pixel-wise non-independent and identically distributed (non-i.i.d) Gaussian predicted by a pre-trained neural network, and the prior for the HSI image is designed by incorporating the deep image prior (DIP) with total variation (TV) and spatio-spectral TV. To solve the resulted maximum a posteriori (MAP) estimation problem, we design a Monte Carlo Expectation–Maximization (MCEM) algorithm, in which the stochastic gradient Langevin dynamics (SGLD) method is used for computing the E-step, and the alternative direction method of multipliers (ADMM) is adopted for solving the optimization in the M-step. Experiments on both synthetic and real noisy HSI datasets have been conducted to verify the effectiveness of the proposed method.

Multilevel Monte Carlo EM scheme for MV-SDEs with small noise

Multilevel Monte Carlo EM scheme for MV-SDEs with small noise

A Laplace-based model with flexible tail behavior

The proposed multiple scaled contaminated asymmetric Laplace (MSCAL) distribution is an extension of the multivariate asymmetric Laplace distribution to allow for a different excess kurtosis on each dimension and for more flexible shapes of the hyper-contours. These peculiarities are obtained by working on the principal component (PC) space. The structure of the MSCAL distribution has the further advantage of allowing for automatic PC-wise outlier detection – i.e., detection of outliers separately on each PC – when convenient constraints on the parameters are imposed. The MSCAL is fitted using a Monte Carlo expectation-maximization (MCEM) algorithm that uses a Monte Carlo method to estimate the orthogonal matrix of eigenvectors. A simulation study is used to assess the proposed MCEM in terms of computational efficiency and parameter recovery. In a real data application, the MSCAL is fitted to a real data set containing the anthropometric measurements of monozygotic/dizygotic twins. Both a skewed bivariate subset of the full data, perturbed by some outlying points, and the full data are considered.

A Proportional Intensity Model with Frailty for Missing Recurrent Failure Data

In some practical circumstances, data are recorded after the systems have begun operations, and data collection is stopped at a predetermined time or after a predetermined number of failures. In such circumstances, incompleteness of various types exists in the aspect of the missing number of failures and their occurrence times beyond the duration of the pilot study. Additionally, multiple repairable systems may present system-to-system variability caused by differences in the operating environments or working loads of individual systems. With respect to left-truncated and right-censored recurrent failure data from multiple repairable systems, we propose a reliability model based on a proportional intensity model with frailty. The frailty model explicitly models unobserved heterogeneity among systems. Covariates incorporated into the proportional intensity model additionally account for the heterogeneity between different operating conditions. To estimate the model parameters for the left-truncated and right-censored recurrent failure data, a Monte Carlo expectation maximization algorithm is proposed. Details of the estimation of the model parameters and the construction of their confidence intervals are examined. A real-world example and simulation studies under various scenarios show prominent applications of the proportional intensity model with frailty to left-truncated and right-censored multiple repairable systems for reliability prediction.

Statistical inference of multi-state transition model for longitudinal data with measurement error and heterogeneity

. Multi-state transition model is typically used to analyze longitudinal data in medicine and sociology. Moreover, variables in longitudinal studies usually are error-prone, and random effects are heterogeneous, which will result in biased estimates of the interest parameters. This article is intended to estimate the parameters of the multi-state transition model for longitudinal data with measurement error and heterogeneous random effects and further consider the covariate related to the covariance matrix of random effects is also error-prone when the covariate in the transition model is error-prone. We model the covariance matrix of random effects through the modified Cholesky decomposition and propose a pseudo-likelihood method based on the Monte Carlo expectation-maximization algorithm and the Bayesian method based on Markov Chain Monte Carlo to infer and calculate the whole estimates. Meanwhile, we obtain the asymptotic properties and evaluate the finite sample performance of the proposed method by simulation, which is well in terms of Bias, RMSE, and coverage rate of confidence intervals. In addition, we apply the proposed method to the MFUS data.

Misspecification analysis of gamma‐ and inverse Gaussian‐based perturbed degradation processes

AbstractAlbeit not equivalent, in many applications the gamma and the inverse Gaussian processes are treated as if they were. This circumstance makes the misspecification problem of these models interesting and important, especially when data are affected by measurement errors, since noisy/perturbed data do not allow to verify whether the selected model is actually able to adequately fit the real (hidden) degradation process. Motivated by the above considerations, in this paper we conduct a large Monte Carlo study to evaluate whether and how the presence of measurement errors affects this misspecification issue. The study is performed considering as reference models a perturbed gamma process recently proposed in the literature and a new perturbed inverse Gaussian process that share the same non‐Gaussian distributed error term. As an alternative option, we also analyze the more classical case where the error term is Gaussian distributed. We consider both the situation where the true model is the perturbed gamma and the one where it is the perturbed inverse Gaussian. Model parameters are estimated from perturbed data using the maximum likelihood method. Estimates are retrieved by using a new sequential Monte Carlo EM algorithm, which use allows to hugely mitigate the severe numerical issues posed by the direct maximization of the likelihood. The risk of incurring in a misspecification is evaluated as percentage of times the Akaike information criterion leads to select the wrong model. The severity of a misspecification is evaluated in terms of its impact on maximum likelihood estimate of the mean remaining useful life.

A composite Bayesian approach for quantile curve fitting with non-crossing constraints

To fit a set of quantile curves, Bayesian simultaneous quantile curve fitting methods face some challenges in properly specifying a feasible formulation and efficiently accommodating the non-crossing constraints. In this article, we propose a new minimization problem and develop its corresponding Bayesian analysis. The new minimization problem imposes two penalties to control not only the smoothness of fitted quantile curves but also the differences between quantile curves. This enables a direct inference on differences of quantile curves and facilitates improved information sharing among quantiles. After adopting B-spline approximation for the positive smoothing functions in the minimization problem, we specify the pseudo composite asymmetric Laplace likelihood and derive its priors. The computation algorithm, including partially collapsed Gibbs sampling for model parameters and Monte Carlo Expectation-Maximization algorithm for penalty parameters, are provided to carry out the proposed approach. The extensive simulation studies show that, compared with other candidate methods, the proposed approach yields more robust estimation. More advantages of the proposed approach are observed for the extreme quantiles, heavy-tailed random errors, and inference on the differences of quantiles. We also demonstrate the relative performances of the proposed approach and other competing methods through two real data analyses.

Spatio-Temporal Analysis and Prediction of Mass Telecommunication Base Station Failure Events

Large-scale telecommunication systems are lifeline infrastructure in modern society. A telecommunication system typically consists of a huge number of base stations with diverse geographical locations across a country, which highly complicates maintenance operations. To allocate maintenance resource properly, it is important to have a good understanding on the failure pattern of these base stations. Statistical inference of recurrent failures of these base stations is challenging because of the large number of base stations and the spatial correlation of their failure processes. Based on eight-month failure data of telecommunication base stations in Harbin, China, we propose a customized nonhomogeneous Poisson process (NHPP) model for recurrent failure data from telecommunication systems. The model consists of two layers, where the temporal layer applies an NHPP with station-specific frailty for failures of each base station, and the spatial layer uses a multivariate lognormal distribution to characterize the correlation among the frailties. The Monte Carlo EM (MCEM) algorithm is applied to estimate parameters included in the proposed model. We demonstrate the proposed model using the Harbin telecommunication system example with 7725 base stations and 4615 failure records.

Variable selection for joint models of multivariate skew-normal longitudinal and survival data.

Many joint models of multivariate skew-normal longitudinal and survival data have been presented to accommodate for the non-normality of longitudinal outcomes in recent years. But existing work did not consider variable selection. This article investigates simultaneous parameter estimation and variable selection in joint modeling of longitudinal and survival data. The penalized splines technique is used to estimate unknown log baseline hazard function, the rectangle integral method is adopted to approximate conditional survival function. Monte Carlo expectation-maximization algorithm is developed to estimate model parameters. Based on local linear approximations to conditional expectation of likelihood function and penalty function, a one-step sparse estimation procedure is proposed to circumvent the computationally challenge in optimizing the penalized conditional expectation of likelihood function, which is utilized to select significant covariates and trajectory functions, and identify the departure from normality of longitudinal data. The conditional expectation of likelihood function-based Bayesian information criterion is developed to select the optimal tuning parameter. Simulation studies and a real example from the clinical trial are used to illustrate the proposed methodologies.

Modeling and maximum likelihood based inference ofinterval-censored data with unknown upper limits andtime-dependent covariates.

Due to the nature of study design or other reasons, the upper limits of the interval-censored data with multiple visits are unknown. A naïve approach is to treat the last observed time as the exact event time, which may induce biased estimators of the model parameters. In this paper, we first develop a Cox model with time-dependent covariates for the event time and a proportional hazards model with frailty for the gap time. We then construct the upper limits using the latent gap times to resolve the issue of interval-censored event time data with unknown upper limits. A data-augmentation technique and a Monte Carlo EM (MCEM) algorithm are developed to facilitate computation. Theoretical properties of the computational algorithm are also investigated. Additionally, new model comparison criteria are developed to assess the fit of the gap time data as well as the fit of the event time data conditional on the gap time data. Our proposed method compares favorably with competing methods in both simulation study and real data analysis.

A general algorithm for error-in-variables regression modelling using Monte Carlo expectation maximization.

In regression modelling, measurement error models are often needed to correct for uncertainty arising from measurements of covariates/predictor variables. The literature on measurement error (or errors-in-variables) modelling is plentiful, however, general algorithms and software for maximum likelihood estimation of models with measurement error are not as readily available, in a form that they can be used by applied researchers without relatively advanced statistical expertise. In this study, we develop a novel algorithm for measurement error modelling, which could in principle take any regression model fitted by maximum likelihood, or penalised likelihood, and extend it to account for uncertainty in covariates. This is achieved by exploiting an interesting property of the Monte Carlo Expectation-Maximization (MCEM) algorithm, namely that it can be expressed as an iteratively reweighted maximisation of complete data likelihoods (formed by imputing the missing values). Thus we can take any regression model for which we have an algorithm for (penalised) likelihood estimation when covariates are error-free, nest it within our proposed iteratively reweighted MCEM algorithm, and thus account for uncertainty in covariates. The approach is demonstrated on examples involving generalized linear models, point process models, generalized additive models and capture-recapture models. Because the proposed method uses maximum (penalised) likelihood, it inherits advantageous optimality and inferential properties, as illustrated by simulation. We also study the model robustness of some violations in predictor distributional assumptions. Software is provided as the refitME package on R, whose key function behaves like a refit() function, taking a fitted regression model object and re-fitting with a pre-specified amount of measurement error.

A Hamiltonian Monte Carlo EM algorithm for generalized linear mixed models with spatial skew latent variables

A Hamiltonian Monte Carlo EM algorithm for generalized linear mixed models with spatial skew latent variables

Modeling left-truncated degradation data using random drift-diffusion Wiener processes

ABSTRACT For products whose performance characteristic (PC) gradually degrades with time, one usually observes its degradation levels repeatedly to predict its remaining useful life (RUL). Due to the limited storage space of the server and the low resolution of a measurement instrument, we seldom record the low-magnitude degradation values at the early degradation stage in applications. Such observation setting introduces left-truncated degradation data, in which the data collection starts later than the unit’s installation. This brings sampling biases and complicates the degradation data analysis. Moreover, due to the uncontrollable factors in applications, the degradation drift and the degradation diffusion may differ among various units. Motivated by an application of high-speed train bearings, we propose a Wiener process model for the left-truncated degradation data and jointly consider the drift-diffusion random effects. Closed-form formulas are available in the expectation-maximization (EM) algorithm for estimating the model parameters. We derive the RUL distribution in closed form. We also extend the proposed model to the multivariate degradation process. The parameters are estimated with the help of the Monte Carlo EM (MCEM) algorithm. An additional laser application illustrates the performance of the proposed model in RUL prediction, which may help to design a predictive maintenance strategy

Small area estimation of general finite-population parameters based on grouped data

This paper proposes a new model-based approach to small area estimation of general finite-population parameters based on grouped data or frequency data, often available from sample surveys. Grouped data contains information on frequencies of some pre-specified groups in each area, for example, the numbers of households in the income classes. Thus, grouped data provide more detailed insight into small areas than area-level aggregated data. A direct application of the widely used small area methods, such as the Fay–Herriot model for area-level data and nested error regression model for unit-level data, is not appropriate since they are not designed for grouped data. Our novel method adopts the multinomial likelihood function for the grouped data. In order to connect the group probabilities of the multinomial likelihood and the auxiliary variables within the framework of small area estimation, we introduce the unobserved unit-level quantities of interest. They follow a linear mixed model with random intercepts and dispersions after some transformation. Then the probabilities that a unit belongs to the groups can be derived and are used to construct the likelihood function for the grouped data given the random effects. The unknown model parameters (hyperparameters) are estimated by a newly developed Monte Carlo EM algorithm which uses an efficient importance sampling. The empirical best predicts (empirical Bayes estimates) of small area parameters are calculated by a simple Gibbs sampling algorithm. The numerical performance of the proposed method is illustrated based on the model-based and design-based simulations. In the application to the city-level grouped income data of Japan, we complete the patchy maps of the Gini coefficient as well as mean income across the country.

Dynamical Non-compensatory Multidimensional IRT Model Using Variational Approximation.

Multidimensional item response theory (MIRT) is a statistical test theory that precisely estimates multiple latent skills of learners from the responses in a test. Both compensatory and non-compensatory models have been proposed for MIRT: the former assumes that each skill can complement other skills, whereas the latter assumes they cannot. This non-compensatory assumption is convincing in many tests that measure multiple skills; therefore, applying non-compensatory models to such data is crucial for achieving unbiased and accurate estimation. In contrast to tests, latent skills will change over time in daily learning. To monitor the growth of skills, dynamical extensions of MIRT models have been investigated. However, most of them assumed compensatory models, and a model that can reproduce continuous latent states of skills under the non-compensatory assumption has not been proposed thus far. To enable accurate skill tracing under the non-compensatory assumption, we propose a dynamical extension of non-compensatory MIRT models by combining a linear dynamical system and a non-compensatory model. This results in a complicated posterior of skills, which we approximate with a Gaussian distribution by minimizing the Kullback-Leibler divergence between the approximated posterior and the true posterior. The learning algorithm for the model parameters is derived through Monte Carlo expectation maximization. Simulation studies verify that the proposed method is able to reproduce latent skills accurately, whereas the dynamical compensatory model suffers from significant underestimation errors. Furthermore, experiments on an actual data set demonstrate that our dynamical non-compensatory model can infer practical skill tracing and clarify differences in skill tracing between non-compensatory and compensatory models.

Multivariate claim count regression model with varying dispersion and dependence parameters

AbstractThe aim of this paper is to present a regression model for multivariate claim frequency data with dependence structures across the claim count responses, which may be of different sign and range, and overdispersion from the unobserved heterogeneity due to systematic effects in the data. For illustrative purposes, we consider the bivariate Poisson-lognormal regression model with varying dispersion. Maximum likelihood estimation of the model parameters is achieved through a novel Monte Carlo expectation–maximization algorithm, which is shown to have a satisfactory performance when we exemplify our approach to Local Government Property Insurance Fund data from the state of Wisconsin.

Exact Inference with Approximate Computation for Differentially Private Data via Perturbations

This paper discusses how two classes of approximate computation algorithms can be adapted, in a modular fashion, to achieve exact statistical inference from differentially private data products. Considered are approximate Bayesian computation for Bayesian inference, and Monte Carlo Expectation-Maximization for likelihood inference. Up to Monte Carlo error, inference from these algorithms is exact with respect to the joint specification of both the analyst's original data model, and the curator's differential privacy mechanism. Highlighted is a duality between approximate computation on exact data, and exact computation on approximate data, which can be leveraged by a well-designed computational procedure for statistical inference.

A parameter estimation method for multivariate binned Hawkes processes

It is often assumed that events cannot occur simultaneously when modelling data with point processes. This raises a problem as real-world data often contains synchronous observations due to aggregation or rounding, resulting from limitations on recording capabilities and the expense of storing high volumes of precise data. In order to gain a better understanding of the relationships between processes, we consider modelling the aggregated event data using multivariate Hawkes processes, which offer a description of mutually-exciting behaviour and have found wide applications in areas including seismology and finance. Here we generalise existing methodology on parameter estimation of univariate aggregated Hawkes processes to the multivariate case using a Monte Carlo expectation–maximization (MC-EM) algorithm and through a simulation study illustrate that alternative approaches to this problem can be severely biased, with the multivariate MC-EM method outperforming them in terms of MSE in all considered cases.

Estimating three- and four-parameter MIRT models with importance-weighted sampling enhanced variational auto-encoder.

Multidimensional Item Response Theory (MIRT) is widely used in educational and psychological assessment and evaluation. With the increasing size of modern assessment data, many existing estimation methods become computationally demanding and hence they are not scalable to big data, especially for the multidimensional three-parameter and four-parameter logistic models (i.e., M3PL and M4PL). To address this issue, we propose an importance-weighted sampling enhanced Variational Autoencoder (VAE) approach for the estimation of M3PL and M4PL. The key idea is to adopt a variational inference procedure in machine learning literature to approximate the intractable marginal likelihood, and further use importance-weighted samples to boost the trained VAE with a better log-likelihood approximation. Simulation studies are conducted to demonstrate the computational efficiency and scalability of the new algorithm in comparison to the popular alternative algorithms, i.e., Monte Carlo EM and Metropolis-Hastings Robbins-Monro methods. The good performance of the proposed method is also illustrated by a NAEP multistage testing data set.