Maximum likelihood estimation of multinomial probit factor analysis models for multivariate t-distribution

We propose a model for multinomial probit factor analysis by assuming t-distribution error in probit factor analysis. To obtain maximum likelihood estimation, we use the Monte Carlo expectation maximization algorithm with its M-step greatly simplified under conditional maximization and its E-step made feasible by Monte Carlo simulation. Standard errors are calculated by using Louis’s method. The methodology is illustrated with numerical simulations.     

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.     

Application of the particle-in-cell method for numerical simulation of sheath plasma

We give a survey of papers on the numerical simulation of the sheath plasma using the particle-in-cell method. We study the problem of the behavior of a plasma bounded in the longitudinal direction of an absorbing wall. The model studied contains charged particles (electrons and ions) that move subject to a self-consistent electrostatic field. New pairs of particles are generated in the region of a distributed source. As a numerical model we use the electrostatic “particle-in-cell” method supplemented with the Emmert model for a bulk source and the algorithm of binary Coulomb collisions using a Monte Carlo method. Only electron-ion collisions are taken into account. Translated fromMetody Matematicheskogo Modelirovaniya, 1998, pp. 101–131.     

Implementable Algorithm for Stochastic Optimization Using Sample Average Approximations

We develop an implementable algorithm for stochastic optimization problems involving probability functions. Such problems arise in the design of structural and mechanical systems. The algorithm consists of a nonlinear optimization algorithm applied to sample average approximations and a precision-adjustment rule. The sample average approximations are constructed using Monte Carlo simulations or importance sampling techniques. We prove that the algorithm converges to a solution with probability one and illustrate its use by an example involving a reliability-based optimal design.     

An interacting multiple model algorithm with a switching Markov chain

A Markov chain plays an important role in an interacting multiple model (IMM) algorithm which has been shown to be effective for target tracking systems. Such systems are described by a mixing of continuous states and discrete modes. The switching between system modes is governed by a Markov chain. In real world applications, this Markov chain may change or needs to be changed. Therefore, one may be concerned about a target tracking algorithm with the switching of a Markov chain. This paper concentrates on fault-tolerant algorithm design and algorithm analysis of IMM estimation with the switching of a Markov chain. Monte Carlo simulations are carried out and several conclusions are given.     

Parallel Speed-Up of Monte Carlo Methods for Global Optimization

We introduce the notion of expected hitting time to a goal as a measure of the convergence rate of a Monte Carlo optimization method. The techniques developed apply to simulated annealing, genetic algorithms, and other stochastic search schemes. The expected hitting time can itself be calculated from the more fundamental complementary hitting time distribution (CHTD) which completely characterizes a Monte Carlo method. The CHTD is asymptotically a geometric series, (1/s)/(1 − λ), characterized by two parameters, s, λ, related to the search process in a simple way. The main utility of the CHTD is in comparing Monte Carlo algorithms. In particular we show that independent, identical Monte Carlo algorithms run in parallel, IIP parallelism, and exhibit superlinear speedup. We give conditions under which this occurs and note that equally likely search is linearly sped up. Further we observe that a serial Monte Carlo search can have an infinite expected hitting time, but the same algorithm when parallelized can have a finite expected hitting time. One consequence of the observed superlinear speedup is an improved uniprocessor algorithm by the technique of in-code parallelism.     

Inchworm Monte Carlo Method for Open Quantum Systems

We investigate in this work a recently proposed diagrammatic quantum Monte Carlo method—the inchworm Monte Carlo method—for open quantum systems. We establish its validity rigorously based on resummation of Dyson series. Moreover, we introduce an integro-differential equation formulation for open quantum systems, which illuminates the mathematical structure of the inchworm algorithm. This new formulation leads to an improvement of the inchworm algorithm by introducing classical deterministic time-integration schemes. The numerical method is validated by applications to the spin-boson model. © 2020 Wiley Periodicals, Inc.     

Computational results of an O(n) volume algorithm

Recently an O^∗(n⁴) volume algorithm has been presented for convex bodies by Lovász and Vempala, where n is the number of dimensions of the convex body. Essentially the algorithm is a series of Monte Carlo integrations. In this paper we describe a computer implementation of the volume algorithm, where we improved the computational aspects of the original algorithm by adding variance decreasing modifications: a stratified sampling strategy, double point integration and orthonormalised estimators. Formulas and methodology were developed so that the errors in each phase of the algorithm can be controlled. Some computational results for convex bodies in dimensions ranging from 2 to 10 are presented as well.     

Monte Carlo EM加速算法

