Sequential Monte Carlo EM for multivariate probit models

Multivariate probit models (MPM) have the appealing feature of capturing some of the dependence structure between the components of multidimensional binary responses. The key for the dependence modelling is the covariance matrix of an underlying latent multivariate Gaussian. Most approaches to MLE in multivariate probit regression rely on MCEM algorithms to avoid computationally intensive evaluations of multivariate normal orthant probabilities. As an alternative to the much used Gibbs sampler a new SMC sampler for truncated multivariate normals is proposed. The algorithm proceeds in two stages where samples are first drawn from truncated multivariate Student $t$ distributions and then further evolved towards a Gaussian. The sampler is then embedded in a MCEM algorithm. The sequential nature of SMC methods can be exploited to design a fully sequential version of the EM, where the samples are simply updated from one iteration to the next rather than resampled from scratch. Recycling the samples in this manner significantly reduces the computational cost. An alternative view of the standard conditional maximisation step provides the basis for an iterative procedure to fully perform the maximisation needed in the EM algorithm. The identifiability of MPM is also thoroughly discussed. In particular, the likelihood invariance can be embedded in the EM algorithm to ensure that constrained and unconstrained maximisation are equivalent. A simple iterative procedure is then derived for either maximisation which takes effectively no computational time. The method is validated by applying it to the widely analysed Six Cities dataset and on a higher dimensional simulated example. Previous approaches to the Six Cities overly restrict the parameter space but, by considering the correct invariance, the maximum likelihood is quite naturally improved when treating the full unrestricted model.

💡 Research Summary

This paper addresses the computational challenges of maximum‑likelihood estimation in multivariate probit models, where the dependence among binary (or multinomial) responses is captured through the covariance matrix of an underlying latent Gaussian vector. Traditional approaches rely on Monte‑Carlo EM (MCEM) with a Gibbs sampler to draw from truncated multivariate normal distributions in the E‑step. However, Gibbs sampling becomes inefficient in moderate to high dimensions because of slow mixing and strong autocorrelation.

The authors propose a novel Sequential Monte‑Carlo (SMC) sampler specifically designed for truncated multivariate normals. The sampler operates in two stages. First, particles are drawn from a truncated multivariate Student‑t distribution. The heavy‑tailed t‑distribution allows rapid exploration of the constrained region, especially when the truncation cuts off a substantial portion of the Gaussian support. Second, the particles are gradually transformed toward the target truncated Gaussian by increasing the degrees of freedom and shifting the location, while adaptively resampling and applying Metropolis‑Hastings moves to maintain particle diversity. The incremental importance weights are computed using the optimal backward kernel, and the scaling of the Metropolis proposal is tuned online to achieve a desired acceptance rate.

Embedding this SMC sampler within an MCEM framework yields a “sequential EM” algorithm. Instead of generating a fresh MCMC chain at each EM iteration, the particle cloud from the previous iteration is simply re‑weighted and propagated after the M‑step update. This recycling dramatically reduces the computational burden of the E‑step while preserving the unbiasedness of the Monte‑Carlo approximation.

The M‑step itself is tackled by extending the conditional maximisation (CM) strategy of Meng and Rubin (1993). The expected complete‑data log‑likelihood can be expressed as
( Q(\psi,\psi^{(m)}) = -N\log|\Sigma| - N\operatorname{tr}(\Sigma^{-1}S) ),
where (S) is the Monte‑Carlo estimate of the conditional second‑moment matrix of the latent variables. The authors show that, for the multivariate probit, the maximisation over (\Sigma) has a closed‑form solution (\Sigma^{(m+1)} = S/N). The regression coefficients (\beta) are updated by solving a weighted least‑squares problem: (\beta^{(m+1)} = (\sum X_j^\top \Sigma^{-1} X_j)^{-1}\sum X_j^\top \Sigma^{-1} \mathbb{E}