Subregion-Adaptive Integration of Functions Having a Dominant Peak

Abstract

Many statistical multiple integration problems involve integrands that have a dominant peak. In applying numerical methods to solve these problems, statisticians have paid relatively little attention to existing quadrature methods and available software developed in the numerical analysis literature. One reason these methods have been largely overlooked, even though they are known to be more efficient than Monte Carlo for well-behaved problems of low dimensionality, may be that when applied naively they are poorly suited for peaked-integrand problems. In this article we use transformations based on “split t” distributions to allow the integrals to be efficiently computed using a subregion-adaptive numerical integration algorithm. Our split t distributions are modifications of those suggested by Geweke and may also be used to define Monte Carlo importance functions. We then compare our approach to Monte Carlo. In the several examples we examine here, we find subregion-adaptive integration to be substantially more efficient than importance sampling.     

Computing steady-state queueing-time distributions of single-server queues:GI X /M/1

This paper presents a computationally efficient method to find the steady-state distributions of actual queueing times of the first customer, as well as of a randomly selected customer, of an arrival group for the queueing systemGI ^X /M/1, and hence the queueing-time distribution of a customer for the systemGI/E _X /1. The distribution of virtual queueing time is also obtained. Approximate analysis based on one or more roots is also discussed. Though the exact detailed as well as approximate computations for a variety of interarrival-time distributions such as generalized Erlang, mixed generalized Erlang, hyperexponential, generalized hyperexponential, and deterministic have been carried out, only representative results in the form of tables have been appended. The results obtained should prove useful to queueing theorists, practitioners, and others.     

Multilevel sequential Monte Carlo: Mean square error bounds under verifiable conditions

In this article, we consider the multilevel sequential Monte Carlo (MLSMC) method of Beskos et al. (Stoch. Proc. Appl. [to appear]). This is a technique designed to approximate expectations w.r.t. probability laws associated to a discretization. For instance, in the context of inverse problems, where one discretizes the solution of a partial differential equation. The MLSMC approach is especially useful when independent, coupled sampling is not possible. Beskos et al. show that for MLSMC the computational effort to achieve a given error, can be less than independent sampling. In this article we significantly weaken the assumptions of Beskos et al., extending the proofs to non-compact state-spaces. The assumptions are based upon multiplicative drift conditions as in Kontoyiannis and Meyn (Electron. J. Probab. 10 [2005]: 61–123). The assumptions are verified for an example.     

A real delivery problem dealt with Monte Carlo Techniques

In this paper we use Monte Carlo Techniques to deal with a real world delivery problem of a food company in Valencia (Spain). The problem is modeled as a set of 11 instances of the well known Vehicle Routing Problem, VRP, with additional time constraints. Given that VRP is a NP-hard problem, a heuristic algorithm, based on Monte Carlo techniques, is implemented. The solution proposed by this heuristic algorithm reaches distance and money savings of about 20% and 5% respectively. This work has been partially supported by thePlan de Incentivo a la Investigación/98 of the Universidad Politécnica de Valencia, under the project “Técnicas Monte Carlo aplicadas a Problemas de Rutas de Vehículos”.     

A direct approach to sojourn times in a busy period of an M/M/1 queue

In this paper we present a direct approach to obtaining joint distributions of various quantities of interest in a busy period in an M/M/1 queue. These quantities are: the sojourn times and waiting times of all the customers in the busy period, the busy period length and the number of customers served in a busy period. Since the evolution of the total workload process between two successive customer arrivals is deterministic, this work gives statistic of the complete evolution of the workload process within a busy period. This work was done when the author was post doctoral fellow with the MAESTRO group at INRIA, Sophia Antipolis, France, and was supported by project no. 2900-IT-1 from the Centre Franco-Indien pour la Promotion de la Recherche Avancee (CEFIPRA).     

Tail asymptotics for M/G/1 type queueing processes with subexponential increments

Bivariate regenerative Markov modulated queueing processes {I _n ,L _n } are described. {I _n } is the phase process, and {L _n } is the level process. Increments in the level process have subexponential distributions. A general boundary behavior at the level 0 is allowed. The asymptotic tail of the cycle maximum, , during a regenerative cycle, , and the asymptotic tail of the stationary random variable L _∞, respectively, of the level process are given and shown to be subexponential with L _∞ having the heavier tail. This revised version was published online in June 2006 with corrections to the Cover Date.     

Sequential Monte Carlo Samplers: Error Bounds and Insensitivity to Initial Conditions