EM算法是近年来常用的求后验众数的估计的一种数据增广算法, 但由于求出其E步中积分的显示表达式有时很困难, 甚至不可能, 限制了其应用的广泛性. 而Monte Carlo EM算法很好地解决了这个问题, 将EM算法中E步的积分用Monte Carlo模拟来有效实现, 使其适用性大大增强. 但无论是EM算法, 还是Monte Carlo EM算法, 其收敛速度都是线性的, 被缺损信息的倒数所控制, 当缺损数据的比例很高时, 收敛速度就非常缓慢. 而Newton-Raphson算法在后验众数的附近具有二次收敛速率. 本文提出Monte Carlo EM加速算法, 将Monte Carlo EM算法与Newton-Raphson算法结合, 既使得EM算法中的E步用Monte Carlo模拟得以实现, 又证明了该算法在后验众数附近具有二次收敛速度. 从而使其保留了Monte Carlo EM算法的优点, 并改进了Monte Carlo EM算法的收敛速度. 本文通过数值例子, 将Monte Carlo EM加速算法的结果与EM算法、Monte Carlo EM算法的结果进行比较, 进一步说明了Monte Carlo EM加速算法的优良性.     

A deflated conjugate gradient method for multiple right hand sides and multiple shifts

We consider the task of computing solutions of linear systems that only differ by a shift with the identity matrix as well as linear systems with several different right-hand sides. In the past, Krylov subspace methods have been developed which exploit either the need for solutions to multiple right-hand sides (e.g. deflation type methods and block methods) or multiple shifts (e.g. shifted CG) with some success. In this paper we present a block Krylov subspace method which, based on a block Lanczos process, exploits both features—shifts and multiple right-hand sides—at once. Such situations arise, for example, in lattice quantum chromodynamics (QCD) simulations within the Rational Hybrid Monte Carlo (RHMC) algorithm. We present numerical evidence that our method is superior compared to applying other iterative methods to each of the systems individually as well as, in typical situations, to shifted or block Krylov subspace methods.     

Hamiltonian structure and statistically relevant conserved quantities for the truncated Burgers‐Hopf equation

In dynamical systems with intrinsic chaos, many degrees of freedom, and many conserved quantities, a fundamental issue is the statistical relevance of suitable subsets of these conserved quantities in appropriate regimes. The Galerkin truncation of the Burgers‐Hopf equation has been introduced recently as a prototype model with solutions exhibiting intrinsic stochasticity and a wide range of correlation scaling behavior that can be predicted successfully by simple scaling arguments. Here it is established that the truncated Burgers‐Hopf model is a Hamiltonian system with Hamiltonian given by the integral of the third power. This additional conserved quantity, beyond the energy, has been ignored in previous statistical mechanics studies of this equation. Thus, the question arises of the statistical significance of the Hamiltonian beyond that of the energy. First, an appropriate statistical theory is developed that includes both the energy and Hamiltonian. Then a convergent Monte Carlo algorithm is developed for computing equilibrium statistical distributions. The probability distribution of the Hamiltonian on a microcanonical energy surface is studied through the Monte‐Carlo algorithm and leads to the concept of statistically relevant and irrelevant values for the Hamiltonian. Empirical numerical estimates and simple analysis are combined to demonstrate that the statistically relevant values of the Hamiltonian have vanishingly small measure as the number of degrees of freedom increases with fixed mean energy. The predictions of the theory for relevant and irrelevant values for the Hamiltonian are confirmed through systematic numerical simulations. For statistically relevant values of the Hamiltonian, these simulations show a surprising spectral tilt rather than equipartition of energy. This spectral tilt is predicted and confirmed independently by Monte Carlo simulations based on equilibrium statistical mechanics together with a heuristic formula for the tilt. On the other hand, the theoretically predicted correlation scaling law is satisfied both for statistically relevant and irrelevant values of the Hamiltonian with excellent accuracy. The results established here for the Burgers‐Hopf model are a prototype for similar issues with significant practical importance in much more complex geophysical applications. Several interesting mathematical problems suggested by this study are mentioned in the final section. © 2002 Wiley Periodicals, Inc.     

Discrete Malliavin calculus and computations of greeks in the binomial tree

This paper proposes new methods for computation of greeks using the binomial tree and the discrete Malliavin calculus. In the last decade, the Malliavin calculus has come to be considered as one of the main tools in financial mathematics. It is particularly important in the computation of greeks using Monte Carlo simulations. In previous studies, greeks were usually represented by expectation formulas that are derived from the Malliavin calculus and these expectations are computed using Monte Carlo simulations. On the other hand, the binomial tree approach can also be used to compute these expectations. In this article, we employ the discrete Malliavin calculus to obtain expectation formulas for greeks by the binomial tree method. All the results are obtained in an elementary manner.     

