indii.org / research
http://www.indii.org/
Lawrence Murray: software, research, photography.en-usTue, 02 Jan 2018 21:15:22 +0100Tue, 02 Jan 2018 21:15:22 +0100Delayed Sampling
http://www.indii.org/research/delayed-sampling-and-automatic-rao-blackwellization-of-probabilistic-programs/index.html
Fri, 22 Dec 2017 00:00:00 +0100Lawrence Murrayhttp://www.indii.org/research/delayed-sampling-and-automatic-rao-blackwellization-of-probabilistic-programs/index.htmlWe introduce a dynamic mechanism for the solution of analytically-tractable substructure in probabilistic programs, to reduce variance in Monte Carlo estimators. For inference with Sequential Monte Carlo, it yields improvements such as locally-optimal proposals and Rao-Blackwellization, with little modification to model code necessary. A directed graph is maintained alongside the running program, evolving dynamically as the program triggers operations upon it. Nodes of the graph represent random variables, and edges the analytically-tractable relationships between them (e.g. conjugate priors and affine transformations). Each random variable is held in the graph for as long as possible, sampled only when used by the program in a context that cannot be resolved analytically. This allows it to be analytically conditioned on as many observations as possible before sampling. We have implemented the approach in both Anglican and a new probabilistic programming language named Birch. We demonstrate it on a number of small examples, and a larger mixed linear-nonlinear state-space model.
Better together?
http://www.indii.org/research/better-together-statistical-learning-in-models-made-of-modules/index.html
Tue, 29 Aug 2017 00:00:00 +0200Lawrence Murrayhttp://www.indii.org/research/better-together-statistical-learning-in-models-made-of-modules/index.htmlIn modern applications, statisticians are faced with integrating heterogeneous data modalities relevant for an inference, prediction, or decision problem. In such circumstances, it is convenient to use a graphical model to represent the statistical dependencies, via a set of connected “modules”, each relating to a specific data modality, and drawing on specific domain expertise in their development. In principle, given data, the conventional statistical update then allows for coherent uncertainty quantification and information propagation through and across the modules. However, misspecification of any module can contaminate the estimate and update of others, often in unpredictable ways. In various settings, particularly when certain modules are trusted more than others, practitioners have preferred to avoid learning with the full model in favor of approaches that restrict the information propagation between modules, for example by restricting propagation to only particular directions along the edges of the graph. In this article, we investigate why these modular approaches might be preferable to the full model in misspecified settings. We propose principled criteria to choose between modular and full-model approaches. The question arises in many applied settings, including large stochastic dynamical systems, meta-analysis, epidemiological models, air pollution models, pharmacokinetics-pharmacodynamics, and causal inference with propensity scores.
Probabilistic learning of nonlinear dynamical systems using sequential Monte Carlo
http://www.indii.org/research/probabilistic-learning-of-nonlinear-dynamical-systems-using-sequential-monte-carlo/index.html
Tue, 07 Mar 2017 00:00:00 +0100Lawrence Murrayhttp://www.indii.org/research/probabilistic-learning-of-nonlinear-dynamical-systems-using-sequential-monte-carlo/index.htmlProbabilistic modeling provides the capability to represent and manipulate uncertainty in data, models, predictions and decisions. We are concerned with the problem of learning probabilistic models of dynamical systems from measured data. Specifically, we consider learning of probabilistic nonlinear state-space models. There is no closed-form solution available for this problem, implying that we are forced to use approximations. In this tutorial we will provide a self-contained introduction to one of the state-of-the-art methods—the particle Metropolis–Hastings algorithm—which has proven to offer a practical approximation. This is a Monte Carlo based method, where the particle filter is used to guide a Markov chain Monte Carlo method through the parameter space. One of the key merits of the particle Metropolis–Hastings algorithm is that it is guaranteed to converge to the “true solution” under mild assumptions, despite being based on a particle filter with only a finite number of particles. We will also provide a motivating numerical example illustrating the method using a modeling language tailored for sequential Monte Carlo methods. The intention of modeling languages of this kind is to open up the power of sophisticated Monte Carlo methods—including particle Metropolis–Hastings—to a large group of users without requiring them to know all the underlying mathematical details.
Anytime Monte Carlo
http://www.indii.org/research/anytime-monte-carlo/index.html
Sat, 10 Dec 2016 00:00:00 +0100Lawrence Murrayhttp://www.indii.org/research/anytime-monte-carlo/index.htmlA Monte Carlo algorithm typically simulates some prescribed number of samples, taking some random real time to complete the computations necessary. This work considers the converse: to impose a real-time budget on the computation, so that the number of samples simulated is random. To complicate matters, the real time taken for each simulation may depend on the sample produced, so that the samples themselves are not independent of their number, and a length bias with respect to compute time is apparent. This is especially problematic when a Markov chain Monte Carlo (MCMC) algorithm is used and the final state of the Markov chain—rather than an average over all states—is required. The length bias does not diminish with the compute budget in this case. It occurs, for example, in sequential Monte Carlo (SMC) algorithms. We propose an anytime framework to address the concern, using a continuous-time Markov jump process to study the progress of the computation in real time. We show that the length bias can be eliminated for any MCMC algorithm by using a multiple chain construction. The utility of this construction is demonstrated on a large-scale SMC$^2$ implementation, using four billion particles distributed across a cluster of 128 graphics processing units on the Amazon EC2 service. The anytime framework imposes a real-time budget on the MCMC move steps within SMC$^2$, ensuring that all processors are simultaneously ready for the resampling step, demonstrably reducing wait times and providing substantial control over the total compute budget.
Comparative Analysis of Dengue and Zika Outbreaks Reveals Differences by Setting and Virus
http://www.indii.org/research/comparative-analysis-of-dengue-and-zika-outbreaks-reveals-differences-by-setting-and-virus/index.html
Wed, 07 Dec 2016 00:00:00 +0100Lawrence Murrayhttp://www.indii.org/research/comparative-analysis-of-dengue-and-zika-outbreaks-reveals-differences-by-setting-and-virus/index.htmlThe pacific islands of Micronesia have experienced several outbreaks of mosquito-borne diseases over the past decade. In outbreaks on small islands, the susceptible population is usually well defined, and there is no co-circulation of pathogens. Because of this, analysing such outbreaks can be useful for understanding the transmission dynamics of the pathogens involved, and particularly so for yet understudied pathogens such as Zika virus. Here, we compared three outbreaks of dengue and Zika virus in two different island settings in Micronesia, the Yap Main Islands and Fais, using a mathematical model of transmission dynamics and making full use of commonalities in disease and setting between the outbreaks. We found that the estimated reproduction numbers for Zika and dengue were similar when considered in the same setting, but that, conversely, reproduction number for the same disease can vary considerably by setting. On the Yap Main Islands, we estimated a reproduction number of 8.0–16 (95% Credible Interval (CI)) for the dengue outbreak and 4.8–14 (95% CI) for the Zika outbreak, whereas for the dengue outbreak on Fais our estimate was 28–102 (95% CI). We further found that the proportion of cases of Zika reported was smaller (95% CI 1.4%–1.9%) than that of dengue (95% CI: 47%–61%). We confirmed these results in extensive sensitivity analysis. They suggest that models for dengue transmission can be useful for estimating the predicted dynamics of Zika transmission, but care must be taken when extrapolating findings from one setting to another.
Sequential Monte Carlo with Highly Informative Observations
http://www.indii.org/research/sequential-monte-carlo-with-highly-informative-observations/index.html
Sun, 01 Mar 2015 00:00:00 +0100Lawrence Murrayhttp://www.indii.org/research/sequential-monte-carlo-with-highly-informative-observations/index.htmlWe propose sequential Monte Carlo (SMC) methods for sampling the posterior distribution of state-space models under highly informative observation regimes, a situation in which standard SMC methods can perform poorly. A special case is simulating bridges between given initial and final values. The basic idea is to introduce a schedule of intermediate weighting and resampling times between observation times, which guide particles towards the final state. This can always be done for continuous-time models, and may be done for discrete-time models under sparse observation regimes; our main focus is on continuous-time diffusion processes. The methods are broadly applicable in that they support multivariate models with partial observation, do not require simulation of the backward transition (which is often unavailable), and, where possible, avoid pointwise evaluation of the forward transition. When simulating bridges, the last characteristic cannot be avoided entirely without concessions, and we suggest an $\epsilon$-ball approach (reminiscent of approximate Bayesian computation) as a workaround. Compared to the bootstrap particle filter, the new methods deliver substantially reduced mean squared error in normalizing constant estimates, even after accounting for execution time. The methods are demonstrated for state estimation with two toy examples, and for parameter estimation (within a particle marginal Metropolis–Hastings sampler) with three applied examples in econometrics, epidemiology, and marine biogeochemistry.
Bayesian State-Space Modelling on High-Performance Hardware Using LibBi
http://www.indii.org/research/bayesian-state-space-modelling-on-high-performance-hardware-using-libbi/index.html
Sun, 01 Feb 2015 00:00:00 +0100Lawrence Murrayhttp://www.indii.org/research/bayesian-state-space-modelling-on-high-performance-hardware-using-libbi/index.htmlLibBi is a software package for state-space modelling and Bayesian inference on modern computer hardware, including multi-core central processing units (CPUs), many-core graphics processing units (GPUs) and distributed-memory clusters of such devices. The software parses a domain-specific language for model specification, then optimises, generates, compiles and runs code for the given model, inference method and hardware platform. In presenting the software, this work serves as an introduction to state-space models and the specialised methods developed for Bayesian inference with them. The focus is on sequential Monte Carlo (SMC) methods such as the particle filter for state estimation, and the particle Markov chain Monte Carlo (PMCMC) and SMC$^2$ methods for parameter estimation. All are well-suited to current computer hardware. Two examples are given and developed throughout, one a linear three-element windkessel model of the human arterial system, the other a nonlinear Lorenz ’96 model. These are specified in the prescribed modelling language, and LibBi demonstrated by performing inference with them. Empirical results are presented, including a performance comparison of the software with different hardware configurations.
Parallel Resampling in the Particle Filter
http://www.indii.org/research/parallel-resampling-in-the-particle-filter/index.html
Thu, 01 Jan 2015 00:00:00 +0100Lawrence Murrayhttp://www.indii.org/research/parallel-resampling-in-the-particle-filter/index.htmlModern parallel computing devices, such as the graphics processing
unit (GPU), have gained significant traction in scientific and
statistical computing. They are particularly well-suited to
data-parallel algorithms such as the particle filter, or more
generally Sequential Monte Carlo (SMC), which are increasingly used in
statistical inference. SMC methods carry a set of weighted
particles through repeated propagation, weighting and
resampling steps. The propagation and weighting steps are
straightforward to parallelise, as they require only independent
operations on each particle. The resampling step is more difficult, as
standard schemes require a collective operation, such as a sum, across
particle weights. Focusing on this resampling step, we analyse two
alternative schemes that do not involve a collective operation
(Metropolis and rejection resamplers), and compare them to standard
schemes (multinomial, stratified and systematic resamplers). We find
that, in certain circumstances, the alternative resamplers can perform
significantly faster on a GPU, and to a lesser extent on a CPU, than
the standard approaches. Moreover, in single precision, the standard
approaches are numerically biased for upwards of hundreds of thousands
of particles, while the alternatives are not. This is particularly
important given greater single- than double-precision throughput on
modern devices, and the consequent temptation to use single precision
with a greater number of particles. Finally, we provide auxiliary
functions useful for implementation, such as for the permutation of
ancestry vectors to enable in-place propagation. Supplementary
materials are available online.
Path Storage in the Particle Filter
http://www.indii.org/research/path-storage-in-the-particle-filter/index.html
Mon, 01 Dec 2014 00:00:00 +0100Lawrence Murrayhttp://www.indii.org/research/path-storage-in-the-particle-filter/index.htmlThis article considers the problem of storing the paths generated by a particle filter and more generally by a sequential Monte Carlo algorithm. It provides a theoretical result bounding the expected memory cost by $T+CN\log N$ where $T$ is the time horizon, $N$ is the number of particles and $C$ is a constant, as well as an efficient algorithm to realise this. The theoretical result and the algorithm are illustrated with numerical experiments.
Rethinking soil carbon modelling: a stochastic approach to quantify uncertainties
http://www.indii.org/research/rethinking-soil-carbon-modelling/index.html
Sat, 01 Feb 2014 00:00:00 +0100Lawrence Murrayhttp://www.indii.org/research/rethinking-soil-carbon-modelling/index.htmlThe benefits of sequestering carbon are many, including improved crop productivity, reductions in greenhouse gases, and financial gains through the sale of carbon credits. Achieving better understanding of the sequestration process has motivated many deterministic models of soil carbon dynamics, but none of these models address uncertainty in a comprehensive manner. Uncertainty arises in many ways—around the model inputs, parameters, and dynamics, and subsequently model predictions. In this paper, these uncertainties are addressed in concert by incorporating a physical-statistical model for carbon dynamics within a Bayesian hierarchical modelling framework. This comprehensive approach to accounting for uncertainty in soil carbon modelling has not been attempted previously. This paper demonstrates proof-of-concept based on a one-pool model and identifies requirements for extension to multi-pool carbon modelling. Our model is based on the soil carbon dynamics in Tarlee, South Australia. We specify the model conditionally through its parameters, soil carbon input and decay processes and observations of those processes. We use a particle marginal Metropolis–Hastings approach specified using the LibBi modelling language. We highlight how samples from the posterior distribution can be used to summarise our knowledge about model parameters, to estimate the probabilities of sequestering carbon and to forecast changes in carbon stocks under crop rotations not represented explicitly in the original field trials.