This article addresses finite sample stability properties of sequential Monte Carlo methods for approximating sequences of probability distributions. The results presented herein are applicable in the scenario where the start and end distributions in the sequence are fixed and the number of intermediate steps is a parameter of the algorithm. Under assumptions which hold on noncompact spaces, it is shown that the effect of the initial distribution decays exponentially fast in the number of intermediate steps and the corresponding stochastic error is stable in 𝕃 _p norm.     

Scrambling Sobol' and Niederreiter–Xing Points

Hybrids of equidistribution and Monte Carlo methods of integration can achieve the superior accuracy of the former while allowing the simple error estimation methods of the latter. In particular, randomized (0, m, s)-nets in basebproduce unbiased estimates of the integral, have a variance that tends to zero faster than 1/nfor any square integrable integrand and have a variance that for finitenis never more thane?2.718 times as large as the Monte Carlo variance. Lower bounds thaneare known for special cases. Some very important (t, m, s)-nets havet>0. The widely used Sobol' sequences are of this form, as are some recent and very promising nets due to Niederreiter and Xing. Much less is known about randomized versions of these nets, especially ins>1 dimensions. This paper shows that scrambled (t, m, s)-nets enjoy the same properties as scrambled (0, m, s)-nets, except the sampling variance is guaranteed only to be belowb^t[(b+1)/(b−1)]^stimes the Monte Carlo variance for a least-favorable integrand and finiten.     

Concave-Convex Adaptive Rejection Sampling

We describe a method for generating independent samples from univariate density functions using adaptive rejection sampling without the log-concavity requirement. The method makes use of the fact that many functions can be expressed as a sum of concave and convex functions. Using a concave-convex decomposition, we bound the log-density by separately bounding the concave and convex parts using piecewise linear functions. The upper bound can then be used as the proposal distribution in rejection sampling. We demonstrate the applicability of the concave-convex approach on a number of standard distributions and describe an application to the efficient construction of sequential Monte Carlo proposal distributions for inference over genealogical trees. Computer code for the proposed algorithms is available online.     

Services within a Busy Period of an M/M/1 Queue and Dyck Paths

We analyze the service times of customers in a stable M/M/1 queue in equilibrium depending on their position in a busy period. We give the law of the service of a customer at the beginning, at the end, or in the middle of the busy period. It enables as a by-product to prove that the process of instants of beginning of services is not Poisson. We then proceed to a more precise analysis. We consider a family of polynomial generating series associated with Dyck paths of length 2n and we show that they provide the correlation function of the successive services in a busy period with n+1 customers.     

Dynamically Rescaled Hamiltonian Monte Carlo for Bayesian Hierarchical Models

Dynamically rescaled Hamiltonian Monte Carlo is introduced as a computationally fast and easily implemented method for performing full Bayesian analysis in hierarchical statistical models. The method relies on introducing a modified parameterization so that the reparameterized target distribution has close to constant scaling properties, and thus is easily sampled using standard (Euclidian metric) Hamiltonian Monte Carlo. Provided that the parameterizations of the conditional distributions specifying the hierarchical model are “constant information parameterizations” (CIPs), the relation between the modified- and original parameterization is bijective, explicitly computed, and admit exploitation of sparsity in the numerical linear algebra involved. CIPs for a large catalogue of statistical models are presented, and from the catalogue, it is clear that many CIPs are currently routinely used in statistical computing. A relation between the proposed methodology and a class of explicitly integrated Riemann manifold Hamiltonian Monte Carlo methods is discussed. The methodology is illustrated on several example models, including a model for inflation rates with multiple levels of nonlinearly dependent latent variables. Supplementary materials for this article are available online.     

The Interarrival Hyperexponential Queues:Hk/M/c/N With Balking And Reneging

This paper treats the analytical solution of the truncated interarrival hyperexponential queue. H_k/M/c/N with balking and reneging for general values of k,c and N. The discipline considered here is FIFO. Some previously published results are shown to be special cases of the present results     

Partition-Weighted Monte Carlo Estimation

Although various efficient and sophisticated Markov chain Monte Carlo sampling methods have been developed during the last decade, the sample mean is still a dominant in computing Bayesian posterior quantities. The sample mean is simple, but may not be efficient. The weighted sample mean is a natural generalization of the sample mean. In this paper, a new weighted sample mean is proposed by partitioning the support of posterior distribution, so that the same weight is assigned to observations that belong to the same subset in the partition. A novel application of this new weighted sample mean in computing ratios of normalizing constants and necessary theory are provided. Illustrative examples are given to demonstrate the methodology.     

Honest Importance Sampling With Multiple Markov Chains

