Source-linked AI summary
Fast MCMC sampling for Markov jump processes and extensions
Vinayak Rao, Yee Whye Teh
TL;DR
Posterior inference over MJP paths from partial and noisy observations is difficult to perform analytically. The paper uses an auxiliary-variable Gibbs sampler based on uniformization, alternating virtual-time and path updates. The sampler is exact and is reported to be significantly more efficient than specialized samplers for MMPP and CTBN models.
Problem
Posterior inference over MJP trajectories given noisy observations typically cannot be performed analytically, while existing methods can be costly, approximate, or difficult to extend to complex likelihoods.
Source: §1 · p. 2Method
An auxiliary-variable Gibbs sampler uses uniformization to alternate sampling virtual jump times given the current path with FFBS resampling of the path given the times.
Source: §1 · p. 2, §3.2 · pp. 7–8Results
The sampler converges to the true posterior without approximation and is reported as significantly more efficient than specialized samplers for MMPPs and CTBNs.
Source: §1 · p. 2Takeaways & Limitations
The sampler extends to Markov-modulated Poisson processes and continuous-time Bayesian networks.
Source: §1 · p. 2In its present form, the sampler cannot handle the infinite state space of the Lotka–Volterra example.
Source: §5.3.1 · p. 19
Abstract
from arXiv · showhide
Markov jump processes (or continuous-time Markov chains) are a simple and important class of continuous-time dynamical systems. In this paper, we tackle the problem of simulating from the posterior distribution over paths in these models, given partial and noisy observations. Our approach is an auxiliary variable Gibbs sampler, and is based on the idea of uniformization. This sets up a Markov chain over paths by alternately sampling a finite set of virtual jump times given the current path and then sampling a new path given the set of extant and virtual jump times using a standard hidden Markov model forward filtering-backward sampling algorithm. Our method is exact and does not involve approximations like time-discretization. We demonstrate how our sampler extends naturally to MJP-based models like Markov-modulated Poisson processes and continuous-time Bayesian networks and show significant computational benefits over state-of-the-art MCMC samplers for these models.
1. Introduction
Markov jump processes are widely used to model piecewise-constant system evolution, but posterior inference from noisy observations is generally not analytic. The paper introduces an exact MCMC sampler that avoids costly computations and extends to more complex MJP-based models.
1. Introduction
- MJP paths serve as Bayesian priors for piecewise-constant system evolution across applications including chemistry, biology, genetics, queuing, and finance.The paper identifies the Markov property as supporting realistic models and convenient approximations for complex phenomena.Source: §1 · pp. 1–2
- Existing posterior methods may be generic, require expensive matrix operations or root-finding, introduce bias through discretization or independence assumptions, or struggle with complex likelihoods.The passage cites sampling-based and deterministic approaches as examples.Source: §1 · p. 2
- The sampler avoids expensive computations and approximations, converges to the true posterior, and is reported as significantly more efficient than specialized samplers for MMPPs and CTBNs.It extends to both model classes.Source: §1 · p. 2
- The sampler alternates between virtual-variable sampling conditioned on the current trajectory and trajectory resampling with discrete-time forward-filtering backward-sampling.The auxiliary-variable distribution is independent of observations, which enter through the trajectory-resampling step.Source: §1 · p. 2
2. Markov Jump Processes (MJPs)
This section defines an MJP through its continuous-time, piecewise-constant paths and describes its rate-matrix representation and path-generation procedure. A path can also be represented with a uniformized set of times.
2. Markov Jump Processes (MJPs)
- An MJP has right-continuous, piecewise-constant paths over a countable state space; the paper assumes a finite state space for typical applications.A path is characterized by its initial state, subsequent states, and jump times.Source: §2 · pp. 2, 4
- The uniformized representation pairs a discrete-time state sequence with a set of times, including virtual jumps, to represent the MJP path.Figure 1 contrasts the MJP path with this uniformized representation.Source: Figure 1 · p. 3, §2 · p. 4
- For a homogeneous MJP, the generator A determines the transition matrices P_t; its off-diagonal entries are transition rates and its columns sum to zero.The diagonal entry for each state is the negative total rate of leaving that state.Source: §2 · p. 3
- Gillespie’s algorithm generates a path by drawing exponential waiting times and then selecting new states according to rates from the current state.It stops when the next waiting time would place an event beyond the interval endpoint.Source: §2 · pp. 3–4
3. MCMC Inference via Uniformization
The section addresses posterior sampling of MJP paths under noisy observations, where time discretization introduces bias and can require expensive computation. It presents an exact uniformization-based sampler and discusses its use for MJP parameter inference and experiments.
3.2 The MCMC Algorithm
- Unlike time-discretization, which introduces bias, the sampler targets the true posterior; its stationary distribution is the posterior, and it is irreducible when Ω > max_s |A_s|.Equality is insufficient for irreducibility.Source: §3 · p. 5, §3.2 · p. 9
3.1 Uniformization
- Uniformization represents an MJP using Poisson event times and a subordinated Markov chain whose self-transitions act as virtual jumps.For Ω ≥ max_s |A_s|, this representation defines the same MJP as the original trajectory.Source: §3.1 · pp. 5–6
3.2 The MCMC Algorithm
- The Gibbs sampler resamples virtual jump times given the current path, then uses FFBS to resample states at all Poisson event times and discards virtual jumps.The FFBS step costs O(N^2|V|), and can exploit structure in the rate matrix inherited by the transition matrix.Source: §3.2 · pp. 7–8, Figure 2 · p. 7
3.3 Previous Posterior Sampling Schemes
- Compared with posterior uniformization methods that require observation-specific sampling of Poisson events and often O(N^3) matrix multiplications, this dependent-sample MCMC scales as O(N^2) per random discretization and exploits rate-matrix structure.When off-diagonal transition rates remain O(1) as N grows, the overall sampler scales as O(N^3), versus O(N^4) for matrix exponentiation; the paper reports it is an order of magnitude faster in either rate regime.Source: §3.3 · p. 10
3.4 Bayesian Inference on the MJP Parameters
- For Bayesian parameter inference, the sampler alternates trajectory updates with updates of A using trajectory sufficient statistics; fixing π0 avoids the O(N^3) stationary-distribution calculation.The conjugate prior uses time spent in each state and transition counts; setting π0 to A’s stationary distribution makes the overall sampler cubic despite the quadratic path update.Source: §3.4 · pp. 10–11
3.5 Experiments
- For fixed A, increasing Ω raised mixing but its computational cost outweighed the benefit; with Bayesian inference on A, ESS was fairly insensitive to Ω.The transition count reached equilibrium within a few iterations using Ω = 2 max_s |A_s|.Source: §3.5 · p. 11, Figure 3 · p. 12, Figure 4 · p. 12
4. Markov-Modulated Poisson Processes
The sampler extends to MMPPs by using observed-event counts to compute emission likelihoods. In experiments, it outperformed the Fearnhead and Sherlock sampler across tested settings, with efficiency depending on MJP dimension and initialization.
4. Markov-Modulated Poisson Processes
- For MMPPs, auxiliary jump-event resampling is unchanged, and emission likelihoods require only counting observed Poisson events between successive uniformization times.Fearnhead and Sherlock (2006) instead use a dynamic program with matrix exponentiation over observed events and backward sampling.Source: §4 · p. 13
4.1 Experiments
- Across MMPP experiments, our sampler significantly outperformed Fearnhead and Sherlock (2006); it was largely unaffected by event counts and up to two orders of magnitude faster in other scenarios.The reported exceptions were low-dimensional problems with few Poisson observations over long intervals and one very unstable state.Source: §4.1 · pp. 13–14, Figure 5 · p. 14
- As MJP dimension increased, both samplers scaled cubically, while fixing the initial state distribution made our sampler scale quadratically.With the equilibrium distribution as the initial distribution, the cubic scaling came from computing a rate-matrix eigenvector, not from the trajectory sampler.Source: §4.1 · pp. 14–15, Figure 6 · p. 14
5. Continuous-Time Bayesian Networks (CTBNs)
CTBNs represent structured MJPs through interacting components whose transition rates depend on parent states. The section develops exact posterior sampling for CTBN trajectories and reports computational advantages over the El-Hay sampler across several experiments.
5. Continuous-Time Bayesian Networks (CTBNs)
- A CTBN represents a structured rate matrix with smaller conditional rate matrices; its graph encodes component dependencies and the rate matrices quantify their dynamics.Each node’s rate matrix depends on its current parent configuration.Source: §5 · pp. 15–16
5.1 Inference in CTBNs
- Conditioning a node trajectory on its Markov blanket is difficult because child trajectories act as continuous observations, making the conditional process piecewise inhomogeneous.The resulting inference problem is described as a generalization of MMPP inference.Source: §5.1 · pp. 17–18
5.2 Auxiliary Variable Gibbs Sampling for CTBNs
- The CTBN node update uses a piecewise-constant dominating rate adapted to parent configurations, then forward-filtering backward-sampling with inhomogeneous transition matrices and child-trajectory likelihoods.The auxiliary-variable Gibbs sampler converges to the posterior distribution over CTBN sample paths.Source: §5.2 · pp. 18–19
5.4 Average Relative Error vs Number of Samples
- About 12 times faster, uniformization produced lower average relative error for the same computational effort in this 1,000-run comparison.Figure 10 reports average relative error versus sample count with burn-in of 200.Source: Figure 10 · p. 20, §5.4 · p. 21
5.3.1 THE LOTKA-VOLTERRA PROCESS
- For Lotka–Volterra paths, MCMC better reflected the true predator trajectory than mean-field in the unobserved second half, while mean-field underestimated path variance, especially with few observations.The sampler was tested with each population capped at 200 states because it cannot handle the model’s infinite state space in its present form.Source: §5.3.1 · pp. 19–20
5.4 Average Relative Error vs Number of Samples
- Across the remaining CTBN experiments, uniformization was significantly faster than El-Hay et al.’s sampler, even faster than its matrix-exponentiation time alone.Matrix exponentiation accounted for 10% to 70% of El-Hay et al.’s total running time.Source: §5.4 · p. 20
5.5 Time Requirements for the Chain-Shaped CTBN
- For chain length, uniformization’s time grew linearly and stayed below matrix-exponentiation time; as node states increased, it remained faster and grew more slowly.The authors expect asymptotic costs of O(N^2) for uniformization and O(N^3) for El-Hay et al.Source: Figure 11 · p. 21, §5.5 · p. 22
5.6 Time Requirements for the Drug-Effect CTBN
- As the drug-effect CTBN interval length increased, uniformization remained faster and its computational cost increased linearly.The comparison used fixed standard parameters and states observed at both endpoints; the authors note that El-Hay et al.’s temporal precision setting affects speed.Source: Figure 11 · p. 21, §5.6 · p. 22
6. Discussion
The discussion presents an auxiliary-variable Gibbs sampler for MJPs, using uniformization to alternate between resampling the time discretization and the Markov chain. It reports significant speed-ups for MMPPs and CTBNs, while identifying open questions about tuning and extending the method.
6. Discussion
- The sampler achieves significant speed-ups over state-of-the-art samplers for both MMPPs and continuous-time Bayesian networks.The discussion reports experiments for both model extensions.Source: §6 · p. 22
- Uniformization subordinates a Markov chain to a Poisson process, yielding a random time discretization for a blocked Gibbs sampler that alternates between resampling each component.The method generalizes to Markov-modulated Poisson processes and continuous-time Bayesian networks.Source: §6 · p. 22
- The Poisson rate Ω is a free sampler parameter; the heuristic Ω = max_s 2|A_s| worked in these experiments but may not suit rate matrices with widely varying transition rates.Adaptive MCMC could learn a setting, and the authors leave open how to identify an optimal rate under mixing-speed and computational-cost measures.Source: §6 · p. 22
- The sampler is applicable to inhomogeneous MJPs, including settings where matrix-exponentiation-based techniques cannot be applied.The discussion also proposes extensions to semi-Markov processes with non-exponential holding times and countably infinite state spaces using slice sampling.Source: §6 · pp. 22–23
- For unbounded rate matrices, the predator-prey experiments bounded maximum population sizes; the discussion proposes treating Ω as trajectory-dependent to extend the method beyond that setting.The authors note that trajectories are bounded with probability 1 in practical settings.Source: §6 · p. 23
Appendix A. The Forward-Filtering Backward-Sampling (FFBS) Algorithm
Appendix A describes FFBS for a discrete-time Markov chain with noisy observations. It uses a forward recursion followed by backward sampling to return a posterior sample of the state sequence.
Appendix A. The Forward-Filtering Backward-Sampling (FFBS) Algorithm
- The forward recursion computes filtered quantities across states in O(N^2) per time step, making the full pass O(TN^2).The recursion is defined from the Markov property.Source: Appendix A · p. 23
- After the forward pass samples the final state, the backward pass successively samples states from S_{T-1} to S_0.The result is a sample of the complete state vector.Source: Appendix A · p. 23
- Given an initial distribution, transition matrices, and observation likelihoods, FFBS returns a posterior sample of the Markov-chain state sequence.The chain may be inhomogeneous, with a transition matrix at each time.Source: Appendix A · pp. 23–24