How to Couple from the Past Using a Read-Once Source of Randomness

We give a new method for generating perfectly random samples from the stationary distribution of a Markov chain. The method is related to coupling from the past (CFTP), but only runs the Markov chain forwards in time, and never restarts it at previous times in the past. The method is also related to an idea known as PASTA (Poisson arrivals see time averages) in the operations research literature. Because the new algorithm can be run using a read-once stream of randomness, we call it read-once CFTP. The memory and time requirements of read-once CFTP are on par with the requirements of the usual form of CFTP, and for a variety of applications the requirements may be noticeably less. Some perfect sampling algorithms for point processes are based on an extension of CFTP known as coupling into and from the past; for completeness, we give a read-once version of coupling into and from the past, but it remains unpractical. For these point process applications, we give an alternative coupling method with which read-once CFTP may be efficiently used.Comment: 28 pages, 2 figure

Highly efficient Bayesian joint inversion for receiver-based data and its application to lithospheric structure beneath the southern Korean Peninsula

With the deployment of extensive seismic arrays, systematic and efficient parameter and uncertainty estimation is of increasing importance and can provide reliable, regional models for crustal and upper-mantle structure.We present an efficient Bayesian method for the joint inversion of surface-wave dispersion and receiver-function data that combines trans-dimensional (trans-D) model selection in an optimization phase with subsequent rigorous parameter uncertainty estimation. Parameter and uncertainty estimation depend strongly on the chosen parametrization such that meaningful regional comparison requires quantitative model selection that can be carried out efficiently at several sites. While significant progress has been made for model selection (e.g. trans-D inference) at individual sites, the lack of efficiency can prohibit application to large data volumes or cause questionable results due to lack of convergence. Studies that address large numbers of data sets have mostly ignored model selection in favour of more efficient/simple estimation techniques (i.e. focusing on uncertainty estimation but employing ad-hoc model choices). Our approach consists of a two-phase inversion that combines trans-D optimization to select the most probable parametrization with subsequent Bayesian sampling for uncertainty estimation given that parametrization. The trans-D optimization is implemented here by replacing the likelihood function with the Bayesian information criterion (BIC). The BIC provides constraints on model complexity that facilitate the search for an optimal parametrization. Parallel tempering (PT) is applied as an optimization algorithm. After optimization, the optimal model choice is identified by the minimum BIC value from all PT chains. Uncertainty estimation is then carried out in fixed dimension. Data errors are estimated as part of the inference problem by a combination of empirical and hierarchical estimation. Data covariance matrices are estimated from data residuals (the difference between prediction and observation) and periodically updated. In addition, a scaling factor for the covariance matrix magnitude is estimated as part of the inversion. The inversion is applied to both simulated and observed data that consist of phase- and group-velocity dispersion curves (Rayleigh wave), and receiver functions. The simulation results show that model complexity and important features are well estimated by the fixed dimensional posterior probability density. Observed data for stations in different tectonic regions of the southern Korean Peninsula are considered. The results are consistent with published results, but important features are better constrained than in previous regularized inversions and are more consistent across the stations. For example, resolution of crustal and Moho interfaces, and absolute values and gradients of velocities in lower crust and upper mantle are better constrained

Bayesian Modelling and Inference on Mixtures of Distributions.

bayesian models;

\v{C}ech-Delaunay gradient flow and homology inference for self-maps

We call a continuous self-map that reveals itself through a discrete set of point-value pairs a sampled dynamical system. Capturing the available information with chain maps on Delaunay complexes, we use persistent homology to quantify the evidence of recurrent behavior. We establish a sampling theorem to recover the eigenspace of the endomorphism on homology induced by the self-map. Using a combinatorial gradient flow arising from the discrete Morse theory for \v{C}ech and Delaunay complexes, we construct a chain map to transform the problem from the natural but expensive \v{C}ech complexes to the computationally efficient Delaunay triangulations. The fast chain map algorithm has applications beyond dynamical systems.Comment: 22 pages, 8 figure

Martingale proofs of many-server heavy-traffic limits for Markovian queues

This is an expository review paper illustrating the ``martingale method'' for proving many-server heavy-traffic stochastic-process limits for queueing models, supporting diffusion-process approximations. Careful treatment is given to an elementary model -- the classical infinite-server model $M/M/\infty$, but models with finitely many servers and customer abandonment are also treated. The Markovian stochastic process representing the number of customers in the system is constructed in terms of rate-1 Poisson processes in two ways: (i) through random time changes and (ii) through random thinnings. Associated martingale representations are obtained for these constructions by applying, respectively: (i) optional stopping theorems where the random time changes are the stopping times and (ii) the integration theorem associated with random thinning of a counting process. Convergence to the diffusion process limit for the appropriate sequence of scaled queueing processes is obtained by applying the continuous mapping theorem. A key FCLT and a key FWLLN in this framework are established both with and without applying martingales.Comment: Published in at http://dx.doi.org/10.1214/06-PS091 the Probability Surveys (http://www.i-journals.org/ps/) by the Institute of Mathematical Statistics (http://www.imstat.org

Spontaneous magnetisation in the plane

The Arak process is a solvable stochastic process which generates coloured patterns in the plane. Patterns are made up of a variable number of random non-intersecting polygons. We show that the distribution of Arak process states is the Gibbs distribution of its states in thermodynamic equilibrium in the grand canonical ensemble. The sequence of Gibbs distributions form a new model parameterised by temperature. We prove that there is a phase transition in this model, for some non-zero temperature. We illustrate this conclusion with simulation results. We measure the critical exponents of this off-lattice model and find they are consistent with those of the Ising model in two dimensions.Comment: 23 pages numbered -1,0...21, 8 figure

SimInf: An R package for Data-driven Stochastic Disease Spread Simulations

We present the R package SimInf which provides an efficient and very flexible framework to conduct data-driven epidemiological modeling in realistic large scale disease spread simulations. The framework integrates infection dynamics in subpopulations as continuous-time Markov chains using the Gillespie stochastic simulation algorithm and incorporates available data such as births, deaths and movements as scheduled events at predefined time-points. Using C code for the numerical solvers and OpenMP to divide work over multiple processors ensures high performance when simulating a sample outcome. One of our design goal was to make SimInf extendable and enable usage of the numerical solvers from other R extension packages in order to facilitate complex epidemiological research. In this paper, we provide a technical description of the framework and demonstrate its use on some basic examples. We also discuss how to specify and extend the framework with user-defined models.Comment: The manual has been updated to the latest version of SimInf (v6.0.0). 41 pages, 16 figure

Harris recurrence of Metropolis-within-Gibbs and trans-dimensional Markov chains

A $\phi$-irreducible and aperiodic Markov chain with stationary probability distribution will converge to its stationary distribution from almost all starting points. The property of Harris recurrence allows us to replace ``almost all'' by ``all,'' which is potentially important when running Markov chain Monte Carlo algorithms. Full-dimensional Metropolis--Hastings algorithms are known to be Harris recurrent. In this paper, we consider conditions under which Metropolis-within-Gibbs and trans-dimensional Markov chains are or are not Harris recurrent. We present a simple but natural two-dimensional counter-example showing how Harris recurrence can fail, and also a variety of positive results which guarantee Harris recurrence. We also present some open problems. We close with a discussion of the practical implications for MCMC algorithms.Comment: Published at http://dx.doi.org/10.1214/105051606000000510 in the Annals of Applied Probability (http://www.imstat.org/aap/) by the Institute of Mathematical Statistics (http://www.imstat.org