Importance sampling is a classical Monte Carlo technique in which a random sample from one probability density, π₁, is used to estimate an expectation with respect to another, π. The importance sampling estimator is strongly consistent and, as long as two simple moment conditions are satisfied, it obeys a central limit theorem (CLT). Moreover, there is a simple consistent estimator for the asymptotic variance in the CLT, which makes for routine computation of standard errors. Importance sampling can also be used in the Markov chain Monte Carlo (MCMC) context. Indeed, if the random sample from π₁ is replaced by a Harris ergodic Markov chain with invariant density π₁, then the resulting estimator remains strongly consistent. There is a price to be paid, however, as the computation of standard errors becomes more complicated. First, the two simple moment conditions that guarantee a CLT in the iid case are not enough in the MCMC context. Second, even when a CLT does hold, the asymptotic variance has a complex form and is difficult to estimate consistently. In this article, we explain how to use regenerative simulation to overcome these problems. Actually, we consider a more general setup, where we assume that Markov chain samples from several probability densities, π₁, …, π_k, are available. We construct multiple-chain importance sampling estimators for which we obtain a CLT based on regeneration. We show that if the Markov chains converge to their respective target distributions at a geometric rate, then under moment conditions similar to those required in the iid case, the MCMC-based importance sampling estimator obeys a CLT. Furthermore, because the CLT is based on a regenerative process, there is a simple consistent estimator of the asymptotic variance. We illustrate the method with two applications in Bayesian sensitivity analysis. The first concerns one-way random effect models under different priors. The second involves Bayesian variable selection in linear regression, and for this application, importance sampling based on multiple chains enables an empirical Bayes approach to variable selection.     

Asymptotics of an Efficient Monte Carlo Estimation for the Transition Density of Diffusion Processes

Discretized simulation is widely used to approximate the transition density of discretely observed diffusions. A recently proposed importance sampler, namely modified Brownian bridge, has gained much attention for its high efficiency relative to other samplers. It is unclear for this sampler, however, how to balance the trade-off between the number of imputed values and the number of Monte Carlo simulations under a given computing resource. This paper provides an asymptotically efficient allocation of computing resource to the importance sampling approach with a modified Brownian bridge as importance sampler. The optimal trade-off is established by investigating two types of errors: Euler discretization error and Monte Carlo error. The main results are illustrated with two simulated examples.      

Performance analysis approximation in a queueing system of type M/G/1

In this work, we apply the strong stability method to obtain an estimate for the proximity of the performance measures in the M/G/1 queueing system to the same performance measures in the M/M/1 system under the assumption that the distributions of the service time are close and the arrival flows coincide. In addition to the proof of the stability fact for the perturbed M/M/1 queueing system, we obtain the inequalities of the stability. These results give with precision the error, on the queue size stationary distribution, due to the approximation. For this, we elaborate from the obtained theoretical results, the STR-STAB algorithm which we execute for a determined queueing system: M/Coxian − 2/1. The accuracy of the approach is evaluated by comparison with simulation results.     

Computing Normalizing Constants for Finite Mixture Models via Incremental Mixture Importance Sampling (IMIS)

This article proposes a method for approximating integrated likelihoods in finite mixture models. We formulate the model in terms of the unobserved group memberships, z, and make them the variables of integration. The integral is then evaluated using importance sampling over the z. We propose an adaptive importance sampling function which is itself a mixture, with two types of component distributions, one concentrated and one diffuse. The more concentrated type of component serves the usual purpose of an importance sampling function, sampling mostly group assignments of high posterior probability. The less concentrated type of component allows for the importance sampling function to explore the space in a controlled way to find other, unvisited assignments with high posterior probability. Components are added adaptively, one at a time, to cover areas of high posterior probability not well covered by the current importance sampling function. The method is called incremental mixture importance sampling (IMIS).

IMIS is easy to implement and to monitor for convergence. It scales easily for higher dimensional mixture distributions when a conjugate prior is specified for the mixture parameters. The simulated values on which the estimate is based are independent, which allows for straightforward estimation of standard errors. The self-monitoring aspects of the method make it easier to adjust tuning parameters in the course of estimation than standard Markov chain Monte Carlo algorithms. With only small modifications to the code, one can use the method for a wide variety of mixture distributions of different dimensions. The method performed well in simulations and in mixture problems in astronomy and medical research.     

On the Stationary Distribution of the GI X /M Y /1 Queueing System

For the GI ^X /M/1 queue, it has been recently proved that there exist geometric distributions that are stochastic lower and upper bounds for the stationary distribution of the embedded Markov chain at arrival epochs. In this note we observe that this is also true for the GI ^X /M ^Y /1 queue. Moreover, we prove that the stationary distribution of its embedded Markov chain is asymptotically geometric. It is noteworthy that the asymptotic geometric parameter is the same as the geometric parameter of the upper bound. This fact justifies previous numerical findings about the quality of the bounds.     

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加速算法的优良性.