Clustering of contingency table and mixture model

Basing cluster analysis on mixture models has become a classical and powerful approach. It enables some classical criteria such as the well-known k-means criterion to be explained. To classify the rows or the columns of a contingency table, an adapted version of k-means known as Mndki2, which uses the chi-square distance, can be used. Unfortunately, this simple, effective method which can be used jointly with correspondence analysis based on the same representation of the data, cannot be associated with a mixture model in the same way as the classical k-means algorithm. In this paper we show that the Mndki2 algorithm can be viewed as an approximation of a classifying version of the EM algorithm for a mixture of multinomial distributions. A comparison of the algorithms belonging in this context are experimentally investigated using Monte Carlo simulations.     

A quasi-Monte Carlo scheme for Smoluchowski's coagulation equation

This paper analyzes a Monte Carlo algorithm for solving Smoluchowski's coagulation equation. A finite number of particles approximates the initial mass distribution. Time is discretized and random numbers are used to move the particles according to the coagulation dynamics. Convergence is proved when quasi-random numbers are utilized and if the particles are relabeled according to mass in every time step. The results of some numerical experiments show that the error of the new algorithm is smaller than the error of a standard Monte Carlo algorithm using pseudo-random numbers without reordering the particles.

     

Efficiency of Monte Carlo EM and Simulated Maximum Likelihood in Two-Stage Hierarchical Models

Likelihood estimation in hierarchical models is often complicated by the fact that the likelihood function involves an analytically intractable integral. Numerical approximation to this integral is an option but it is generally not recommended when the integral dimension is high. An alternative approach is based on the ideas of Monte Carlo integration, which approximates the intractable integral by an empirical average based on simulations. This article investigates the efficiency of two Monte Carlo estimation methods, the Monte Carlo EM (MCEM) algorithm and simulated maximum likelihood (SML). We derive the asymptotic Monte Carlo errors of both methods and show that, even under the optimal SML importance sampling distribution, the efficiency of SML decreases rapidly (relative to that of MCEM) as the missing information about the unknown parameter increases. We illustrate our results in a simple mixed model example and perform a simulation study which shows that, compared to MCEM, SML can be extremely inefficient in practical applications.     

Monte Carlo evaluation of biological variation: Random generation of correlated non-Gaussian model parameters

When modelling the behaviour of horticultural products, demonstrating large sources of biological variation, we often run into the issue of non-Gaussian distributed model parameters. This work presents an algorithm to reproduce such correlated non-Gaussian model parameters for use with Monte Carlo simulations. The algorithm works around the problem of non-Gaussian distributions by transforming the observed non-Gaussian probability distributions using a proposed SKN-distribution function before applying the covariance decomposition algorithm to generate Gaussian random co-varying parameter sets. The proposed SKN-distribution function is based on the standard Gaussian distribution function and can exhibit different degrees of both skewness and kurtosis. This technique is demonstrated using a case study on modelling the ripening of tomato fruit evaluating the propagation of biological variation with time.     

On the partial synchronization of iterative methods

The rapidly growing field of parallel computing systems promotes the study of parallel algorithms, with the Monte Carlo method and asynchronous iterations being among the most valuable types. These algorithms have a number of advantages. There is no need for a global time in a parallel system (no need for synchronization), and all computational resources are efficiently loaded (the minimum processor idle time). The method of partial synchronization of iterations for systems of equations was proposed by the authors earlier. In this article, this method is generalized to include the case of nonlinear equations of the form x = F(x), where x is an unknown column vector of length n, and F is an operator from ?ⁿ into ?ⁿ. We consider operators that do not satisfy conditions that are sufficient for the convergence of asynchronous iterations, with simple iterations still converging. In this case, one can specify such an incidence of the operator and such properties of the parallel system that asynchronous iterations fail to converge. Partial synchronization is one of the effective ways to solve this problem. An algorithm is proposed that guarantees the convergence of asynchronous iterations and the Monte Carlo method for the above class of operators. The rate of convergence of the algorithm is estimated. The results can prove useful for solving high-dimensional problems on multiprocessor computational systems.     

Automatic selection of indicators in a fully saturated regression

We consider selecting a regression model, using a variant of the general-to-specific algorithm in PcGets, when there are more variables than observations. We look at the special case where the variables are single impulse dummies, one defined for each observation. We show that this setting is unproblematic if tackled appropriately, and obtain the asymptotic distribution of the mean and variance in a location-scale model, under the null that no impulses matter. Monte Carlo simulations confirm the null distributions and suggest extensions to highly non-normal cases. An erratum to this article can be found at