Introduction
Reaction networks are commonly used to describe dynamic processes in chemistry, epidemiology, ecology and biology. Their evolution can be represented by Markov jump processes (MJPs), which are continuous-time discrete-state stochastic processes. Compared to ordinary differential equations or stochastic differential equations (SDEs), MJPs provide a more realistic description of the evolution of such systems when there is a small number of individuals (or molecules) of each species (or reactant).
Given a sequence of noisy observations of the number of individuals of a subset of the species at a discrete set of times, we focus on performing Bayesian inference for the parameters which control the Markov jump process. In other words, we consider Hidden Markov models (HMMs) whose latent processes are MJPs.
Particle Markov chain Monte Carlo (PMCMC) algorithms [1 ] enable inference on HMMs. However, a major issue with these algorithms is particle degeneracy, where over successive resampling steps in the internal particle filter most of the ancestor trajectories get discarded leaving very few unique particles to represent the historical path of the system. For particle Gibbs algorithms, particle degeneracy leads to poor mixing of the states of the latent process at early time points. This occurs especially in moderate- or high-dimensional state spaces and/or when the time series of observations is long. We explain particle degeneracy in more detail in § 2.5 .
One approach that has been introduced to tackle particle degeneracy is particle Gibbs with backward simulation (PGBS) [21 , 22 , 12 ] . This strategy significantly improves the sample efficiency of particle Gibbs, and was extended to particle Gibbs with ancestor sampling (PGAS) by Lindsten et al. [11 ] . However, PGBS and PGAS are only usable when the transition densities of the MJP are tractable. In addition, they only address the particle degeneracy issue caused by long time series; the issue persists when the dimensionality of the state space is high [7 ] .
Exchangeable particle Gibbs (xPG) [13 ] , iterated random-walk conditional sequential Monte Carlo (i-RW-CSMC) [7 ] and Particle-MALA and Particle-mGRAD [2 ] have been introduced to improve the sample efficiency of particle Gibbs in moderate and high-dimensional state spaces. These approaches modify the conditional particle filter, also known as conditional SMC (CSMC), within particle Gibbs to mitigate the particle degeneracy issue. However, i-RW-CSMC, particle-MALA, and particle-mGRAD all require an analytical form of the transition density in the HMM, whereas xPG does not. This represents a significant advantage when the transition density is intractable. For instance, when the latent process evolves as a diffusion in continuous time, the density of the state at a given observation time (conditional on the previous state) typically lacks a closed-form expression. Malory [13 ] applies xPG to such diffusions. By correlating the proposed particles with those in the reference path, xPG generates particles with higher weights. As a result, it is less likely to resample the same particle excessively, thereby mitigating particle degeneracy.
In this paper we develop xPG for HMMs whose latent processes are Markov jump processes (MJPs). We focus on the tau-leap method [8 ] , a time discretisation of the continuous-time MJP. The MJP transition density is difficult to evaluate for two main reasons:
For some reaction networks, such as the susceptible-infectious-susceptible (SIS) model, there is no one-to-one map between a pair of states and the number of each type of reaction required to move from one state to the other.
Even if the count of each type of reaction is uniquely determined when a pair of states is given, evaluating the transition mass function at consecutive observation times requires considering all possible intermediate values between the two states, which can make the computation prohibitively expensive.
We apply xPG to simulated data for different HMMs and compare its performance against PG and, where tractable transitions permit, PGBS and xPG with backward simulation (xPGBS). The mathematical formulation of the models, the particle Gibbs algorithm for MJPs, and the issue of particle degeneracy are introduced in § 2 . § 3 presents the construction of xPG for MJPs, together with a tuning strategy and a backward sampling step for our algorithm. We prove that our algorithm satisfies detailed balance with respect to the intended posterior, and describe the properties relevant to tuning. The algorithm has two key parameters: δ \delta δ M M M δ \delta δ M M M § 3.4.1 explains why we choose this to quantify sample efficiency. In § 4 , we compare the sample efficiency of xPG with PG on data simulated from SIS, susceptible-infectious-recovered (SIR), and susceptible-exposed-infectious-recovered (SEIR) models with univariate and multivariate state spaces, and data simulated from the gene autoregulatory model from Golightly and Wilkinson [10 ] . In the simulation studies of § 4.2 and § 4.3 , xPG produces much larger improvements in efficiency compared with PG than backwards sampling does, but xPGBS is the most efficient of the four algorithms. Section 5 concludes the paper and points out some further research directions.
Notation
We employ the standard notation: a : b : = { a , a + 1 , … , b } a:b := \{a, a+1, \ldots, b\} a : b := { a , a + 1 , … , b } b ≥ a b \ge a b ≥ a a : b : = ∅ a:b := \varnothing a : b := ∅ x a : b x_{a:b} x a : b { x c : c ∈ a : b } \{x_c : c \in a:b\} { x c : c ∈ a : b } Pois ( x ; μ ) \operatorname{Pois}(x;\mu) Pois ( x ; μ ) Bin ( x ; n , p ) \operatorname{Bin}(x;n,p) Bin ( x ; n , p ) μ \mu μ n n n p p p x x x
η , τ \eta, \tau η , τ
T T T T = n × τ T = n \times\tau T = n × τ
L , t 0 : L L, t_{0:L} L , t 0 : L t 0 = 0 t_0 = 0 t 0 = 0 t L = T t_L = T t L = T
M M M
δ ∈ [ 0 , 1 ] \delta\in[0,1] δ ∈ [ 0 , 1 ] δ \delta δ δ = 1 \delta= 1 δ = 1
N N N
N i t e r N_{\mathrm{iter}} N iter
B ℓ m B_\ell^m B ℓ m t ℓ t_\ell t ℓ m m m B L m = m B_L^m = m B L m = m
A ℓ m A_\ell^m A ℓ m t ℓ t_\ell t ℓ x ~ t ℓ m \tilde{x}_{t_\ell}^m x ~ t ℓ m t ℓ t_\ell t ℓ t ℓ t_\ell t ℓ x t ℓ + 1 m x_{t_{\ell+1}}^m x t ℓ + 1 m a ℓ b ℓ + 1 = b ℓ a_\ell^{b_{\ell+1}} = b_\ell a ℓ b ℓ + 1 = b ℓ ℓ = 1 , … , L − 1 \ell= 1,\ldots,L-1 ℓ = 1 , … , L − 1
X t , j ( m ) , i X_{t,j}^{(m),i} X t , j ( m ) , i m m m j j j t t t i i i
− m -m − m { 0 , … , M } ∖ { m } \{0,\ldots,M\} \setminus\{m\} { 0 , … , M } ∖ { m }
a 1 : ℓ − 1 − b 2 : ℓ a_{1:\ell-1}^{-b_{2:\ell}} a 1 : ℓ − 1 − b 2 : ℓ a 1 − b 2 , a 2 − b 3 , … , a ℓ − 1 − b ℓ a_1^{-b_2}, a_2^{-b_3}, \ldots, a_{\ell-1}^{-b_\ell} a 1 − b 2 , a 2 − b 3 , … , a ℓ − 1 − b ℓ
U U U
U − U_- U − U U U
Background
Hidden Markov model
The Hidden Markov model (HMM) is a mathematical framework modelling cases where there are noisy and/or partial observations of an underlying Markov process. The underlying Markov process is referred to as the state process or latent process , while the observations are referred to as the observation process .
If the state process evolves in continuous time, the observations can be made at any time t ∈ [ 0 , T ] t \in[0,T] t ∈ [ 0 , T ] L L L { y 1 , … , y L } \{y_1,\ldots,y_L\} { y 1 , … , y L } t 1 , … , t L t_1,\ldots,t_L t 1 , … , t L t 1 > 0 t_1>0 t 1 > 0 t L = T t_L=T t L = T X 0 ∼ p 0 ( ⋅ ∣ θ ) X_0 \sim p_0(\cdot\mid\theta) X 0 ∼ p 0 ( ⋅ ∣ θ )
X ( t ℓ − 1 , t ℓ ] ∣ ( x [ 0 , t ℓ − 1 ] , y 1 : ℓ − 1 , θ ) ∼ p ( ⋅ ∣ x t ℓ − 1 , θ ) , X_{(t_{\ell-1},t_\ell]} \mid\left(x_{[0,t_{\ell-1}]},y_{1:\ell-1},\theta\right) \sim p\left(\cdot\mid x_{t_{\ell-1}},\theta\right), X ( t ℓ − 1 , t ℓ ] ∣ ( x [ 0 , t ℓ − 1 ] , y 1 : ℓ − 1 , θ ) ∼ p ( ⋅ ∣ x t ℓ − 1 , θ ) , Y ℓ ∣ ( x [ 0 , t ℓ ] , y 1 : ℓ − 1 , θ ) ∼ g ( ⋅ ∣ x t ℓ , θ ) . Y_\ell\mid\left(x_{[0,t_\ell]},y_{1:\ell-1},\theta\right) \sim g\left(\cdot\mid x_{t_\ell},\theta\right). Y ℓ ∣ ( x [ 0 , t ℓ ] , y 1 : ℓ − 1 , θ ) ∼ g ( ⋅ ∣ x t ℓ , θ ) . where ℓ = 1 , … , L \ell=1,\ldots,L ℓ = 1 , … , L p 0 p_0 p 0 θ \theta θ ( t ℓ − 1 , t ℓ ] (t_{\ell-1},t_\ell] ( t ℓ − 1 , t ℓ ] p p p t ℓ − 1 t_{\ell-1} t ℓ − 1 θ \theta θ x t ℓ x_{t_\ell} x t ℓ g g g θ \theta θ
Reaction networks
Reaction networks have applications across a variety of settings. For generality, they will be explained in terms of computational biology, using the terms “molecule ” and “species ”. For an epidemic model, these could become “individual ” and “infection status ”.
Suppose there are u u u v v v R i \mathcal{R}_i R i
∑ j = 1 u a i j X j → h i ∑ j = 1 u b i j X j , i = 1 , … , v . \sum_{j=1}^{u} a_{ij}X_j \xrightarrow{h_i} \sum_{j=1}^{u} b_{ij}X_j,\qquad i=1,\ldots,v. j = 1 ∑ u a ij X j h i j = 1 ∑ u b ij X j , i = 1 , … , v . Here, X j \mathcal{X}_j X j j j j a i j a_{ij} a ij b i j b_{ij} b ij stoichiometric coefficients , representing the number of molecules of X j \mathcal{X}_j X j R i \mathcal{R}_i R i h i h_i h i hazard rate , which governs the frequency at which individual events of reaction R i \mathcal{R}_i R i R i \mathcal{R}_i R i t t t X j \mathcal{X}_j X j b i j − a i j b_{ij}-a_{ij} b ij − a ij j = 1 , … , u j=1,\ldots,u j = 1 , … , u X t = X t − + Γ i X_t=X_{t-}+\Gamma^i X t = X t − + Γ i X t − X_{t-} X t − t t t Γ i = ( b i 1 − a i 1 , … , b i u − a i u ) T \Gamma^i=(b_{i1}-a_{i1},\ldots,b_{iu}-a_{iu})^T Γ i = ( b i 1 − a i 1 , … , b i u − a i u ) T u × v u \times v u × v Γ \Gamma Γ i i i Γ i \Gamma^i Γ i stoichiometry matrix of the reaction network. We now describe the specific reaction networks that will be used in this article; the corresponding stoichiometry matrices are provided in Appendix A.
Susceptible-Infectious-Susceptible (SIS) model
The SIS model assumes two categories of individuals in a population: susceptible (S) and infectious (I). There are two reaction constants, β > 0 \beta>0 β > 0 γ > 0 \gamma>0 γ > 0
R 1 : S + I → β S I / N 2 I (Infection) , R 2 : I → γ I S (Recovery) , \mathcal{R}_1:\mathrm{S}+\mathrm{I}\xrightarrow{\beta SI/N}2\mathrm{I}\quad\text{(Infection)},\qquad\mathcal{R}_2:\mathrm{I}\xrightarrow{\gamma I}\mathrm{S}\quad\text{(Recovery)}, R 1 : S + I β S I / N 2 I (Infection) , R 2 : I γ I S (Recovery) , where S S S I I I
Susceptible-Infectious-Removed (SIR) model
The SIR model incorporates a third species: recovered (R), and has two reaction types:
R 1 : S + I → β S I / N 2 I (Infection) , R 2 : I → γ I R (Recovery) . \mathcal{R}_1:\mathrm{S}+\mathrm{I}\xrightarrow{\beta SI/N}2\mathrm{I}\quad\text{(Infection)},\qquad\mathcal{R}_2:\mathrm{I}\xrightarrow{\gamma I}\mathrm{R}\quad\text{(Recovery)}. R 1 : S + I β S I / N 2 I (Infection) , R 2 : I γ I R (Recovery) . Susceptible-Exposed-Infectious-Recovered (SEIR) model
The SEIR model includes an additional species, exposed (E) individuals, who are infected but not yet infectious:
R 1 : S + I → β S I / N E + I (Infection) , R 2 : E → η E I (Progression to infectious) , \mathcal{R}_1:\mathrm{S}+\mathrm{I}\xrightarrow{\beta SI/N}\mathrm{E}+\mathrm{I}\quad\text{(Infection)},\qquad\mathcal{R}_2:\mathrm{E}\xrightarrow{\eta E}\mathrm{I}\quad\text{(Progression to infectious)}, R 1 : S + I β S I / N E + I (Infection) , R 2 : E η E I (Progression to infectious) , R 3 : I → γ I R (Recovery) . \mathcal{R}_3:\mathrm{I}\xrightarrow{\gamma I}\mathrm{R}\quad\text{(Recovery)}. R 3 : I γ I R (Recovery) . Coupled epidemic models
Suppose a pandemic is spreading across d d d j j j i i i d × d d \times d d × d C C C C i j C_{ij} C ij j j j i i i i , j = 1 , … , d i,j=1,\ldots,d i , j = 1 , … , d i i i
h i ( S , I ) = β S i ∑ j = 1 d C i j I j N i , i = 1 , … , d . h_i(S,I)=\beta S_i\frac{\sum_{j=1}^{d}C_{ij}I_j}{N_i}, \qquad i=1,\ldots,d. h i ( S , I ) = β S i N i ∑ j = 1 d C ij I j , i = 1 , … , d . where S j S_j S j I j I_j I j j j j N j N_j N j j j j j = 1 , … , d j=1,\ldots,d j = 1 , … , d
Autoregulatory model
Golightly and Wilkinson [10 ] models a gene regulation system in which DNA makes RNA, and RNA makes a protein P. Two protein P molecules may form a protein pair, P 2 P_2 P 2 P 2 P_2 P 2 § 6
Markov jump processes
A Markov jump process (MJP) is a continuous-time, discrete-state-space process. It models the stochastic evolution of a reaction network.
An MJP may be simulated exactly using Gillespie’s direct method [8 ] . We focus on a popular approximation of the evolution of an MJP: the tau-leap discretisation [9 ] , given in Algorithm 1. After dividing the continuous time window into small time steps of length τ \tau τ t j t_j t j τ \tau τ R i \mathcal{R}_i R i h i h_i h i R i \mathcal{R}_i R i τ \tau τ
Algorithm 1 Tau-leap method
Divide the full observation interval into time steps of length τ \tau τ
Set t = 0 t=0 t = 0 x 0 = ( x 0 , 1 , … , x 0 , u ) ′ x_0=(x_{0,1},\ldots,x_{0,u})^\prime x 0 = ( x 0 , 1 , … , x 0 , u ) ′
Calculate the hazard rate of each reaction, h i ( x t ) h_i(x_t) h i ( x t ) i = 1 , … , v i=1,\ldots,v i = 1 , … , v
Simulate the number of each reaction R i \mathcal{R}_i R i [ t , t + τ ] [t,t+\tau] [ t , t + τ ] N R i ∼ Poisson ( h i ( x t ) τ ) N^{\mathcal{R}_i}\sim\operatorname{Poisson}(h_i(x_t)\tau) N R i ∼ Poisson ( h i ( x t ) τ )
Set x t + τ : = x t + ∑ i = 1 v N R i s i x_{t+\tau}:=x_t+\sum_{i=1}^{v}N^{\mathcal{R}_i}s_i x t + τ := x t + ∑ i = 1 v N R i s i
Set t : = t + τ t:=t+\tau t := t + τ x t , N R 1 , … , N R v x_t,N^{\mathcal{R}_1},\ldots,N^{\mathcal{R}_v} x t , N R 1 , … , N R v
If t < T t<T t < T
The step size τ \tau τ τ \tau τ τ \tau τ K K K t ℓ = ℓ K τ t_\ell=\ell K\tau t ℓ = ℓ K τ T = L K τ T=LK\tau T = L K τ
Particle Gibbs for MJPs
In this context, a path refers to the collection of latent states { x k τ } k = 0 n \{x_{k\tau}\}_{k=0}^{n} { x k τ } k = 0 n n n n n = K L n=KL n = K L conditional particle filter is conditioned, and from which proposed paths are generated through the resampling steps.
Particle Gibbs alternates sampling the latent process and the model parameter from their conditional posterior distributions, p ( x ( 0 : n ) τ ∣ y 1 : L , θ ) p(x_{(0:n)\tau}\mid y_{1:L},\theta) p ( x ( 0 : n ) τ ∣ y 1 : L , θ ) p ( θ ∣ y 1 : L , x ( 0 : n ) τ ) p(\theta\mid y_{1:L},x_{(0:n)\tau}) p ( θ ∣ y 1 : L , x ( 0 : n ) τ ) θ \theta θ
In each iteration, it first generates M M M { x k τ ( 0 ) } k = 0 n \{x_{k\tau}^{(0)}\}_{k=0}^{n} { x k τ ( 0 ) } k = 0 n M M M M + 1 M+1 M + 1
The conditional particle filter samples M M M x 0 ( 1 : M ) x_0^{(1:M)} x 0 ( 1 : M ) q 0 q_0 q 0 x 0 ( m ) x_0^{(m)} x 0 ( m ) t 1 = K τ t_1=K\tau t 1 = K τ τ \tau τ q τ q_\tau q τ
At each observation time t ℓ t_\ell t ℓ ℓ = 1 , … , L \ell=1,\ldots,L ℓ = 1 , … , L
w ℓ ( m ) = [ p 0 ( x 0 ( m ) ∣ θ ) q 0 ( x 0 ( m ) ∣ θ ) ] I { ℓ = 1 } p ( x t ℓ − 1 + τ ( 1 : K ) ( m ) ∣ x ~ t ℓ − 1 , θ ) g ( y ℓ ∣ x t ℓ − 1 + τ ( 1 : K ) ( m ) , θ ) q ( x t ℓ − 1 + τ ( 1 : K ) ( m ) ∣ x ~ t ℓ − 1 , θ ) , m = 0 , … , M , (1) w_\ell^{(m)} = \left[\frac{p_0(x_0^{(m)}\mid\theta)}{q_0(x_0^{(m)}\mid\theta)}\right]^{\mathbb{I}\{\ell=1\}}\frac{p\left(x_{t_{\ell-1}+\tau(1:K)}^{(m)}\mid\widetilde{x}_{t_{\ell-1}},\theta\right)g\left(y_\ell\mid x_{t_{\ell-1}+\tau(1:K)}^{(m)},\theta\right)}{q\left(x_{t_{\ell-1}+\tau(1:K)}^{(m)}\mid\widetilde{x}_{t_{\ell-1}},\theta\right)},\quad m=0,\ldots,M, \tag*{(1)} w ℓ ( m ) = [ q 0 ( x 0 ( m ) ∣ θ ) p 0 ( x 0 ( m ) ∣ θ ) ] I { ℓ = 1 } q ( x t ℓ − 1 + τ ( 1 : K ) ( m ) ∣ x t ℓ − 1 , θ ) p ( x t ℓ − 1 + τ ( 1 : K ) ( m ) ∣ x t ℓ − 1 , θ ) g ( y ℓ ∣ x t ℓ − 1 + τ ( 1 : K ) ( m ) , θ ) , m = 0 , … , M , (1) where t 0 = 0 t_0=0 t 0 = 0 x ~ 0 ( m ) = x 0 ( m ) \widetilde{x}_0^{(m)}=x_0^{(m)} x 0 ( m ) = x 0 ( m ) I { ℓ = 1 } \mathbb{I}\{\ell=1\} I { ℓ = 1 } ℓ = 1 \ell=1 ℓ = 1
q ( x t ℓ − 1 + τ ( 1 : K ) ∣ x t ℓ − 1 , θ ) = ∏ k = 1 K q τ ( x t ℓ − 1 + k τ ∣ x t ℓ − 1 + ( k − 1 ) τ , θ ) q\left(x_{t_{\ell-1}+\tau(1:K)}\mid x_{t_{\ell-1}},\theta\right)=\prod_{k=1}^{K}q_\tau\left(x_{t_{\ell-1}+k\tau}\mid x_{t_{\ell-1}+(k-1)\tau},\theta\right) q ( x t ℓ − 1 + τ ( 1 : K ) ∣ x t ℓ − 1 , θ ) = k = 1 ∏ K q τ ( x t ℓ − 1 + k τ ∣ x t ℓ − 1 + ( k − 1 ) τ , θ ) is the product of the proposal probability mass functions for the steps of size τ \tau τ q 0 q_0 q 0
Resampling according to the normalised weights yields a sample of unweighted particles { x ~ t ℓ ( m ) } m = 1 M \{\widetilde{x}_{t_\ell}^{(m)}\}_{m=1}^{M} { x t ℓ ( m ) } m = 1 M t L t_L t L
In addition, at resampling time t ℓ t_\ell t ℓ x ~ t ℓ ( m ) \widetilde{x}_{t_\ell}^{(m)} x t ℓ ( m ) a ℓ m a_\ell^m a ℓ m x ~ t ℓ ( m ) = x t ℓ ( a ℓ m ) \widetilde{x}_{t_\ell}^{(m)}=x_{t_\ell}^{(a_\ell^m)} x t ℓ ( m ) = x t ℓ ( a ℓ m ) m = 1 , … , M m=1,\ldots,M m = 1 , … , M ℓ = 1 , … , L \ell=1,\ldots,L ℓ = 1 , … , L a ℓ 0 = 0 a_\ell^0=0 a ℓ 0 = 0 ℓ = 1 , … , L \ell=1,\ldots,L ℓ = 1 , … , L M + 1 M+1 M + 1 x t L ( m ) x_{t_L}^{(m)} x t L ( m )
x ( 0 : n ) τ ( b m 0 : L ) = ( x 0 ( b m 0 ) , x t 0 + τ ( b m 0 ) , … , x t 1 ( b m 1 ) , … , x t ℓ − 1 + τ ( b m ℓ − 1 ) , … , x t ℓ ( b m ℓ ) , … , x t L − 1 + τ ( b m L − 1 ) , … , x t L ( b m L ) ) , x_{(0:n)\tau}^{(b_m^{0:L})}=\left(x_0^{(b_m^0)},x_{t_0+\tau}^{(b_m^0)},\ldots,x_{t_1}^{(b_m^1)},\ldots,x_{t_{\ell-1}+\tau}^{(b_m^{\ell-1})},\ldots,x_{t_\ell}^{(b_m^\ell)},\ldots,x_{t_{L-1}+\tau}^{(b_m^{L-1})},\ldots,x_{t_L}^{(b_m^L)}\right), x ( 0 : n ) τ ( b m 0 : L ) = ( x 0 ( b m 0 ) , x t 0 + τ ( b m 0 ) , … , x t 1 ( b m 1 ) , … , x t ℓ − 1 + τ ( b m ℓ − 1 ) , … , x t ℓ ( b m ℓ ) , … , x t L − 1 + τ ( b m L − 1 ) , … , x t L ( b m L ) ) , where b ℓ m b_\ell^m b ℓ m t ℓ t_\ell t ℓ m m m b L m = m b_L^m=m b L m = m
Within the conditional particle filter, if we initialise from the prior and propagate forward by directly sampling from the transition distribution p p p (1) reduces to those of the bootstrap particle filter, the likelihood of the corresponding observation. Consequently, the closed form of the transition density p ( x t ℓ ∣ x t ℓ − 1 ) p(x_{t_\ell}\mid x_{t_{\ell-1}}) p ( x t ℓ ∣ x t ℓ − 1 ) § 3 . The particle Gibbs algorithm is shown in Algorithm 2.
Particle degeneracy
The multiple proposed paths generated by the conditional particle filter share fewer and fewer distinct ancestor trajectories when traced backwards from the terminal time T T T particle degeneracy . It is illustrated in Figure 5 in Appendix C, and becomes particularly problematic in high-dimensional state spaces and/or when the time horizon is long. In the plot, we generate 100 paths based on the reference path in a single run of the conditional particle filter. However, due to the resampling steps, when tracing these 100 paths backwards from the particles at the terminal time, most share the same ancestors as the reference path during the early part of the time interval. We describe this phenomenon as the paths coalescing backwards in time with the reference path.
Figure 5. Trajectories traced back from the particles at the terminal time in a single run of the standard conditional particle filter with M = 100 M=100 M = 100
Particle Gibbs with Backward Sampling
Backward sampling is an extension to particle Gibbs which seeks to counter particle degeneracy by allowing a point at time t t t t − 1 t-1 t − 1 b L m = m b_L^m=m b L m = m
Algorithm 2 Particle Gibbs
Require: the prior distribution of θ \theta θ π 0 \pi_0 π 0 θ \theta θ p ( θ ∣ y t 1 : L , x ( 0 : n ) τ ( 0 ) ) p(\theta\mid y_{t_{1:L}}, x_{(0:n)\tau}^{(0)}) p ( θ ∣ y t 1 : L , x ( 0 : n ) τ ( 0 ) ) q ∗ q_* q ∗ U − U_- U − N i t e r N_{\mathrm{iter}} N iter
Initialise x ( 0 : n ) τ ( 0 ) , 0 x_{(0:n)\tau}^{(0),0} x ( 0 : n ) τ ( 0 ) , 0
for i = 1 i = 1 i = 1 N i t e r N_{\mathrm{iter}} N iter do
Sample θ ∼ p ( ⋅ ∣ y t 1 : L , x ( 0 : n ) τ ( 0 ) , i − 1 ) \theta\sim p(\cdot\mid y_{t_{1:L}},x_{(0:n)\tau}^{(0),i-1}) θ ∼ p ( ⋅ ∣ y t 1 : L , x ( 0 : n ) τ ( 0 ) , i − 1 )
Sample u − ∼ q ∗ ( ⋅ ∣ θ , x ( 0 : n ) τ ( 1 ) , i − 1 ) u_- \sim q_*(\cdot\mid\theta,x_{(0:n)\tau}^{(1),i-1}) u − ∼ q ∗ ( ⋅ ∣ θ , x ( 0 : n ) τ ( 1 ) , i − 1 ) p ( x ( 0 : n ) τ ∣ θ , y t 1 : L ) p(x_{(0:n)\tau}\mid\theta,y_{t_{1:L}}) p ( x ( 0 : n ) τ ∣ θ , y t 1 : L ) x ( 0 : n ) τ ( 1 : M ) x_{(0:n)\tau}^{(1:M)} x ( 0 : n ) τ ( 1 : M ) a 1 : ( L − 1 ) 1 : M a_{1:(L-1)}^{1:M} a 1 : ( L − 1 ) 1 : M w L 0 : M w_L^{0:M} w L 0 : M
Set u = { u − , x ( 0 : n ) τ ( 0 ) , i − 1 } u = \{u_-,x_{(0:n)\tau}^{(0),i-1}\} u = { u − , x ( 0 : n ) τ ( 0 ) , i − 1 }
Sample an index m m m { 0 , 1 , … , M } \{0,1,\ldots,M\} { 0 , 1 , … , M } w L m w_L^m w L m
Set θ i = θ \theta^i=\theta θ i = θ x ( 0 : n ) τ ( 0 ) , i = x ( 0 : n ) τ ( b 1 m : L ) x_{(0:n)\tau}^{(0),i}=x_{(0:n)\tau}^{(b_1^{m:L})} x ( 0 : n ) τ ( 0 ) , i = x ( 0 : n ) τ ( b 1 m : L )
end for
return { θ i } i = 1 N i t e r \{\theta^i\}_{i=1}^{N_{\mathrm{iter}}} { θ i } i = 1 N iter { x ( 0 : n ) τ ( 0 ) , i } i = 0 N i t e r \{x_{(0:n)\tau}^{(0),i}\}_{i=0}^{N_{\mathrm{iter}}} { x ( 0 : n ) τ ( 0 ) , i } i = 0 N iter
obtained by sampling the index backward as follows. For ℓ = L − 1 , … , 1 \ell= L - 1,\ldots,1 ℓ = L − 1 , … , 1 b ℓ m b_\ell^m b ℓ m P ( B ℓ m = i ) ∝ w ~ ℓ i \mathbb{P}(B_\ell^m=i)\propto\widetilde{w}_\ell^i P ( B ℓ m = i ) ∝ w ℓ i
w ~ ℓ i = w ℓ i ⋅ p ( x t ℓ + 1 ( b ℓ + 1 m ) ∣ x t ℓ ( i ) , θ ) , i = 0 , … , M . (2) \widetilde{w}_\ell^i=w_\ell^i\cdot p\left(x_{t_{\ell+1}}^{(b_{\ell+1}^m)}\mid x_{t_\ell}^{(i)},\theta\right),\quad i=0,\ldots,M. \tag*{(2)} w ℓ i = w ℓ i ⋅ p ( x t ℓ + 1 ( b ℓ + 1 m ) ∣ x t ℓ ( i ) , θ ) , i = 0 , … , M . (2) Then we obtain a new reference path, x ( 0 : n ) τ ( b 1 m : L ) x_{(0:n)\tau}^{(b_1^{m:L})} x ( 0 : n ) τ ( b 1 m : L )
By doing this, the particle degeneracy issue is largely alleviated. However, this requires the transition probabilities from one observation time to the next to be tractable. They are generally intractable when the corresponding reaction counts are not uniquely determined by the latent state at the two times, or when the interval between observations contains multiple intermediate time steps.
Exchangeable particle Gibbs for diffusions
To address the particle degeneracy issue, Malory [13 ] introduced a method called exchangeable particle Gibbs (xPG). This method is intended for discrete-time models with continuous state spaces.
Firstly, the latent processes were intractable, multidimensional diffusions { X t } t ∈ [ 0 , T ] \{X_t\}_{t\in[0,T]} { X t } t ∈ [ 0 , T ] X t X_t X t t t t X ⊆ R d \mathcal{X}\subseteq\mathbb{R}^d X ⊆ R d d X t = μ ( X t ) d t + Σ ( X t ) d W t dX_t=\mu(X_t)dt+\sqrt{\Sigma(X_t)}dW_t d X t = μ ( X t ) d t + Σ ( X t ) d W t X t = x X_t=x X t = x
X t + τ = x + μ ( x ) τ + Σ ( x ) Δ W t (3) X_{t+\tau}=x+\mu(x)\tau+\sqrt{\Sigma(x)}\Delta W_t \tag*{(3)} X t + τ = x + μ ( x ) τ + Σ ( x ) Δ W t (3) where Δ W t = τ Z t \Delta W_t=\sqrt{\tau}Z_t Δ W t = τ Z t Z t ∼ N ( 0 , 1 ) Z_t\sim N(0,1) Z t ∼ N ( 0 , 1 )
The key idea of xPG applied to SDEs is to correlate both the initial state X 0 X_0 X 0 π ( x [ 0 , T ] ∣ y 1 : T , θ ) \pi(x_{[0,T]}\mid y_{1:T},\theta) π ( x [ 0 , T ] ∣ y 1 : T , θ ) M M M π ( x [ 0 , T ] ∣ θ ) \pi(x_{[0,T]}\mid\theta) π ( x [ 0 , T ] ∣ θ ) M + 1 M+1 M + 1 M + 1 M+1 M + 1 θ \theta θ
Given the initial state of the reference path x 0 ( 0 ) x_0^{(0)} x 0 ( 0 )
q ~ 0 ( x 0 ( 1 : M ) ∣ x 0 ( 0 ) ) = q 0 ( M ) ( x 0 ( 0 : M ) ) q 0 ( x 0 ( 0 ) ) , (4) \widetilde{q}_0(x_0^{(1:M)}\mid x_0^{(0)})=\frac{q_0^{(M)}(x_0^{(0:M)})}{q_0(x_0^{(0)})}, \tag*{(4)} q 0 ( x 0 ( 1 : M ) ∣ x 0 ( 0 ) ) = q 0 ( x 0 ( 0 ) ) q 0 ( M ) ( x 0 ( 0 : M ) ) , (4) where q 0 ( M ) q_0^{(M)} q 0 ( M ) q 0 ≡ p 0 q_0\equiv p_0 q 0 ≡ p 0 σ \sigma σ 0 : M 0:M 0 : M
q 0 ( M ) ( x 0 ( 0 : M ) ) = q 0 ( M ) ( x 0 ( σ ( 0 : M ) ) ) . (5) q_0^{(M)}\left(x_0^{(0:M)}\right)=q_0^{(M)}\left(x_0^{(\sigma(0:M))}\right). \tag*{(5)} q 0 ( M ) ( x 0 ( 0 : M ) ) = q 0 ( M ) ( x 0 ( σ ( 0 : M )) ) . (5) Moreover, for the simulation from time t − τ t-\tau t − τ t t t t = τ , 2 τ , … , K τ t=\tau,2\tau,\ldots,K\tau t = τ , 2 τ , … , K τ M M M x t ( 1 ) , … , x t ( M ) x_t^{(1)},\ldots,x_t^{(M)} x t ( 1 ) , … , x t ( M ) q ~ \widetilde{q} q x t − τ ( 0 : M ) x_{t-\tau}^{(0:M)} x t − τ ( 0 : M ) x t ( 0 ) x_t^{(0)} x t ( 0 )
q ~ ( x t ( 1 : M ) ∣ x t ( 0 ) , x t − τ ( 0 : M ) ) = q ( M ) ( x t ( 0 : M ) ∣ x t − τ ( 0 : M ) ) q ( x t ( 0 ) ∣ x t − τ ( 0 ) ) (6) \widetilde{q}\left(x_t^{(1:M)}\mid x_t^{(0)},x_{t-\tau}^{(0:M)}\right)=\frac{q^{(M)}\left(x_t^{(0:M)}\mid x_{t-\tau}^{(0:M)}\right)}{q\left(x_t^{(0)}\mid x_{t-\tau}^{(0)}\right)} \tag*{(6)} q ( x t ( 1 : M ) ∣ x t ( 0 ) , x t − τ ( 0 : M ) ) = q ( x t ( 0 ) ∣ x t − τ ( 0 ) ) q ( M ) ( x t ( 0 : M ) ∣ x t − τ ( 0 : M ) ) (6) where q ( M ) ( ⋅ ∣ x t − τ ( 0 : M ) ) q^{(M)}(\cdot\mid x_{t-\tau}^{(0:M)}) q ( M ) ( ⋅ ∣ x t − τ ( 0 : M ) ) M + 1 M+1 M + 1 { ( x t ( 0 ) , x t − τ ( 0 ) ) , … , ( x t ( M ) , x t − τ ( M ) ) } \{(x_t^{(0)},x_{t-\tau}^{(0)}),\ldots,(x_t^{(M)},x_{t-\tau}^{(M)})\} {( x t ( 0 ) , x t − τ ( 0 ) ) , … , ( x t ( M ) , x t − τ ( M ) )} q ( x t ∣ x t − τ ) ≡ p τ ( x t ∣ x t − τ ) q(x_t\mid x_{t-\tau})\equiv p_\tau(x_t\mid x_{t-\tau}) q ( x t ∣ x t − τ ) ≡ p τ ( x t ∣ x t − τ )
Sampling from both distributions requires a mechanism for obtaining Z = ( Z ( 0 ) , … , Z ( M ) ) ⊤ \mathbf{Z}=(Z^{(0)},\ldots,Z^{(M)})^\top Z = ( Z ( 0 ) , … , Z ( M ) ) ⊤ N ( μ M + 1 , Σ M + 1 ) N(\boldsymbol{\mu}_{M+1},\boldsymbol{\Sigma}_{M+1}) N ( μ M + 1 , Σ M + 1 ) μ M + 1 = 0 ( M + 1 ) d \boldsymbol{\mu}_{M+1}=\mathbf{0}_{(M+1)d} μ M + 1 = 0 ( M + 1 ) d Σ ∈ R ( M + 1 ) d × ( M + 1 ) d \boldsymbol{\Sigma}\in\mathbb{R}^{(M+1)d\times(M+1)d} Σ ∈ R ( M + 1 ) d × ( M + 1 ) d I d I_d I d ( M + 1 ) × ( M + 1 ) (M+1)\times(M+1) ( M + 1 ) × ( M + 1 ) ρ \rho ρ
At the centre of the mechanism is the methodology of Tjelmeland [20 ] for generating correlated standard Gaussian random variables. Imagine Z ( 0 ) ∼ N ( 0 , 1 ) Z^{(0)}\sim N(0,1) Z ( 0 ) ∼ N ( 0 , 1 ) Z ∗ = 1 − δ Z ( 0 ) + Z δ ( 0 ) Z^*=\sqrt{1-\delta}Z^{(0)}+Z_\delta^{(0)} Z ∗ = 1 − δ Z ( 0 ) + Z δ ( 0 ) Z δ ( 0 ) Z_\delta^{(0)} Z δ ( 0 ) N ( 0 , δ 2 ) N(0,\delta^2) N ( 0 , δ 2 ) m = 1 , … , M m=1,\ldots,M m = 1 , … , M Z ( m ) = Z ∗ + Z δ ( m ) Z^{(m)}=Z^*+Z_\delta^{(m)} Z ( m ) = Z ∗ + Z δ ( m ) Z δ ( 1 : M ) Z_\delta^{(1:M)} Z δ ( 1 : M ) ∼ N ( 0 , δ 2 ) \sim N(0,\delta^2) ∼ N ( 0 , δ 2 ) x t ( 1 : M ) x_t^{(1:M)} x t ( 1 : M ) x t ( 0 ) x_t^{(0)} x t ( 0 ) x t − τ ( 0 : M ) x_{t-\tau}^{(0:M)} x t − τ ( 0 : M )
In Malory [13 ] , this methodology was applied to the setting in which the latent process is a diffusion, driven by Gaussian increments. Our contribution is the far-from trivial extension of this method to the case when the latent process is a discretised Markov jump process, driven by Poisson increments.
xPGibbs for MJPs
Main ideas
In Malory [13 ] , the reference and proposed initial values and the Gaussian Brownian increments were made exchangeable with respect to their prior distribution by correlating the unit Gaussians on which they were based. This meant that the paths were exchangeable with marginals given by their prior distribution. Here, we describe methodology for making the initial values and the Poisson increments over each time step exchangeable between the reference path and the proposals. As in Malory [13 ] the marginals are given by the priors; that is, in (4) and (6) , q 0 ( x 0 ) = p 0 ( x 0 ) q_0(x_0)=p_0(x_0) q 0 ( x 0 ) = p 0 ( x 0 ) q ( x t ∣ x t − τ ) = p τ ( x t ∣ x t − τ ) q(x_t\mid x_{t-\tau})=p_\tau(x_t\mid x_{t-\tau}) q ( x t ∣ x t − τ ) = p τ ( x t ∣ x t − τ )
The key steps in this process are jointly simulating:
Initial values for the M M M q 0 ( x 0 ) = p 0 ( x 0 ) q_0(x_0)=p_0(x_0) q 0 ( x 0 ) = p 0 ( x 0 )
The number of events of a particular reaction over a single interval of length τ \tau τ M M M M M M
The second task is the more difficult, and so we tackle that first.
Simulating event numbers
We divide the continuous time interval from 0 0 0 T T T τ \tau τ R r \mathcal{R}_r R r m m m [ ( k − 1 ) τ , k τ ] [(k-1)\tau,k\tau] [( k − 1 ) τ , k τ ] N k ( m ) , R r N_k^{(m),\mathcal{R}_r} N k ( m ) , R r N k ( m ) , R r N_k^{(m),\mathcal{R}_r} N k ( m ) , R r { N k ( m ) , R s ; s ≠ r } \{N_k^{(m),\mathcal{R}_s};s\ne r\} { N k ( m ) , R s ; s = r } R r \mathcal{R}_r R r μ k ( m ) \mu_k^{(m)} μ k ( m ) R r \mathcal{R}_r R r k k k
At time point ( k − 1 ) τ (k-1)\tau ( k − 1 ) τ N k ( 1 ) , … , N k ( M ) ∣ N k ( 0 ) , μ k ( 0 ) , … , μ k ( M ) N_k^{(1)},\ldots,N_k^{(M)}\mid N_k^{(0)},\mu_k^{(0)},\ldots,\mu_k^{(M)} N k ( 1 ) , … , N k ( M ) ∣ N k ( 0 ) , μ k ( 0 ) , … , μ k ( M ) q ~ ( x k ( 1 : M ) ∣ x k ( 0 ) , x ( k − 1 ) τ ( 0 : M ) ) \widetilde{q}(x_k^{(1:M)}\mid x_k^{(0)},x_{(k-1)\tau}^{(0:M)}) q ( x k ( 1 : M ) ∣ x k ( 0 ) , x ( k − 1 ) τ ( 0 : M ) ) (6) .
We first thin N k ( 0 ) N_k^{(0)} N k ( 0 ) N ~ k ( 0 ) \widetilde{N}_k^{(0)} N k ( 0 )
N ~ k ( 0 ) ∣ N k ( 0 ) ∼ Bin ( N k ( 0 ) , 1 − δ ) (7) \widetilde{N}_k^{(0)}\mid N_k^{(0)}\sim\operatorname{Bin}(N_k^{(0)},1-\delta) \tag*{(7)} N k ( 0 ) ∣ N k ( 0 ) ∼ Bin ( N k ( 0 ) , 1 − δ ) (7) If N k ( 0 ) ∼ Pois ( μ k ( 0 ) ) N_k^{(0)} \sim\operatorname{Pois}(\mu_k^{(0)}) N k ( 0 ) ∼ Pois ( μ k ( 0 ) ) N ~ k ( 0 ) ∼ Pois ( ( 1 − δ ) μ k ( 0 ) ) \widetilde{N}_k^{(0)} \sim\operatorname{Pois}((1-\delta)\mu_k^{(0)}) N k ( 0 ) ∼ Pois (( 1 − δ ) μ k ( 0 ) ) N ~ k ( 0 ) \widetilde{N}_k^{(0)} N k ( 0 ) N ~ k ( 0 ) ∼ Pois ( μ k ( 0 ) ) \widetilde{N}_k^{(0)} \sim\operatorname{Pois}(\mu_k^{(0)}) N k ( 0 ) ∼ Pois ( μ k ( 0 ) ) N ~ k ( 0 : M ) \widetilde{N}_k^{(0:M)} N k ( 0 : M ) 1 − δ 1-\delta 1 − δ μ k ( 0 : M ) \mu_k^{(0:M)} μ k ( 0 : M ) μ k ( s − l o ) ≤ ⋯ ≤ μ k ( s − 1 ) ≤ μ k ( s 0 ) ≤ μ k ( s 1 ) ≤ ⋯ ≤ μ k ( s h i ) \mu_k^{(s_{-lo})} \le\cdots\le\mu_k^{(s_{-1})} \le\mu_k^{(s_0)} \le\mu_k^{(s_1)} \le\cdots\le\mu_k^{(s_{hi})} μ k ( s − l o ) ≤ ⋯ ≤ μ k ( s − 1 ) ≤ μ k ( s 0 ) ≤ μ k ( s 1 ) ≤ ⋯ ≤ μ k ( s hi ) s 0 = 0 s_0=0 s 0 = 0
N ~ k ( s i + 1 ) ∣ N ~ k ( s i ) = N ~ k ( s i ) + Pois ( ( 1 − δ ) ( μ k ( s i + 1 ) − μ k ( s i ) ) ) , (8) \widetilde{N}_k^{(s_{i+1})} \mid\widetilde{N}_k^{(s_i)} = \widetilde{N}_k^{(s_i)} + \operatorname{Pois}\left((1-\delta)\left(\mu_k^{(s_{i+1})}-\mu_k^{(s_i)}\right)\right), \tag*{(8)} N k ( s i + 1 ) ∣ N k ( s i ) = N k ( s i ) + Pois ( ( 1 − δ ) ( μ k ( s i + 1 ) − μ k ( s i ) ) ) , (8) i = 0 , 1 , … , h i − 1 i=0,1,\ldots,hi-1 i = 0 , 1 , … , hi − 1
N ~ k ( s i − 1 ) ∣ N ~ k ( s i ) ∼ Bin ( N ~ k ( s i ) , μ k ( s i − 1 ) μ k ( s i ) ) , (9) \widetilde{N}_k^{(s_{i-1})} \mid\widetilde{N}_k^{(s_i)} \sim\operatorname{Bin}\left(\widetilde{N}_k^{(s_i)},\frac{\mu_k^{(s_{i-1})}}{\mu_k^{(s_i)}}\right), \tag*{(9)} N k ( s i − 1 ) ∣ N k ( s i ) ∼ Bin ( N k ( s i ) , μ k ( s i ) μ k ( s i − 1 ) ) , (9) i = 0 , − 1 , … , − l o + 1 i=0,-1,\ldots,-lo+1 i = 0 , − 1 , … , − l o + 1
N k ( j ) ∣ N ~ k ( j ) = N ~ k ( j ) + Pois ( δ μ k ( j ) ) j = 1 , … , M . (10) N_k^{(j)} \mid\widetilde{N}_k^{(j)} = \widetilde{N}_k^{(j)}+\operatorname{Pois}\left(\delta\mu_k^{(j)}\right) \qquad j=1,\ldots,M. \tag*{(10)} N k ( j ) ∣ N k ( j ) = N k ( j ) + Pois ( δ μ k ( j ) ) j = 1 , … , M . (10) Therefore, marginally N k ( j ) ∼ Pois ( μ k ( j ) ) N_k^{(j)} \sim\operatorname{Pois}(\mu_k^{(j)}) N k ( j ) ∼ Pois ( μ k ( j ) ) N ~ k ( 0 : M ) ∣ μ k ( 0 : M ) \widetilde{N}_k^{(0:M)}\mid\mu_k^{(0:M)} N k ( 0 : M ) ∣ μ k ( 0 : M ) N ~ k ( σ ( 0 : M ) ) ∣ μ k ( σ ( 0 : M ) ) \widetilde{N}_k^{(\sigma(0:M))}\mid\mu_k^{(\sigma(0:M))} N k ( σ ( 0 : M )) ∣ μ k ( σ ( 0 : M )) μ k ( 0 : M ) \mu_k^{(0:M)} μ k ( 0 : M ) 1 − δ 1-\delta 1 − δ N ~ k ( 0 : M ) \widetilde{N}_k^{(0:M)} N k ( 0 : M ) N k ( 0 : M ) ∣ μ k ( 0 : M ) N_k^{(0:M)}\mid\mu_k^{(0:M)} N k ( 0 : M ) ∣ μ k ( 0 : M ) N k ( σ ( 0 : M ) ) ∣ μ k ( σ ( 0 : M ) ) N_k^{(\sigma(0:M))}\mid\mu_k^{(\sigma(0:M))} N k ( σ ( 0 : M )) ∣ μ k ( σ ( 0 : M )) μ k ( m ) \mu_k^{(m)} μ k ( m ) x ( k − 1 ) τ ( m ) x_{(k-1)\tau}^{(m)} x ( k − 1 ) τ ( m ) x k τ ( m ) x_{k\tau}^{(m)} x k τ ( m ) N k ( m ) N_k^{(m)} N k ( m ) q ( M ) q^{(M)} q ( M )
The correlation between N k ( s i ) N_k^{(s_i)} N k ( s i ) N k ( s j ) N_k^{(s_j)} N k ( s j ) i < j i<j i < j δ \delta δ δ = 1 \delta=1 δ = 1 N k ( 0 ) , … , N k ( M ) N_k^{(0)},\ldots,N_k^{(M)} N k ( 0 ) , … , N k ( M ) δ = 0 \delta=0 δ = 0
Simulating initial values
We now turn to simulating exchangeably from the prior at time 0, conditional on the value in the reference path; i.e., q ~ 0 \widetilde{q}_0 q 0 p 0 p_0 p 0 X 0 = x 0 X_0=x_0 X 0 = x 0 X 0 ∼ Pois ( a 0 ) X_0\sim\operatorname{Pois}(a_0) X 0 ∼ Pois ( a 0 ) X 0 ∼ Bin ( N , p 0 ) X_0\sim\operatorname{Bin}(N,p_0) X 0 ∼ Bin ( N , p 0 )
Case (a) is trivial. For case (b) we first thin X 0 ( 0 ) X_0^{(0)} X 0 ( 0 ) X ~ 0 ∼ Bin ( X 0 , 1 − δ ) \widetilde{X}_0\sim\operatorname{Bin}(X_0,1-\delta) X 0 ∼ Bin ( X 0 , 1 − δ ) X 0 ( 1 : M ) X_0^{(1:M)} X 0 ( 1 : M )
X 0 ( m ) = X ~ 0 + Pois ( δ a 0 ) m = 1 , … , M . (11) X_0^{(m)}=\widetilde{X}_0+\operatorname{Pois}(\delta a_0) \qquad m=1,\ldots,M. \tag*{(11)} X 0 ( m ) = X 0 + Pois ( δ a 0 ) m = 1 , … , M . (11) Clearly X 0 ( 0 : M ) X_0^{(0:M)} X 0 ( 0 : M ) Pois ( a 0 ) \operatorname{Pois}(a_0) Pois ( a 0 )
For case (c), the core is generated by sampling from a hypergeometric distribution
X ~ 0 ∼ Hypergeometric ( N , X 0 , N 0 ) (12) \widetilde{X}_0\sim\operatorname{Hypergeometric}(N,X_0,N_0) \tag*{(12)} X 0 ∼ Hypergeometric ( N , X 0 , N 0 ) (12) where N 0 = ⌊ ( 1 − δ ) N ⌋ N_0=\lfloor(1-\delta)N\rfloor N 0 = ⌊( 1 − δ ) N ⌋ N 0 N_0 N 0 N N N X 0 X_0 X 0 X ~ 0 \widetilde{X}_0 X 0 N 0 N_0 N 0 X ~ 0 ∼ Bin ( N 0 , p 0 ) \widetilde{X}_0\sim\operatorname{Bin}(N_0,p_0) X 0 ∼ Bin ( N 0 , p 0 )
X 0 ( m ) = X ~ 0 + Bin ( N − N 0 , p 0 ) , m = 1 , … , M . X_0^{(m)}=\widetilde{X}_0+\operatorname{Bin}(N-N_0,p_0), \qquad m=1,\ldots,M. X 0 ( m ) = X 0 + Bin ( N − N 0 , p 0 ) , m = 1 , … , M . Therefore, X 0 ( 0 : M ) X_0^{(0:M)} X 0 ( 0 : M ) Bin ( N , p 0 ) \operatorname{Bin}(N,p_0) Bin ( N , p 0 )
In cases (b) and (c) the pairwise correlation is 0 when δ = 1 \delta=1 δ = 1 δ = 0 \delta=0 δ = 0
Exchangeable conditional particle filter for the tau-leap
Given a sequence of observations y t 1 , … , y t L y_{t_1},\ldots,y_{t_L} y t 1 , … , y t L t 1 > 0 t_1>0 t 1 > 0 t L = T t_L=T t L = T q 0 ( x 0 ) = p 0 ( x 0 ) q_0(x_0)=p_0(x_0) q 0 ( x 0 ) = p 0 ( x 0 ) q ( x t ∣ x t − τ ) = p ( x t ∣ x t − τ ) q(x_t\mid x_{t-\tau})=p(x_t\mid x_{t-\tau}) q ( x t ∣ x t − τ ) = p ( x t ∣ x t − τ ) x 0 ( 0 ) , … , x 0 ( M ) x_0^{(0)},\ldots,x_0^{(M)} x 0 ( 0 ) , … , x 0 ( M ) q ~ 0 ( ⋅ ∣ x 0 ( 0 ) ) \widetilde{q}_0(\cdot\mid x_0^{(0)}) q 0 ( ⋅ ∣ x 0 ( 0 ) ) § 3.1.2 , ensuring that x 0 ( 0 ) , x 0 ( 1 ) , … , x 0 ( M ) x_0^{(0)},x_0^{(1)},\ldots,x_0^{(M)} x 0 ( 0 ) , x 0 ( 1 ) , … , x 0 ( M ) q ~ 0 ( M ) \widetilde{q}_0^{(M)} q 0 ( M ) (5) . Starting from the initial states, we update the number of each reaction R i \mathcal{R}_i R i i = 1 , … , v i=1,\ldots,v i = 1 , … , v [ 0 , τ ] [0,\tau] [ 0 , τ ] § 3.1.1 . We then update the states for all the proposed paths at the starting point of the next time step [ τ , 2 τ ] [\tau,2\tau] [ τ , 2 τ ] t 1 t_1 t 1 t 1 t_1 t 1 g ( y t 1 ∣ x t 1 , θ ) g(y_{t_1}\mid x_{t_1},\theta) g ( y t 1 ∣ x t 1 , θ ) M M M { x t 1 ( m ) } m = 0 M \{x_{t_1}^{(m)}\}_{m=0}^{M} { x t 1 ( m ) } m = 0 M T T T
To obtain exchangeable particle Gibbs, we replace the conditional particle filter in Algorithm 2 with the exchangeable conditional particle filter presented in Algorithm 3.
Algorithm 3 Exchangeable conditional particle filter
Require: the parameter, θ \theta θ x ( 0 : n ) τ ( 0 ) x_{(0:n)\tau}^{(0)} x ( 0 : n ) τ ( 0 ) y t 1 : L y_{t_{1:L}} y t 1 : L q 0 ( ⋅ ) q_0(\cdot) q 0 ( ⋅ ) X 0 X_0 X 0 q ( ⋅ ∣ x t − 1 ) q(\cdot\mid x_{t-1}) q ( ⋅ ∣ x t − 1 ) X t X_t X t X t − 1 = x t − 1 X_{t-1}=x_{t-1} X t − 1 = x t − 1
Sample x 0 ( 1 : M ) x_0^{(1:M)} x 0 ( 1 : M ) q ~ 0 ( ⋅ ∣ x 0 ( 0 ) , θ ) \widetilde{q}_0(\cdot\mid x_0^{(0)},\theta) q 0 ( ⋅ ∣ x 0 ( 0 ) , θ )
For k = 1 k=1 k = 1 K K K x k τ ( 1 : M ) x_{k\tau}^{(1:M)} x k τ ( 1 : M ) q ~ ( ⋅ ∣ x k τ ( 0 ) , x ( k − 1 ) τ ( 1 : M ) , θ ) \widetilde{q}(\cdot\mid x_{k\tau}^{(0)},x_{(k-1)\tau}^{(1:M)},\theta) q ( ⋅ ∣ x k τ ( 0 ) , x ( k − 1 ) τ ( 1 : M ) , θ )
Compute the weights of x t 1 ( 0 ) , … , x t 1 ( M ) x_{t_1}^{(0)},\ldots,x_{t_1}^{(M)} x t 1 ( 0 ) , … , x t 1 ( M )
for ℓ = 1 \ell=1 ℓ = 1 L L L do
Set x ~ t ℓ − 1 ( 0 ) = x t ℓ − 1 ( 0 ) \widetilde{x}_{t_{\ell-1}}^{(0)}=x_{t_{\ell-1}}^{(0)} x t ℓ − 1 ( 0 ) = x t ℓ − 1 ( 0 )
Sample M M M { x t ℓ − 1 ( m ) } m = 0 M \{x_{t_{\ell-1}}^{(m)}\}_{m=0}^{M} { x t ℓ − 1 ( m ) } m = 0 M x t ℓ − 1 ( m ) x_{t_{\ell-1}}^{(m)} x t ℓ − 1 ( m ) w t ℓ − 1 ( m ) w_{t_{\ell-1}}^{(m)} w t ℓ − 1 ( m ) { x ~ t ℓ − 1 ( m ) } m = 1 M \{\widetilde{x}_{t_{\ell-1}}^{(m)}\}_{m=1}^{M} { x t ℓ − 1 ( m ) } m = 1 M
for k = 1 k=1 k = 1 K K K do
Sample x t ℓ − 1 + k τ ( 1 : M ) x_{t_{\ell-1}+k\tau}^{(1:M)} x t ℓ − 1 + k τ ( 1 : M ) q ~ ( ⋅ ∣ x t ℓ − 1 + k τ ( 0 ) , x ~ t ℓ − 1 + ( k − 1 ) τ ( 1 : M ) , θ ) \widetilde{q}(\cdot\mid x_{t_{\ell-1}+k\tau}^{(0)},\widetilde{x}_{t_{\ell-1}+(k-1)\tau}^{(1:M)},\theta) q ( ⋅ ∣ x t ℓ − 1 + k τ ( 0 ) , x t ℓ − 1 + ( k − 1 ) τ ( 1 : M ) , θ )
If k < K k<K k < K x ~ t ℓ − 1 + k τ ( m ) = x t ℓ − 1 + k τ ( m ) \widetilde{x}_{t_{\ell-1}+k\tau}^{(m)}=x_{t_{\ell-1}+k\tau}^{(m)} x t ℓ − 1 + k τ ( m ) = x t ℓ − 1 + k τ ( m ) m = 1 , … , M m=1,\ldots,M m = 1 , … , M
end for
Compute the weights of x t ℓ ( 0 ) , … , x t ℓ ( M ) x_{t_\ell}^{(0)},\ldots,x_{t_\ell}^{(M)} x t ℓ ( 0 ) , … , x t ℓ ( M )
end for
return { ( x ( 0 : n ) τ ( m ) , a 1 : ( L − 1 ) m , w L ( m ) ) } m = 0 M \{(x_{(0:n)\tau}^{(m)},a_{1:(L-1)}^m,w_L^{(m)})\}_{m=0}^{M} {( x ( 0 : n ) τ ( m ) , a 1 : ( L − 1 ) m , w L ( m ) ) } m = 0 M
Theorem 1, below, is proved in Appendix § 11 .§ 1 and shows that the exchangeable particle Gibbs step targets the correct conditional distribution. For simplicity in the statement and proof, we take observations to be at integer times, t ℓ = ℓ t_\ell=\ell t ℓ = ℓ ℓ = 1 , … , L \ell=1,\ldots,L ℓ = 1 , … , L
Theorem 1. For any i = 1 , … , M i=1,\ldots,M i = 1 , … , M p ( x ( 0 : n ) τ ∣ y 1 : L , θ ) p(x_{(0:n)\tau}\mid y_{1:L},\theta) p ( x ( 0 : n ) τ ∣ y 1 : L , θ )
p ( x ( 0 : n ) τ ( 0 ) ∣ y 1 : L , θ ) q ~ ( u ∣ { x ( 0 : n ) τ ( 0 ) } ∣ x ( 0 : n ) τ ( 0 ) , θ ) × q ∗ ( i ∣ u , θ ) = p ( x ( 0 : n ) τ ( i ) ∣ y 1 : L , θ ) q ~ ( u ∣ { x ( 0 : n ) τ ( i ) } ∣ x ( 0 : n ) τ ( i ) , θ ) × q ∗ ( 0 ∣ u , θ ) . \begin{aligned}
p\left(x_{(0:n)\tau}^{(0)}\mid y_{1:L},\theta\right)\widetilde{q}\left(u\mid\left\{x_{(0:n)\tau}^{(0)}\right\}\mid x_{(0:n)\tau}^{(0)},\theta\right)\times q_*(i\mid u,\theta) \\
&= p\left(x_{(0:n)\tau}^{(i)}\mid y_{1:L},\theta\right)\widetilde{q}\left(u\mid\left\{x_{(0:n)\tau}^{(i)}\right\}\mid x_{(0:n)\tau}^{(i)},\theta\right)\times q_*(0\mid u,\theta).
\end{aligned} p ( x ( 0 : n ) τ ( 0 ) ∣ y 1 : L , θ ) q ( u ∣ { x ( 0 : n ) τ ( 0 ) } ∣ x ( 0 : n ) τ ( 0 ) , θ ) × q ∗ ( i ∣ u , θ ) = p ( x ( 0 : n ) τ ( i ) ∣ y 1 : L , θ ) q ( u ∣ { x ( 0 : n ) τ ( i ) } ∣ x ( 0 : n ) τ ( i ) , θ ) × q ∗ ( 0 ∣ u , θ ) . Backward Sampling for exchangeable particle Gibbs
We now describe an approach for incorporating backward sampling into both particle Gibbs and exchangeable particle Gibbs for MJPs, yielding particle Gibbs with backward sampling (PGBS) and exchangeable particle Gibbs with backward sampling (xPGBS), respectively. This method is only practicable for reaction networks in which the reaction counts are uniquely determined by a given pair of states, such as the SIR and SEIR models.
For the backward sampling, we take the state at observation time t 1 t_1 t 1 x 1 = ( x 0 , x τ , … , x t ℓ ) T x_1=(x_0,x_\tau,\ldots,x_{t_\ell})^{\mathsf{T}} x 1 = ( x 0 , x τ , … , x t ℓ ) T t ℓ t_\ell t ℓ x ℓ = ( x t ℓ − 1 + τ , x t ℓ − 1 + 2 τ , … , x t ℓ ) T x_\ell=(x_{t_{\ell-1}+\tau},x_{t_{\ell-1}+2\tau},\ldots,x_{t_\ell})^{\mathsf{T}} x ℓ = ( x t ℓ − 1 + τ , x t ℓ − 1 + 2 τ , … , x t ℓ ) T ℓ = 2 , … , L \ell=2,\ldots,L ℓ = 2 , … , L B ~ ℓ \widetilde{B}_\ell B ℓ x ℓ x_\ell x ℓ B ~ L = m \widetilde{B}_L=m B L = m ℓ = L − 1 , … , 1 \ell=L-1,\ldots,1 ℓ = L − 1 , … , 1 b ~ ℓ m \widetilde{b}_\ell^m b ℓ m P ( B ~ ℓ m = i ) ∝ w ~ ℓ i \mathbb{P}(\widetilde{B}_\ell^m=i)\propto\widetilde{w}_\ell^i P ( B ℓ m = i ) ∝ w ℓ i
w ~ ℓ i = w ℓ i ⋅ p ( x ℓ + 1 ( b ~ ℓ + 1 ) ∣ x ℓ ( i ) , θ ) = w ℓ i ⋅ p ( x t ℓ + τ ( b ~ ℓ + 1 ) ∣ x t ℓ ( i ) , θ ) , i = 0 , … , M . \widetilde{w}_\ell^i=w_\ell^i\cdot p\left(x_{\ell+1}^{(\widetilde{b}_{\ell+1})}\mid x_\ell^{(i)},\theta\right)=w_\ell^i\cdot p\left(x_{t_{\ell}+\tau}^{(\widetilde{b}_{\ell+1})}\mid x_{t_\ell}^{(i)},\theta\right),\qquad i=0,\ldots,M. w ℓ i = w ℓ i ⋅ p ( x ℓ + 1 ( b ℓ + 1 ) ∣ x ℓ ( i ) , θ ) = w ℓ i ⋅ p ( x t ℓ + τ ( b ℓ + 1 ) ∣ x t ℓ ( i ) , θ ) , i = 0 , … , M . The validity of this approach is proved in Appendix § 12 . Under this construction, we only need to compute the transition probability between states across one time step.
Properties
Efficiency metric
For some test function h : X → R h:\mathcal{X}\to\mathbb{R} h : X → R effective sample size (ESS) (e.g. Fearnhead et al. [4 ] ) is the number n n n Z 1 , … , Z n ∼ π Z_1,\ldots,Z_n\sim\pi Z 1 , … , Z n ∼ π
Var [ 1 n ∑ i = 1 n h ( Z i ) ] = Var [ I ^ N i t e r ( h ) ] , where I ^ N i t e r ( h ) = 1 N i t e r ∑ i = 1 N i t e r h ( X i ) , \operatorname{Var}\left[\frac{1}{n}\sum_{i=1}^{n}h(Z_i)\right]=\operatorname{Var}\left[\hat{I}_{N_{\mathrm{iter}}}(h)\right],
\qquad\text{where}\qquad\hat{I}_{N_{\mathrm{iter}}}(h)=\frac{1}{N_{\mathrm{iter}}}\sum_{i=1}^{N_{\mathrm{iter}}}h(X_i), Var [ n 1 i = 1 ∑ n h ( Z i ) ] = Var [ I ^ N iter ( h ) ] , where I ^ N iter ( h ) = N iter 1 i = 1 ∑ N iter h ( X i ) , when X 0 ∼ π X_0\sim\pi X 0 ∼ π E S S = N i t e r / I A C h \mathrm{ESS}=N_{\mathrm{iter}}/\mathrm{IAC}_h ESS = N iter / IAC h I A C h = 1 + 2 ∑ k = 1 ∞ ρ k \mathrm{IAC}_h=1+2\sum_{k=1}^{\infty}\rho_k IAC h = 1 + 2 ∑ k = 1 ∞ ρ k integrated autocorrelation time and ρ k : = Corr [ h ( X 0 ) , h ( X k ) ] = Corr [ h ( X j ) , h ( X j + k ) ] \rho_k:=\operatorname{Corr}[h(X_0),h(X_k)]=\operatorname{Corr}[h(X_j),h(X_{j+k})] ρ k := Corr [ h ( X 0 ) , h ( X k )] = Corr [ h ( X j ) , h ( X j + k )] k k k
Define σ h 2 = Var ( Z ) \sigma_h^2=\operatorname{Var}(Z) σ h 2 = Var ( Z ) Z ∼ π Z\sim\pi Z ∼ π I ^ n ( h ) \hat{I}_n(h) I ^ n ( h ) lim n → ∞ n Var [ I ^ n ( h ) ] \lim_{n\to\infty}n\operatorname{Var}[\hat{I}_n(h)] lim n → ∞ n Var [ I ^ n ( h )] ∑ k = 1 ∞ ∣ ρ k ∣ < ∞ \sum_{k=1}^{\infty}|\rho_k|<\infty ∑ k = 1 ∞ ∣ ρ k ∣ < ∞ I ^ n ( h ) \hat{I}_n(h) I ^ n ( h ) σ h 2 I A C h \sigma_h^2\mathrm{IAC}_h σ h 2 IAC h [4 ] ).
Another commonly used measure of Markov chain mixing is the expected squared jump distance (ESJD) (e.g. Sherlock and Roberts [18 ] ). For a test function h h h
E S J D = E X j ∼ π [ { h ( X j + 1 ) − h ( X j ) } 2 ] . \mathrm{ESJD}=\mathbb{E}_{X_j\sim\pi}\left[\{h(X_{j+1})-h(X_j)\}^2\right]. ESJD = E X j ∼ π [ { h ( X j + 1 ) − h ( X j ) } 2 ] . Maximising ESJD minimises the lag-1 autocorrelation (e.g., Sherlock et al. [16 ] ). In the high-dimensional limit, for certain RWM and MALA algorithms (e.g. Roberts and Rosenthal [15 ] ), ESJD is the speed of the limiting diffusion for the process, so maximising ESJD maximises the speed of the diffusion, which is equivalent to minimising the integrated auto-correlation time and, hence, maximising the (theoretical) ESS.
Consider two successive reference paths generated by xPG without backward sampling, x [ 0 , T ] x_{[0,T]} x [ 0 , T ] x [ 0 , T ] ′ x'_{[0,T]} x [ 0 , T ] ′ x [ 0 , T ] ′ = x [ 0 , T ] x'_{[0,T]}=x_{[0,T]} x [ 0 , T ] ′ = x [ 0 , T ] t ∗ ∈ { t 1 , … , t L } t^*\in\{t_1,\ldots,t_L\} t ∗ ∈ { t 1 , … , t L } x [ 0 , t ∗ ] ′ = x [ 0 , t ∗ ] x'_{[0,t^*]}=x_{[0,t^*]} x [ 0 , t ∗ ] ′ = x [ 0 , t ∗ ]
If x t ′ ≠ x t x'_t\ne x_t x t ′ = x t t t t x [ t , T ] ′ x'_{[t,T]} x [ t , T ] ′ t t t x [ t , T ] x_{[t,T]} x [ t , T ] x t x_t x t x t ′ x'_t x t ′ s s s s < t s<t s < t x [ 0 , s ] ′ = x [ 0 , s ] x'_{[0,s]}=x_{[0,s]} x [ 0 , s ] ′ = x [ 0 , s ]
Remark 1. For an MJP, the state space is discrete; therefore, even though a new path might not be the reference path—in that the x T x_T x T x T ( 0 ) x_T^{(0)} x T ( 0 )
Poor mixing of certain components in the reaction network may lead to incomplete exploration of the posterior of θ \theta θ X t X_t X t
In addition, for any algorithm, we can increase the ESS by a factor of k k k k k k k k k M M M δ \delta δ δ = 1 \delta=1 δ = 1 δ \delta δ k k k k k k k k k
E f f = min ℓ , d E S S ( X t ℓ d ) M . (13) \mathrm{Eff}=\frac{\min_{\ell,d}\mathrm{ESS}(X_{t_\ell}^{d})}{M}. \tag*{(13)} Eff = M min ℓ , d ESS ( X t ℓ d ) . (13) where X t ℓ d X_{t_\ell}^{d} X t ℓ d d d d t ℓ t_\ell t ℓ
Tradeoff between smaller jumps and higher acceptance rate
For Metropolis–Hastings algorithms such as the random-walk Metropolis (RWM), the trade off between smaller jumps and higher acceptance rate is well known (e.g. Fearnhead et al. [5 ] ). For particle Gibbs and exchangeable particle Gibbs, at any given time, acceptance is equivalent to the new reference path at that time not being the old reference path. To be clear, in light of Remark 1, acceptance can include a new reference path which has the same state as the old reference path at that time point, as long as it is not the reference path itself. With this definition we find a similar trade-off: the smaller δ \delta δ
As in, Roberts et al. [14 ] , for example, in higher dimensions, a smaller jump is required. Contour plots in Section 4.1 show that, for fixed M M M δ \delta δ
Tuning particle Gibbs and exchangeable particle Gibbs
Optimal acceptance rate
We define the acceptance rate for X t X_t X t t t t α ( t ) \alpha(t) α ( t ) t t t x t ( 0 ) x_t^{(0)} x t ( 0 ) δ \delta δ M M M α ( 0 ) ≈ e − 1 ≈ 0.368 \alpha(0) \approx e^{-1} \approx0.368 α ( 0 ) ≈ e − 1 ≈ 0.368
Due to the structure of xPG and particle Gibbs, the empirical acceptance rates are non-decreasing over time; that is, α ( s ) ≤ α ( t ) \alpha(s) \le\alpha(t) α ( s ) ≤ α ( t ) s < t s < t s < t s , t ∈ { 0 , t 1 , … , t L } s,t \in\{0,t_1,\ldots,t_L\} s , t ∈ { 0 , t 1 , … , t L } t 1 > 0 t_1 > 0 t 1 > 0 α ( 0 ) = α ( t 1 ) \alpha(0) = \alpha(t_1) α ( 0 ) = α ( t 1 ) X 0 X_0 X 0 α ( 0 ) \alpha(0) α ( 0 ) α ( 0 ) \alpha(0) α ( 0 ) α \alpha α
Within a single run of the conditional exchangeable particle filter, recall that at observation time t ℓ t_\ell t ℓ x t ℓ ( 0 : M ) x_{t_\ell}^{(0:M)} x t ℓ ( 0 : M ) M M M w ℓ ( 0 : M ) w_\ell^{(0:M)} w ℓ ( 0 : M ) x ~ t ℓ ( 1 : M ) \widetilde{x}_{t_\ell}^{(1:M)} x t ℓ ( 1 : M ) x ~ t ℓ ( 0 ) = x t ℓ ( 0 ) \widetilde{x}_{t_\ell}^{(0)} = x_{t_\ell}^{(0)} x t ℓ ( 0 ) = x t ℓ ( 0 )
Let W ℓ ( 0 ) W_\ell^{(0)} W ℓ ( 0 ) t ℓ t_\ell t ℓ W ‾ ℓ \overline{W}_\ell W ℓ M M M t ℓ t_\ell t ℓ
R ℓ = W ℓ ( 0 ) W ‾ ℓ . R_\ell= \frac{W_\ell^{(0)}}{\overline{W}_\ell}. R ℓ = W ℓ W ℓ ( 0 ) . Further, for ℓ = 1 , … , L \ell= 1,\ldots,L ℓ = 1 , … , L M ℓ G M_\ell^G M ℓ G t ℓ t_\ell t ℓ t j t_j t j j < ℓ j < \ell j < ℓ
Assumption 1. For each observation time t ℓ t_\ell t ℓ ℓ = 1 , … , L \ell= 1,\ldots,L ℓ = 1 , … , L
Assumption 2. For ℓ = 2 , … , L \ell= 2,\ldots,L ℓ = 2 , … , L W ℓ ( 0 ) W_\ell^{(0)} W ℓ ( 0 ) W ‾ ℓ \overline{W}_\ell W ℓ M ℓ − 1 G M_{\ell-1}^G M ℓ − 1 G M M M
Assumption 3. { R ℓ } ℓ ∈ N \{R_\ell\}_{\ell\in\mathbb{N}} { R ℓ } ℓ ∈ N
a stationary distribution π R \pi_R π R E π R [ R 2 ] < ∞ \mathbb{E}_{\pi_R}[R^2] < \infty E π R [ R 2 ] < ∞
a finite asymptotic variance, denoted by v R v_R v R
We do not expect Assumptions 1, 2, and 3 to hold exactly, but to be reasonable approximations; see Appendix § 13 .§ 1 for a rationale and discussion.
Define p ℓ ( G ) p_\ell(G) p ℓ ( G ) t ℓ t_\ell t ℓ [ 0 , t ℓ ] [0,t_\ell] [ 0 , t ℓ ] E [ p ℓ ( G ) ] \mathbb{E}[p_\ell(G)] E [ p ℓ ( G )] M 1 G , … , M ℓ − 1 G M_1^G,\ldots,M_{\ell-1}^G M 1 G , … , M ℓ − 1 G X t 0 : t ℓ ( 0 : M ) X_{t_0:t_\ell}^{(0:M)} X t 0 : t ℓ ( 0 : M ) X 0 ( 0 ) X_0^{(0)} X 0 ( 0 ) X 0 X_0 X 0 § 13 .§ 2 .
Lemma 2. For ℓ = 1 , … , L \ell= 1,\ldots,L ℓ = 1 , … , L
P ℓ ( M ) : = M W ‾ ℓ M W ‾ ℓ + W ℓ ( 0 ) ≡ 1 1 + R ℓ / M . P_\ell^{(M)} := \frac{M\overline{W}_\ell}{M\overline{W}_\ell+ W_\ell^{(0)}} \equiv\frac{1}{1+R_\ell/M}. P ℓ ( M ) := M W ℓ + W ℓ ( 0 ) M W ℓ ≡ 1 + R ℓ / M 1 . If Assumptions 1 and 2 hold then, for each ℓ = 2 , … , L \ell= 2,\ldots,L ℓ = 2 , … , L
E M 1 G : ℓ − 1 [ p ℓ ( G ) ∣ X t 0 : t ℓ ( 0 : M ) ] = ∏ k = 1 ℓ P k ( M ) . \mathbb{E}_{M_1^G:\ell-1}\left[p_\ell(G)\mid X_{t_0:t_\ell}^{(0:M)}\right] = \prod_{k=1}^{\ell} P_k^{(M)}. E M 1 G : ℓ − 1 [ p ℓ ( G ) ∣ X t 0 : t ℓ ( 0 : M ) ] = k = 1 ∏ ℓ P k ( M ) . We are, therefore, interested in
P ( L , M ) : = ∏ ℓ = 1 L P ℓ ( M ) . (14) P^{(L,M)} := \prod_{\ell=1}^{L} P_{\ell}^{(M)}. \tag*{(14)} P ( L , M ) := ℓ = 1 ∏ L P ℓ ( M ) . (14) Defining E ℓ : = E [ R ℓ ] E_{\ell} := \mathbb{E}[R_{\ell}] E ℓ := E [ R ℓ ]
E ( L , M ) : = exp ( − 1 M ∑ ℓ = 1 L E ℓ ) (15) E^{(L,M)} := \exp\left(-\frac{1}{M}\sum_{\ell=1}^{L} E_{\ell}\right) \tag*{(15)} E ( L , M ) := exp ( − M 1 ℓ = 1 ∑ L E ℓ ) (15) is a good approximation to P ( M ) P^{(M)} P ( M ) M → ∞ M \to\infty M → ∞
Lemma 3.
E [ ∣ P ( L , M ) − E ( L , M ) ∣ ] ≤ 1 M 2 E [ ∑ ℓ = 1 L R ℓ 2 ] + 1 M Var [ ∑ ℓ = 1 L R ℓ ] 1 2 . \mathbb{E}\left[\left|P^{(L,M)}-E^{(L,M)}\right|\right]\leq\frac{1}{M^{2}}\mathbb{E}\left[\sum_{\ell=1}^{L}R_{\ell}^{2}\right]+\frac{1}{M}\operatorname{Var}\left[\sum_{\ell=1}^{L}R_{\ell}\right]^{\frac{1}{2}}. E [ P ( L , M ) − E ( L , M ) ] ≤ M 2 1 E [ ℓ = 1 ∑ L R ℓ 2 ] + M 1 Var [ ℓ = 1 ∑ L R ℓ ] 2 1 . Corollary 1. Subject to Assumption 3, if M = ⌈ c L ⌉ M=\lceil cL\rceil M = ⌈ c L ⌉
L 1 / 2 E [ ∣ P ( L , M ) − E ( L , M ) ∣ ] → v R c L^{1/2}\mathbb{E}\left[\left|P^{(L,M)}-E^{(L,M)}\right|\right]\to\frac{\sqrt{v_R}}{c} L 1/2 E [ P ( L , M ) − E ( L , M ) ] → c v R as L → ∞ L\to\infty L → ∞
Proof. Parts 1 and 2 of Assumption 3, respectively, imply that
1 L E [ ∑ ℓ = 1 L R ℓ 2 ] → E π R [ R 2 ] and 1 L Var [ ∑ ℓ = 1 L R ℓ ] → v R , \frac{1}{L}\mathbb{E}\left[\sum_{\ell=1}^{L}R_{\ell}^{2}\right]\to\mathbb{E}_{\pi_R}[R^{2}]
\quad\text{and}\quad
\frac{1}{L}\operatorname{Var}\left[\sum_{\ell=1}^{L}R_{\ell}\right]\to v_R, L 1 E [ ℓ = 1 ∑ L R ℓ 2 ] → E π R [ R 2 ] and L 1 Var [ ℓ = 1 ∑ L R ℓ ] → v R , as L → ∞ L\to\infty L → ∞
E [ ∣ P ( L , M ) − E ( L , M ) ∣ ] ≤ L M 2 E π R [ R 2 ] + L M 1 L Var [ ∑ ℓ = 1 L R ℓ ] 1 2 . \mathbb{E}\left[\left|P^{(L,M)}-E^{(L,M)}\right|\right]\leq\frac{L}{M^{2}}\mathbb{E}_{\pi_R}[R^{2}]+\frac{\sqrt{L}}{M}\frac{1}{\sqrt{L}}\operatorname{Var}\left[\sum_{\ell=1}^{L}R_{\ell}\right]^{\frac{1}{2}}. E [ P ( L , M ) − E ( L , M ) ] ≤ M 2 L E π R [ R 2 ] + M L L 1 Var [ ℓ = 1 ∑ L R ℓ ] 2 1 . So
L E [ ∣ P ( L , M ) − E ( L , M ) ∣ ] ≤ L 3 / 2 M 2 E π R [ R 2 ] + L M 1 L Var [ ∑ ℓ = 1 L R ℓ ] 1 2 . \sqrt{L}\mathbb{E}\left[\left|P^{(L,M)}-E^{(L,M)}\right|\right]\leq\frac{L^{3/2}}{M^{2}}\mathbb{E}_{\pi_R}[R^{2}]+\frac{L}{M}\frac{1}{\sqrt{L}}\operatorname{Var}\left[\sum_{\ell=1}^{L}R_{\ell}\right]^{\frac{1}{2}}. L E [ P ( L , M ) − E ( L , M ) ] ≤ M 2 L 3/2 E π R [ R 2 ] + M L L 1 Var [ ℓ = 1 ∑ L R ℓ ] 2 1 . If we set M = ⌈ c L ⌉ M=\lceil cL\rceil M = ⌈ c L ⌉ L / M → 1 / c L/M\to1/c L / M → 1/ c L → ∞ L\to\infty L → ∞ v R / c \sqrt{v_R}/c v R / c L → ∞ L\to\infty L → ∞
Thus, with M ∝ L M\propto L M ∝ L P ( L , M ) P^{(L,M)} P ( L , M ) E ( L , M ) E^{(L,M)} E ( L , M ) O ( L − 1 / 2 ) = O ( M − 1 / 2 ) O(L^{-1/2})=O(M^{-1/2}) O ( L − 1/2 ) = O ( M − 1/2 ) E ( L , M ) E^{(L,M)} E ( L , M ) E : = ∑ ℓ = 1 L E ℓ E:=\sum_{\ell=1}^{L}E_{\ell} E := ∑ ℓ = 1 L E ℓ α ≈ exp ( − E / M ) \alpha\approx\exp(-E/M) α ≈ exp ( − E / M )
Let J ( δ , M ) J(\delta,M) J ( δ , M ) X 0 X_0 X 0 X 0 X_0 X 0 X 0 ( 0 ) X_0^{(0)} X 0 ( 0 ) δ \delta δ δ \delta δ M M M M M M X 0 X_0 X 0 X 0 X_0 X 0 J ( δ , M ) J(\delta,M) J ( δ , M ) P ( M ) J ( δ , M ) ≈ exp ( − E / M ) J ( δ , M ) P^{(M)}J(\delta,M)\approx\exp(-E/M)J(\delta,M) P ( M ) J ( δ , M ) ≈ exp ( − E / M ) J ( δ , M ) § 3.4.1 , the sample efficiency is approximately equal to
1 M exp ( − E M ) J ( δ , M ) . \frac{1}{M}\exp\left(-\frac{E}{M}\right)J(\delta,M). M 1 exp ( − M E ) J ( δ , M ) . Assumption 2 indicates that E E E M M M δ \delta δ J ( δ , M ) J(\delta,M) J ( δ , M ) M M M M M M M ^ = E \hat{M}=E M ^ = E α ^ : = exp ( − 1 ) ≈ 0.368 \hat{\alpha}:=\exp(-1)\approx0.368 α ^ := exp ( − 1 ) ≈ 0.368
Although the minimum acceptance rate is attained by the first component of the latent process, the minimum ESS may, in some cases, occur at an observation time other than the first. Therefore, minimizing ESS ( X 0 ) / M \operatorname{ESS}(X_0)/M ESS ( X 0 ) / M § 3.4.1 . Nevertheless, the two quantities typically behave similarly, so the advice still provides useful guidance for achieving relatively high sample efficiency. We will delve further into this in § 4 .
Tuning strategy
From approximation of α \alpha α M M M δ \delta δ
M ^ δ = − M ∗ log ( α δ M ∗ ) . (16) \widehat{M}_{\delta}=-M^{*}\log\left(\alpha_{\delta}^{M^{*}}\right). \tag*{(16)} M δ = − M ∗ log ( α δ M ∗ ) . (16) where α δ M ∗ \alpha_{\delta}^{M^{*}} α δ M ∗ δ \delta δ M = M ∗ M=M^{*} M = M ∗ (16) provides an approximation of the optimal value of M M M α / M \alpha/M α / M δ \delta δ α / M \alpha/M α / M α / M \alpha/M α / M δ 2 \delta^{2} δ 2 δ \delta δ
In practice, to maximise efficiency over ( δ , M ) (\delta,M) ( δ , M ) α δ M ≃ α ^ \alpha_{\delta}^{M}\simeq\widehat{\alpha} α δ M ≃ α δ = 0.5 \delta=0.5 δ = 0.5 M = M ∗ M=M^{*} M = M ∗ α 0.5 M ∗ \alpha_{0.5}^{M^{*}} α 0.5 M ∗ α 0.5 M ∗ < α ^ \alpha_{0.5}^{M^{*}}<\widehat{\alpha} α 0.5 M ∗ < α δ \delta δ α δ M ∗ ≃ α ^ \alpha_{\delta}^{M^{*}}\simeq\widehat{\alpha} α δ M ∗ ≃ α α 0.5 M ∗ > α ^ \alpha_{0.5}^{M^{*}}>\widehat{\alpha} α 0.5 M ∗ > α δ \delta δ α δ M ∗ ≃ α ^ \alpha_{\delta}^{M^{*}}\simeq\widehat{\alpha} α δ M ∗ ≃ α δ \delta δ M = M ∗ M=M^{*} M = M ∗ δ \delta δ M M M α ^ \widehat{\alpha} α
Numerical Experiments
In this section, we apply xPG and particle Gibbs to synthetic data. We simulate the true latent process and generate observations based on it. We then run xPG and particle Gibbs given the observations to infer the hidden state and the model parameters. The code and data used in these experiments are available at https://github.com/yanglanya980929/xPGTauLeap.
We conduct the experiments on the cases where the latent process is governed by the following classes of models: (i) The product of multiple independent SIS models, with noisy observations of the number of infectious people in each region. (ii) Multiple coupled SIR models, with noisy observations of the number of infectious people in each region. (iii) Multiple coupled SEIR models, with noisy observations of the number of infectious people in each region. (iv) A single autoregulatory model, with exact observations of the number of molecules of RNA, P and P2 _2 2
In the first three cases, each region has population N = 1000 N=1000 N = 1000
When reporting results we use the following abbreviations: PG for particle Gibbs, xPG for exchangeable particle Gibbs, PGBS for particle Gibbs with backwards sampling and xPGBS for exchangeable particle Gibbs with backwards sampling.
SIS
We run xPG in multi-dimensional space, where the state in each dimension follows an independent SIS model. All the SIS models share the same true values of β = 1.2 \beta=1.2 β = 1.2 γ = 0.6 \gamma=0.6 γ = 0.6 τ = 0.2 \tau=0.2 τ = 0.2 t = 0 , 1 , … , 20 t=0,1,\ldots,20 t = 0 , 1 , … , 20 d = 1 , 3 , 5 d=1,3,5 d = 1 , 3 , 5 7 7 7 t ∈ { 0 , 1 , … , 20 } t\in\{0,1,\ldots,20\} t ∈ { 0 , 1 , … , 20 } I t i I_t^i I t i i i i i ∈ { 1 , … , d } i\in\{1,\ldots,d\} i ∈ { 1 , … , d } Y t i Y_t^i Y t i g ( Y t i ∣ X t i ) g(Y_t^i\mid X_t^i) g ( Y t i ∣ X t i ) Y t i ∈ Z Y_t^i\in\mathbb{Z} Y t i ∈ Z
g ( Y t i ∣ I t i ) = { ϕ ( ∣ Y t i − I t i ∣ ) ∑ k = − 2 σ 2 σ ϕ ( k ) , if ∣ Y t i − I t i ∣ ≤ 2 σ , 0 , otherwise g(Y_t^i\mid I_t^i)=
\begin{cases}
\dfrac{\phi\left(\lvert Y_t^i-I_t^i\rvert\right)}{\sum_{k=-2\sigma}^{2\sigma}\phi(k)}, & \text{if } \lvert Y_t^i-I_t^i\rvert\leq2\sigma,\\
0, & \text{otherwise}
\end{cases} g ( Y t i ∣ I t i ) = ⎩ ⎨ ⎧ ∑ k = − 2 σ 2 σ ϕ ( k ) ϕ ( ∣ Y t i − I t i ∣ ) , 0 , if ∣ Y t i − I t i ∣ ≤ 2 σ , otherwise where ϕ \phi ϕ σ \sigma σ N ( 0 , σ 2 ) N(0,\sigma^2) N ( 0 , σ 2 ) σ \sigma σ Y t = ( Y t 1 , … , Y t d ) Y_t=(Y_t^1,\ldots,Y_t^d) Y t = ( Y t 1 , … , Y t d ) d > 1 d>1 d > 1
I 0 i ∼ Bin ( N , 0.05 ) , i = 1 , … , d . I_0^i\sim\operatorname{Bin}(N,0.05), \qquad i=1,\ldots,d. I 0 i ∼ Bin ( N , 0.05 ) , i = 1 , … , d . For each region i i i
In both data simulation and inference, we set σ \sigma σ
β ∼ Gamma ( 1 , 1 ) and γ ∼ Gamma ( 1 , 1 ) . (17) \beta\sim\operatorname{Gamma}(1,1)\quad\text{and}\quad\gamma\sim\operatorname{Gamma}(1,1). \tag*{(17)} β ∼ Gamma ( 1 , 1 ) and γ ∼ Gamma ( 1 , 1 ) . (17) Their initial values are β 0 = 1 \beta^0 = 1 β 0 = 1 γ 0 = 1 \gamma^0 = 1 γ 0 = 1
We ran xPG and particle Gibbs (equivalent to xPG with δ = 1 \delta= 1 δ = 1 M ∈ { 1 , 5 , 10 , 30 , 50 , 70 , 100 , 120 , 150 , 180 } ∪ { 200 + 50 k : k = 0 , 1 , … , 9 } M \in\{1, 5, 10, 30, 50, 70, 100, 120, 150, 180\} \cup\{200 + 50k : k = 0, 1, \ldots, 9\} M ∈ { 1 , 5 , 10 , 30 , 50 , 70 , 100 , 120 , 150 , 180 } ∪ { 200 + 50 k : k = 0 , 1 , … , 9 } δ ∈ { 10 − 2.5 + 0.2 k : k = 0 , 1 , … , 12 } ∪ { 1 } \delta\in\{10^{-2.5+0.2k} : k = 0, 1, \ldots, 12\} \cup\{1\} δ ∈ { 1 0 − 2.5 + 0.2 k : k = 0 , 1 , … , 12 } ∪ { 1 } S + I = N S + I = N S + I = N S S S I I I
Figure 1 shows the contour plots of the sample efficiency of xPG after running 10 5 10^5 1 0 5 δ \delta δ ( δ , M ) = ( 0.79 , 50 ) (\delta, M) = (0.79, 50) ( δ , M ) = ( 0.79 , 50 ) ( 0.20 , 130 ) (0.20, 130) ( 0.20 , 130 ) ( 0.05 , 150 ) (0.05, 150) ( 0.05 , 150 ) ( 0.02 , 130 ) (0.02, 130) ( 0.02 , 130 )
Figure 1. Contour plots of sample efficiency after running xPG for 10 5 10^5 1 0 5 β = 1.2 \beta= 1.2 β = 1.2 γ = 0.6 \gamma= 0.6 γ = 0.6 t 1 t_1 t 1 exp ( − 1 ) \exp(-1) exp ( − 1 )
In addition, the plots suggest the presence of a tipping point in the relative efficiency of xPG vs particle Gibbs at around three dimensions. Beyond this point, xPG shifts from being only slightly more efficient than particle Gibbs to being substantially better, and further increasing M M M M M M M = 650 M = 650 M = 650
Moreover, as shown in Figure 1 , the red dashed line, which corresponds to the optimal acceptance rate derived in § 3.5 , generally passes through the circle indicating the highest sample efficiency for each fixed value of δ \delta δ M M M δ \delta δ
SIR
We apply xPG to d = 4 d = 4 d = 4 β = 0.6 \beta= 0.6 β = 0.6 γ = 0.2 \gamma= 0.2 γ = 0.2 d × d d \times d d × d
The length of a time step, τ \tau τ 0.2 0.2 0.2 Y t i Y_t^i Y t i t = 1 , … , 20 t = 1, \ldots, 20 t = 1 , … , 20 i i i i ∈ { 1 , 2 , 3 , 4 } i \in\{1, 2, 3, 4\} i ∈ { 1 , 2 , 3 , 4 } I t i I_t^i I t i
Y t i ∼ Bin ( I t i , 0.9 ) . (18) Y_t^i \sim\operatorname{Bin}\left(I_t^i, 0.9\right). \tag*{(18)} Y t i ∼ Bin ( I t i , 0.9 ) . (18) For each region, the initial state at time zero is assumed to be fixed and known. The priors and initial values of β \beta β γ \gamma γ (17) ).
We run PG, xPG, PGBS (equivalent to xPGBS with δ = 1 \delta= 1 δ = 1 M ∈ { 1 , 3 , 5 , 10 , 30 , 50 , 70 , 100 , 130 , 150 , 180 , ( 4 : 10 ) × 50 } M \in\{1, 3, 5, 10, 30, 50, 70, 100, 130, 150, 180, (4 : 10) \times50\} M ∈ { 1 , 3 , 5 , 10 , 30 , 50 , 70 , 100 , 130 , 150 , 180 , ( 4 : 10 ) × 50 } δ ∈ { 0.01 , 0.05 , ( 1 : 5 ) × 0.1 , 1 } \delta\in\{0.01, 0.05, (1 : 5) \times0.1, 1\} δ ∈ { 0.01 , 0.05 , ( 1 : 5 ) × 0.1 , 1 } Figure 2 shows the contour plots of the sample efficiency in terms of susceptible and infectious populations of xPG and xPGBS after running 10 6 10^6 1 0 6 ( δ ^ , M ^ ) (\hat{\delta}, \hat{M}) ( δ ^ , M ^ )
Figure 2. Contour plots of sample efficiency in terms of the number of susceptible and infectious individuals after running xPG and xPGBS for 10 6 10^6 1 0 6 β = 0.6 \beta= 0.6 β = 0.6 γ = 0.2 \gamma= 0.2 γ = 0.2
Table 1 presents the optimal parameter values of the four algorithms and the corresponding sample efficiency for the susceptible state, infectious state, and model parameters β \beta β γ \gamma γ
Table 1. Efficiency of different particle Gibbs algorithms run on simulated data from four coupled SIR models. All runs were for 10 6 10^6 1 0 6 M = M ^ M = \hat{M} M = M ^ δ = δ ^ \delta= \hat{\delta} δ = δ ^
efficiency for the worst mixing component of the state for PG being over eight times higher for xPG—and xPGBS similarly provides a marked improvement over PGBS—the sample efficiency for the worst mixing component for PGBS being approximately seven times higher for xPGBS. Relative to PG, xPG improves sample efficiency more than PGBS.
SEIR
In this section, we consider four coupled SEIR models representing four different regions. We set τ = 0.1 \tau= 0.1 τ = 0.1 § 4.2 , the data are distributed according to (18) at times t = 1 , … , 20 t = 1,\ldots,20 t = 1 , … , 20 § 6 .
The underlying true values of the model parameters are β = 0.6 \beta= 0.6 β = 0.6 η = 0.6 \eta= 0.6 η = 0.6 γ = 0.2 \gamma= 0.2 γ = 0.2 β , η ∼ Gamma ( 1 , 1 ) \beta,\eta\sim\operatorname{Gamma}(1,1) β , η ∼ Gamma ( 1 , 1 ) γ ∼ Gamma ( 20 , 100 ) \gamma\sim\operatorname{Gamma}(20,100) γ ∼ Gamma ( 20 , 100 )
Figures 3a and 3b present the results of applying xPG and xPGBS to the simulated data, respectively. Both plots show the sample efficiency for the state with the lowest ESS across all states. Corresponding plots for the susceptible, exposed and infectious states are provided in Appendix § 17 .
The mixing efficiency of the susceptible state is much lower than that obtained in § 4.2 . In the SIR model, the number of individuals in the susceptible state is relatively tightly constrained given observations of the number of infected individuals. In contrast, in the SEIR model, the inclusion of the exposed state introduces additional flexibility, so there are many more possible combinations of the numbers of susceptible and exposed individuals that are consistent with the observed infectious counts. As a result, the high-probability region of the posterior distribution of susceptible states is wider than that in the SIR model. At a given observation time point, the tight conditional distribution of the observed number infectious (18) requires a small value of δ \delta δ δ \delta δ
As in the previous section, we define the pair of tuning parameters that give the highest minimum sample efficiency across all states as the optimal parameters, ( δ ^ , M ^ ) (\hat{\delta},\hat{M}) ( δ ^ , M ^ ) < < < < < < < < < § 4.2 ). In particular, PGBS is only slightly better for the lowest sample efficiency of the model parameters for PG (η \eta η
In the results for the SIR model (§ 4.2 ) and the SEIR model, the optimal values of M M M δ \delta δ δ \delta δ M M M
For the SIR and SEIR models, the minimum ESS over all latent states and observation time points does not necessarily occur at the first observation time. As a result, the tuning strategy may overestimate the number of particles required. Suppose the minimum ESS occurs at time t ℓ ^ t_{\hat{\ell}} t ℓ ^ ℓ ^ > 1 \hat{\ell} > 1 ℓ ^ > 1 α ( t ℓ ^ ) \alpha(t_{\hat{\ell}}) α ( t ℓ ^ ) α ^ \hat{\alpha} α ^ § 3.5.1 . However, since t ℓ ^ t_{\hat{\ell}} t ℓ ^ t 1 t_1 t 1 α ^ \hat{\alpha} α ^ α ( t ℓ ^ ) ≥ α ( t 1 ) ≈ α ^ \alpha(t_{\hat{\ell}}) \ge\alpha(t_1) \approx\hat{\alpha} α ( t ℓ ^ ) ≥ α ( t 1 ) ≈ α ^ δ \delta δ α ( t 1 ) \alpha(t_1) α ( t 1 ) α ^ \hat{\alpha} α ^ M M M § 16 .
Table 2. Efficiency of different particle Gibbs algorithms run on simulated data from four coupled SEIR models. All runs were for 10 6 10^6 1 0 6
Choosing an M M M M M M M M M Figure 1 , Figure 2 a, Figure 2 c and Figure 3 a all show that, for any fixed value of δ ∈ ( 0 , 1 ] \delta\in(0, 1] δ ∈ ( 0 , 1 ] M M M M M M [19 ] and Doucet et al. [3 ] , where the mixing efficiency is shown to decrease rapidly as σ \sigma σ σ ∝ 1 / M \sigma\propto1/M σ ∝ 1/ M
Figure 3. Contour plots of sample efficiency after running xPG and xPGBS for 10 6 10^6 1 0 6 β = 0.6 \beta= 0.6 β = 0.6 η = 0.6 \eta= 0.6 η = 0.6 γ = 0.2 \gamma= 0.2 γ = 0.2
Autoregulatory model
In this section, xPG is applied to the autoregulatory model with five species and eight types of reaction introduced in § 2.2 . We simulate the data using ( θ 1 , … , θ 8 ) = ( 0.2 , 1.4 , 0.7 , 0.4 , 0.2 , 1.8 , 0.6 , 0.2 ) (\theta_1,\ldots,\theta_8) = (0.2, 1.4, 0.7, 0.4, 0.2, 1.8, 0.6, 0.2) ( θ 1 , … , θ 8 ) = ( 0.2 , 1.4 , 0.7 , 0.4 , 0.2 , 1.8 , 0.6 , 0.2 ) ( X 1 , X 2 , X 3 ) = ( 5 , 8 , 8 ) (X_1, X_2, X_3) = (5, 8, 8) ( X 1 , X 2 , X 3 ) = ( 5 , 8 , 8 ) G = 10 G = 10 G = 10 X 4 X_4 X 4 X 5 X_5 X 5 X 5 = G − X 4 X_5=G-X_4 X 5 = G − X 4 { 0 , 1 , … , G } \{0,1,\ldots,G\} { 0 , 1 , … , G } θ 1 , θ 5 ∼ Gamma ( 100 , 500 ) \theta_1,\theta_5\sim\operatorname{Gamma}(100,500) θ 1 , θ 5 ∼ Gamma ( 100 , 500 ) θ i ∼ Gamma ( 1 , 1 ) \theta_i\sim\operatorname{Gamma}(1,1) θ i ∼ Gamma ( 1 , 1 ) i ∈ { 2 , 3 , 4 , 6 , 7 , 8 } i\in\{2,3,4,6,7,8\} i ∈ { 2 , 3 , 4 , 6 , 7 , 8 } X 1 X_1 X 1 X 2 X_2 X 2 X 3 X_3 X 3 t = 0 , 1 , … , 20 t=0,1,\ldots,20 t = 0 , 1 , … , 20 M ∈ { 200 , 500 , 1000 } M\in\{200,500,1000\} M ∈ { 200 , 500 , 1000 } δ ∈ { 10 − 0.2 k , k = 0 , 1 , … , 12 } \delta\in\{10^{-0.2k},k=0,1,\ldots,12\} δ ∈ { 1 0 − 0.2 k , k = 0 , 1 , … , 12 } Figure 4 shows the sample efficiency for X 4 X_4 X 4 X 5 X_5 X 5 log 10 δ \log_{10}\delta log 10 δ δ \delta δ M M M
Figure 4. Sample efficiency of xPG run for 10 5 10^{5} 1 0 5 log 10 δ \log_{10}\delta log 10 δ M M M 2 _2 2 θ 1 : 8 = ( 0.2 , 1.4 , 0.7 , 0.4 , 0.2 , 1.8 , 0.6 , 0.2 ) \theta_{1:8}=(0.2,1.4,0.7,0.4,0.2,1.8,0.6,0.2) θ 1 : 8 = ( 0.2 , 1.4 , 0.7 , 0.4 , 0.2 , 1.8 , 0.6 , 0.2 ) X 1 : 5 = ( 5 , 8 , 8 , 8 , 2 ) X_{1:5}=(5,8,8,8,2) X 1 : 5 = ( 5 , 8 , 8 , 8 , 2 )
At each observation time, only proposals which hit the three-dimensional observation receive non-zero weight and can therefore survive resampling. The standard PG is unlikely to generate such latent states. Consequently, most particles receive zero weight, leading to severe particle degeneracy in PG. In contrast, xPG mitigates this issue by proposing paths that are highly correlated with the reference path. In the simulation study, when M = 1000 M=1000 M = 1000 X 4 X_4 X 4 δ = 0.05 \delta=0.05 δ = 0.05 δ = 0.5 \delta=0.5 δ = 0.5 X 4 X_4 X 4 δ = 1 \delta=1 δ = 1 M M M
Discussion
We have introduced xPG and xPGBS methods for estimating the model parameters and latent states of hidden Markov models whose latent processes are Markov jump processes. Our xPG algorithm builds on the xPG method for diffusions introduced in Malory [13 ] and Finke and Thiery [7 ] , but focuses on Markov jump processes discretised by the popular tau-leap method. On a variety of reaction-networks, we demonstrate that it mitigates the particle degeneracy and mixing issues for which it is designed. The xPGBS algorithm adds an additional backward sampling step to further improve the sample efficiency of both the latent states and the model parameters. However, this approach is only applicable when the number of reactions of each type can be uniquely identified from two states at consecutive time points, as in the SIR and SEIR epidemic models.
In addition, we have suggested a strategy for tuning the two parameters: δ \delta δ δ ⇒ \delta\Rightarrow δ ⇒ M M M M M M Figure 1 ) show that picking δ \delta δ M M M
We applied xPG to simulated data from models composed of multiple independent SIS sub-models; since SIS is one-dimensional, the state dimension is the number of models. In each dimension, inference was performed for each point in a grid of ( δ , M ) (\delta,M) ( δ , M ) δ \delta δ M = 650 M=650 M = 650
We then compared the performance of particle Gibbs and xPG, with and without backwards sampling, on data simulated from four-region coupled SIR and SEIR models. These experiments further illustrate the efficiency gains achieved by our methods. For both the SIR and SEIR models, adding backward sampling to particle Gibbs improves sample efficiency, but much less than using xPG. xPGBS achieves the highest sample efficiency, suggesting that this should be the method of choice when backward sampling is possible.
The benefit of xPG over PG is most pronounced when the observations are highly informative. In such settings, xPG can explore regions of high likelihood more effectively by constraining the proposed particles to remain close to the reference trajectory. In contrast, when the observations are less informative, standard particle Gibbs may already produce a reasonably large proportion of particles with high weights.
To the best of our knowledge, xPG is also the first algorithm to make exact joint inference of the latent states and rate parameters of the five-dimensional autoregulatory model of Golightly and Wilkinson [10 ] from a sequence of exact state observations feasible. This setting is particularly challenging for many particle MCMC methods because the low observation noise can lead to severe particle degeneracy. By contrast, alternative Bayesian inference approaches, such as the nearly minimal extended state-space algorithm (nMESA) [17 ] , are designed primarily for inference on the rate parameters rather than the latent state trajectories. Running xPG with M = 1000 M = 1000 M = 1000 10 5 10^{5} 1 0 5 δ = 0.05 \delta= 0.05 δ = 0.05 [17 ] report that running nMESA for 2 × 10 5 2 \times10^{5} 2 × 1 0 5
Our method is based on the tau-leap approximation of MJPs. Smaller τ \tau τ
Further details on reaction networks
The stoichiometry matrix for the SIS model is
Γ = ( − 1 1 1 − 1 ) . \Gamma= \begin{pmatrix} -1 & 1 \\ 1 & -1 \end{pmatrix}. Γ = ( − 1 1 1 − 1 ) . The stoichiometry matrix for the SIR model is
Γ = ( − 1 0 1 − 1 0 1 ) . \Gamma= \begin{pmatrix} -1 & 0 \\ 1 & -1 \\ 0 & 1 \end{pmatrix}. Γ = − 1 1 0 0 − 1 1 . When considering four coupled SIR models, we use a connectivity matrix of
C = ( 1 0.1 0.1 0.1 0.2 1 0.05 0.1 0.1 0.15 1 0.1 0.1 0.2 0.1 1 ) . C = \begin{pmatrix} 1 & 0.1 & 0.1 & 0.1 \\ 0.2 & 1 & 0.05 & 0.1 \\ 0.1 & 0.15 & 1 & 0.1 \\ 0.1 & 0.2 & 0.1 & 1 \end{pmatrix}. C = 1 0.2 0.1 0.1 0.1 1 0.15 0.2 0.1 0.05 1 0.1 0.1 0.1 0.1 1 . The stoichiometry matrix for the SEIR model is
Γ = ( − 1 0 0 1 − 1 0 0 1 − 1 0 0 1 ) . \Gamma= \begin{pmatrix} -1 & 0 & 0 \\ 1 & -1 & 0 \\ 0 & 1 & -1 \\ 0 & 0 & 1 \end{pmatrix}. Γ = − 1 1 0 0 0 − 1 1 0 0 0 − 1 1 . When considering four coupled SEIR models, we use a connectivity matrix of
C = ( 1 0.2 0.1 0.1 0.2 1 0.15 0.2 0.1 0.15 1 0.1 0.1 0.2 0.1 1 ) . C = \begin{pmatrix} 1 & 0.2 & 0.1 & 0.1 \\ 0.2 & 1 & 0.15 & 0.2 \\ 0.1 & 0.15 & 1 & 0.1 \\ 0.1 & 0.2 & 0.1 & 1 \end{pmatrix}. C = 1 0.2 0.1 0.1 0.2 1 0.15 0.2 0.1 0.15 1 0.1 0.1 0.2 0.1 1 . The reactions in the autoregulatory model are given by
D N A + P 2 ⇌ θ 1 ( G − X 4 ) X 3 θ 2 X 4 D N A ⋅ P 2 , D N A → θ 3 ( G − X 4 ) D N A + R N A , R N A → θ 4 X 1 R N A + P , 2 P ⇌ θ 5 X 2 ( X 2 − 1 ) / 2 θ 6 X 3 P 2 , R N A → θ 7 X 1 ∅ , P → θ 8 X 2 ∅ . \begin{aligned}
\mathrm{DNA} + \mathrm{P}_2 &\underset{\theta_2 X_4}{\overset{\theta_1(G-X_4)X_3}{\rightleftharpoons}} \mathrm{DNA}\cdot\mathrm{P}_2, \\
\mathrm{DNA} &\xrightarrow{\theta_3(G-X_4)} \mathrm{DNA} + \mathrm{RNA}, \\
\mathrm{RNA} &\xrightarrow{\theta_4 X_1} \mathrm{RNA} + \mathrm{P}, \\
2\mathrm{P} &\underset{\theta_6 X_3}{\overset{\theta_5 X_2(X_2-1)/2}{\rightleftharpoons}} \mathrm{P}_2, \\
\mathrm{RNA} &\xrightarrow{\theta_7 X_1} \varnothing, \\
\mathrm{P} &\xrightarrow{\theta_8 X_2} \varnothing.
\end{aligned} DNA + P 2 DNA RNA 2 P RNA P θ 2 X 4 ⇌ θ 1 ( G − X 4 ) X 3 DNA ⋅ P 2 , θ 3 ( G − X 4 ) DNA + RNA , θ 4 X 1 RNA + P , θ 6 X 3 ⇌ θ 5 X 2 ( X 2 − 1 ) /2 P 2 , θ 7 X 1 ∅ , θ 8 X 2 ∅ . Here, X 1 : 5 X_{1:5} X 1 : 5 2 _2 2 ⋅ \cdot ⋅ 2 _2 2 X 4 + X 5 = G X_4 + X_5 = G X 4 + X 5 = G
Conditional Particle filter
Algorithm 4 Conditional Particle Filter
Require: the parameter, θ \theta θ x 0 : n τ ( 0 ) x_{0:n\tau}^{(0)} x 0 : n τ ( 0 ) y t 1 : L y_{t_{1:L}} y t 1 : L q 0 ( ⋅ ∣ θ ) q_0(\cdot\mid\theta) q 0 ( ⋅ ∣ θ ) X 0 X_0 X 0 q τ ( ⋅ ∣ x t − 1 , θ ) q_\tau(\cdot\mid x_{t-1},\theta) q τ ( ⋅ ∣ x t − 1 , θ ) X t X_t X t X t − 1 = x t − 1 X_{t-1}=x_{t-1} X t − 1 = x t − 1
Sample x 0 ( 1 : M ) ∼ i . i . d . q 0 ( x 0 ∣ θ ) x_0^{(1:M)} \mathrel{\overset{\mathrm{i.i.d.}}{\sim}} q_0(x_0\mid\theta) x 0 ( 1 : M ) ∼ i.i.d. q 0 ( x 0 ∣ θ )
For k = 1 k=1 k = 1 K K K x k τ ( m ) ∼ q τ ( ⋅ ∣ x ( k − 1 ) τ ( m ) , θ ) x_{k\tau}^{(m)}\sim q_\tau(\cdot\mid x_{(k-1)\tau}^{(m)},\theta) x k τ ( m ) ∼ q τ ( ⋅ ∣ x ( k − 1 ) τ ( m ) , θ ) m = 1 , … , M m=1,\ldots,M m = 1 , … , M
Compute the weights of x t 1 ( 0 ) , … , x t 1 ( M ) x_{t_1}^{(0)},\ldots,x_{t_1}^{(M)} x t 1 ( 0 ) , … , x t 1 ( M ) x ~ 0 ( m ) = x 0 ( m ) \tilde{x}_0^{(m)}=x_0^{(m)} x ~ 0 ( m ) = x 0 ( m ) m = 0 , … , M m=0,\ldots,M m = 0 , … , M
for ℓ = 2 \ell=2 ℓ = 2 L L L do
Set x ~ t ℓ − 1 ( 0 ) = x t ℓ − 1 ( 0 ) \tilde{x}_{t_{\ell-1}}^{(0)}=x_{t_{\ell-1}}^{(0)} x ~ t ℓ − 1 ( 0 ) = x t ℓ − 1 ( 0 )
Draw M M M { x t ℓ − 1 ( m ) } m = 0 M \{x_{t_{\ell-1}}^{(m)}\}_{m=0}^{M} { x t ℓ − 1 ( m ) } m = 0 M x t ℓ − 1 ( m ) x_{t_{\ell-1}}^{(m)} x t ℓ − 1 ( m ) w ℓ − 1 m w_{\ell-1}^{m} w ℓ − 1 m { x ~ t ℓ − 1 ( m ) } m = 1 M \{\tilde{x}_{t_{\ell-1}}^{(m)}\}_{m=1}^{M} { x ~ t ℓ − 1 ( m ) } m = 1 M a ℓ − 1 1 : M a_{\ell-1}^{1:M} a ℓ − 1 1 : M
for k = 1 k=1 k = 1 K K K do
Sample x t ℓ − 1 + k τ ( m ) ∼ q τ ( ⋅ ∣ x ~ t ℓ − 1 + ( k − 1 ) τ ( m ) , θ ) x_{t_{\ell-1}+k\tau}^{(m)}\sim q_\tau(\cdot\mid\tilde{x}_{t_{\ell-1}+(k-1)\tau}^{(m)},\theta) x t ℓ − 1 + k τ ( m ) ∼ q τ ( ⋅ ∣ x ~ t ℓ − 1 + ( k − 1 ) τ ( m ) , θ ) m = 1 , … , M m=1,\ldots,M m = 1 , … , M
If k < K k<K k < K x ~ t ℓ − 1 + k τ ( m ) = x t ℓ − 1 + k τ ( m ) \tilde{x}_{t_{\ell-1}+k\tau}^{(m)}=x_{t_{\ell-1}+k\tau}^{(m)} x ~ t ℓ − 1 + k τ ( m ) = x t ℓ − 1 + k τ ( m ) m = 1 , … , M m=1,\ldots,M m = 1 , … , M
end for
Compute the weights of x t ℓ ( 0 ) , … , x t ℓ ( M ) x_{t_\ell}^{(0)},\ldots,x_{t_\ell}^{(M)} x t ℓ ( 0 ) , … , x t ℓ ( M ) q ( ⋅ ∣ x t ℓ − 1 , θ ) = ∏ k = 1 K q τ ( x t ℓ − 1 + k τ ∣ x t ℓ − 1 + ( k − 1 ) τ , θ ) q(\cdot\mid x_{t_{\ell-1}},\theta)=\prod_{k=1}^{K}q_\tau(x_{t_{\ell-1}+k\tau}\mid x_{t_{\ell-1}+(k-1)\tau},\theta) q ( ⋅ ∣ x t ℓ − 1 , θ ) = ∏ k = 1 K q τ ( x t ℓ − 1 + k τ ∣ x t ℓ − 1 + ( k − 1 ) τ , θ )
end for
return { ( x 0 : n τ ( m ) , a 1 : ( L − 1 ) m , w L m ) } m = 0 M \{(x_{0:n\tau}^{(m)},a_{1:(L-1)}^{m},w_L^m)\}_{m=0}^{M} {( x 0 : n τ ( m ) , a 1 : ( L − 1 ) m , w L m ) } m = 0 M
Particle degeneracy
Exchangeable particle Gibbs in continuous state space
In practice, Algorithm 5 illustrates how we generate x t ( 1 : M ) x_t^{(1:M)} x t ( 1 : M ) x t ( 0 ) x_t^{(0)} x t ( 0 ) x t − τ ( 0 : M ) x_{t-\tau}^{(0:M)} x t − τ ( 0 : M ) d d d z ( 0 ) z^{(0)} z ( 0 ) t − τ t-\tau t − τ t t t (3) , we construct an invertible mapping h : R d × R d → R d h:\mathbb{R}^d\times\mathbb{R}^d\to\mathbb{R}^d h : R d × R d → R d X = h ( Z , x t − τ ) X=h(Z,x_{t-\tau}) X = h ( Z , x t − τ ) q ( ⋅ ∣ x t − τ , θ ) q(\cdot\mid x_{t-\tau},\theta) q ( ⋅ ∣ x t − τ , θ )
h ( Z , x ) = x + μ ( x ) τ + σ ( x ) τ Z , h(Z,x)=x+\mu(x)\tau+\sqrt{\sigma(x)\tau}Z, h ( Z , x ) = x + μ ( x ) τ + σ ( x ) τ Z , therefore,
Z ( 0 ) = h − 1 ( X t ; x t − τ ) ∼ N ( 0 , I d ) Z^{(0)} = h^{-1}(X_t; x_{t-\tau}) \sim N(0, I_d) Z ( 0 ) = h − 1 ( X t ; x t − τ ) ∼ N ( 0 , I d ) We then sample Z ^ ( 0 ) , … , Z ^ ( M ) \hat{Z}^{(0)}, \ldots, \hat{Z}^{(M)} Z ^ ( 0 ) , … , Z ^ ( M ) N d ( 0 , ϵ 2 I d ) N_d(0, \epsilon^2 I_d) N d ( 0 , ϵ 2 I d ) Z ∗ = 1 − ϵ 2 Z ( 0 ) + Z ^ ( 0 ) Z^* = \sqrt{1-\epsilon^2}Z^{(0)} + \hat{Z}^{(0)} Z ∗ = 1 − ϵ 2 Z ( 0 ) + Z ^ ( 0 ) Z ( m ) = 1 − ϵ 2 Z ∗ + Z ^ ( m ) Z^{(m)} = \sqrt{1-\epsilon^2}Z^* + \hat{Z}^{(m)} Z ( m ) = 1 − ϵ 2 Z ∗ + Z ^ ( m ) m = 1 , … , M m = 1, \ldots, M m = 1 , … , M x t ( m ) x_t^{(m)} x t ( m ) x t ( m ) = h ( z ( m ) , x t − τ ) x_t^{(m)} = h(z^{(m)}, x_{t-\tau}) x t ( m ) = h ( z ( m ) , x t − τ ) Cov ( Z ( i ) , Z ( j ) ) = 1 − ϵ 2 \operatorname{Cov}(Z^{(i)}, Z^{(j)}) = 1-\epsilon^2 Cov ( Z ( i ) , Z ( j ) ) = 1 − ϵ 2 i , j ∈ { 0 , … , M } i,j \in\{0,\ldots,M\} i , j ∈ { 0 , … , M } i ≠ j i \ne j i = j ϵ \epsilon ϵ 1 1 1 ϵ \epsilon ϵ 0 0 0
In addition, given x 0 ( 0 ) x_0^{(0)} x 0 ( 0 ) q 0 ( x 0 ∣ θ ) q_0(x_0 \mid\theta) q 0 ( x 0 ∣ θ ) M M M q ~ 0 ( x 0 ( 1 : M ) ∣ x 0 ( 0 ) , θ ) \widetilde{q}_0(x_0^{(1:M)} \mid x_0^{(0)},\theta) q 0 ( x 0 ( 1 : M ) ∣ x 0 ( 0 ) , θ ) h : R d → R d h:\mathbb{R}^d \to\mathbb{R}^d h : R d → R d Z ∼ N d ( 0 , I d ) Z \sim N_d(0,I_d) Z ∼ N d ( 0 , I d ) X = h ( Z ) X=h(Z) X = h ( Z ) q 0 ( ⋅ ∣ θ ) q_0(\cdot\mid\theta) q 0 ( ⋅ ∣ θ ) z 0 ( 0 ) = h − 1 ( x 0 ( 0 ) ) z_0^{(0)} = h^{-1}(x_0^{(0)}) z 0 ( 0 ) = h − 1 ( x 0 ( 0 ) ) x 0 ( 0 ) , … , x 0 ( M ) x_0^{(0)},\ldots,x_0^{(M)} x 0 ( 0 ) , … , x 0 ( M )
Algorithm 5 Exchangeable Proposal q ~ ( ⋅ ∣ x t ( 0 : M ) , x t − τ ( 0 ) , θ ) \widetilde{q}(\cdot\mid x_t^{(0:M)},x_{t-\tau}^{(0)},\theta) q ( ⋅ ∣ x t ( 0 : M ) , x t − τ ( 0 ) , θ ) q q q
Require: An invertible mapping h : R d × R d → R d h:\mathbb{R}^d \times\mathbb{R}^d \to\mathbb{R}^d h : R d × R d → R d h ( Z 1 : d ; x t − τ ) ∼ q ( ⋅ ∣ x t − τ , θ ) h(Z_{1:d};x_{t-\tau}) \sim q(\cdot\mid x_{t-\tau},\theta) h ( Z 1 : d ; x t − τ ) ∼ q ( ⋅ ∣ x t − τ , θ ) Z 1 : d ∼ N ( 0 , I d ) Z_{1:d} \sim N(0,I_d) Z 1 : d ∼ N ( 0 , I d ) ϵ \epsilon ϵ x t ( 0 ) x_t^{(0)} x t ( 0 )
Calculate z ( 0 ) = h − 1 ( x t ( 0 ) , x t − τ ) z^{(0)} = h^{-1}(x_t^{(0)},x_{t-\tau}) z ( 0 ) = h − 1 ( x t ( 0 ) , x t − τ )
Sample z ^ ( 0 ) \hat{z}^{(0)} z ^ ( 0 ) N ( 0 , ϵ 2 I d ) N(0,\epsilon^2 I_d) N ( 0 , ϵ 2 I d )
z ∗ = 1 − ϵ 2 z ( 0 ) + z ^ ( 0 ) z^* = \sqrt{1-\epsilon^2}z^{(0)} + \hat{z}^{(0)} z ∗ = 1 − ϵ 2 z ( 0 ) + z ^ ( 0 )
for m = 1 , … , M m = 1,\ldots,M m = 1 , … , M do
Sample z ^ ( m ) \hat{z}^{(m)} z ^ ( m ) N ( 0 , ϵ 2 I d ) N(0,\epsilon^2 I_d) N ( 0 , ϵ 2 I d )
Set z ( m ) = z ∗ 1 − ϵ 2 + z ^ ( m ) z^{(m)} = z^*\sqrt{1-\epsilon^2} + \hat{z}^{(m)} z ( m ) = z ∗ 1 − ϵ 2 + z ^ ( m )
Set x t ( m ) = h ( z ( m ) , x t − τ ) x_t^{(m)} = h(z^{(m)},x_{t-\tau}) x t ( m ) = h ( z ( m ) , x t − τ )
end for
Therefore, in Algorithm 5, we can prove that, if Z ( 0 ) ∼ N ( 0 , I d ) Z^{(0)} \sim N(0,I_d) Z ( 0 ) ∼ N ( 0 , I d ) Z = ( Z ( 0 ) , … , Z ( M ) ) ⊤ ∼ N ( μ M + 1 , Σ M + 1 ) \mathbf{Z} = (Z^{(0)},\ldots,Z^{(M)})^\top\sim N(\boldsymbol{\mu}_{M+1},\boldsymbol{\Sigma}_{M+1}) Z = ( Z ( 0 ) , … , Z ( M ) ) ⊤ ∼ N ( μ M + 1 , Σ M + 1 ) ρ = 1 − ϵ 2 \rho= 1-\epsilon^2 ρ = 1 − ϵ 2 z ( 0 ) , … , z ( M ) z^{(0)},\ldots,z^{(M)} z ( 0 ) , … , z ( M ) ( z ( 0 ) , … , z ( M ) ) (z^{(0)},\ldots,z^{(M)}) ( z ( 0 ) , … , z ( M ) ) x 0 ( m ) = h ( z 0 ( m ) ) x_0^{(m)} = h(z_0^{(m)}) x 0 ( m ) = h ( z 0 ( m ) ) m = 0 , … , M m = 0,\ldots,M m = 0 , … , M x 0 ( 0 ) , … , x 0 ( M ) x_0^{(0)},\ldots,x_0^{(M)} x 0 ( 0 ) , … , x 0 ( M )
The conditional particle filter within xPG for continuous state spaces follows the standard conditional particle filter, with modifications in steps 1, 2 and 8. In step 1, rather than sampling x 0 ( 1 : M ) x_0^{(1:M)} x 0 ( 1 : M ) (4) using Algorithm 5. Likewise, in step 2 and 8, the states x k τ ( 1 : M ) x_{k\tau}^{(1:M)} x k τ ( 1 : M )
Building on particle Gibbs, the exchangeable particle Gibbs method replaces the standard conditional particle filter with the exchangeable conditional particle filter . The exchangeability ensures that the detailed balance still holds for xPG, meaning this new method remains valid; see Malory [13 ] . The proof of detailed balance for the full algorithm is similar to that presented in Section F.1, so it will not be repeated here.
Exchangeable particle Gibbs
Algorithm 6 Update proposals in each time step
Require: N k ( 0 ) , μ k ( 0 ) , … , μ k ( M ) N_k^{(0)},\mu_k^{(0)},\ldots,\mu_k^{(M)} N k ( 0 ) , μ k ( 0 ) , … , μ k ( M )
N ~ k ( 0 ) ∼ Bin ( N k ( 0 ) , 1 − δ ) \widetilde{N}_k^{(0)} \sim\operatorname{Bin}(N_k^{(0)},1-\delta) N k ( 0 ) ∼ Bin ( N k ( 0 ) , 1 − δ )
Sort 0 : M 0:M 0 : M μ k ( i ) \mu_k^{(i)} μ k ( i ) i = 0 , … , M i=0,\ldots,M i = 0 , … , M S = { s i } i = − l o h i \mathcal{S}=\{s_i\}_{i=-\mathrm{lo}}^{\mathrm{hi}} S = { s i } i = − lo hi s 0 = 0 s_0=0 s 0 = 0
for i = 0 , 1 , … , h i − 1 i=0,1,\ldots,\mathrm{hi}-1 i = 0 , 1 , … , hi − 1 do
generate N ~ k ( s i + 1 ) \widetilde{N}_k^{(s_{i+1})} N k ( s i + 1 ) N ~ k ( s i ) \widetilde{N}_k^{(s_i)} N k ( s i ) Eq 8
end for
for i = 0 , − 1 , … , − l o + 1 i=0,-1,\ldots,-\mathrm{lo}+1 i = 0 , − 1 , … , − lo + 1 do
generate N ~ k ( s i − 1 ) \widetilde{N}_k^{(s_{i-1})} N k ( s i − 1 ) N ~ k ( s i ) \widetilde{N}_k^{(s_i)} N k ( s i ) Eq 9
end for
for j = 1 , … , M j=1,\ldots,M j = 1 , … , M do
Generate N k ( j ) N_k^{(j)} N k ( j ) Eq 10 .
end for
The correlation between N k ( s i ) N_k^{(s_i)} N k ( s i ) N k ( s j ) N_k^{(s_j)} N k ( s j ) i < j i<j i < j
Corr ( N k ( s i ) , N k ( s j ) ) = ( 1 − δ ) μ k ( s i ) μ k ( s i ) μ k ( s j ) . \operatorname{Corr}\left(N_k^{(s_i)},N_k^{(s_j)}\right)=\frac{(1-\delta)\mu_k^{(s_i)}}{\sqrt{\mu_k^{(s_i)}\mu_k^{(s_j)}}}. Corr ( N k ( s i ) , N k ( s j ) ) = μ k ( s i ) μ k ( s j ) ( 1 − δ ) μ k ( s i ) . where i , j = − l o , … , h i − 1 i,j=-\mathrm{lo},\ldots,\mathrm{hi}-1 i , j = − lo , … , hi − 1
Justification for xPG
Proof of Theorem 1
In this section, we prove the exchangeable particle Gibbs for MJPs satisfies detailed balance with respect to the target posterior distribution of the latent states, p ( x 0 : n ( τ ) ∣ y 1 : L , θ ) p(x_{0:n}^{(\tau)}\mid y_{1:L},\theta) p ( x 0 : n ( τ ) ∣ y 1 : L , θ ) t ℓ = ℓ t_\ell=\ell t ℓ = ℓ ℓ = 1 , … , L \ell=1,\ldots,L ℓ = 1 , … , L
We begin by introducing Lemma 4, which plays a central role in proving Theorem 1. In a conditional particle filter, at each observation time t ℓ t_\ell t ℓ M M M x ~ t ℓ ( 1 : M ) \widetilde{x}_{t_\ell}^{(1:M)} x t ℓ ( 1 : M ) t ℓ + 1 t_{\ell+1} t ℓ + 1 segment .
The lemma states that, within the propagation of a single time step of length τ \tau τ n ℓ K + 1 : ℓ K + k ( 0 ) , R 1 : v n_{\ell K+1:\ell K+k}^{(0),R_{1:v}} n ℓ K + 1 : ℓ K + k ( 0 ) , R 1 : v M M M j j j n ℓ K + 1 : ℓ K + k ( j ) , R 1 : v n_{\ell K+1:\ell K+k}^{(j),R_{1:v}} n ℓ K + 1 : ℓ K + k ( j ) , R 1 : v j = 1 , … , M j=1,\ldots,M j = 1 , … , M
Recall that we defined that x t ℓ ( a ℓ j ) = x ~ t ℓ ( j ) x_{t_\ell}^{(a_\ell^j)}=\widetilde{x}_{t_\ell}^{(j)} x t ℓ ( a ℓ j ) = x t ℓ ( j ) x t ℓ ( a ℓ j ) x_{t_\ell}^{(a_\ell^j)} x t ℓ ( a ℓ j ) t ℓ t_\ell t ℓ j j j [ ℓ K τ , ( ℓ + 1 ) K τ ] [\ell K\tau,(\ell+1)K\tau] [ ℓ K τ , ( ℓ + 1 ) K τ ] μ ℓ K + k \mu_{\ell K+k} μ ℓ K + k k k k ℓ \ell ℓ n ℓ K + k n_{\ell K+k} n ℓ K + k
Lemma 4. For any ℓ ∈ { 0 , 1 , … , L − 1 } \ell\in\{0,1,\ldots,L-1\} ℓ ∈ { 0 , 1 , … , L − 1 } k ∈ { 1 , … , K } k\in\{1,\ldots,K\} k ∈ { 1 , … , K } j ∈ { 1 , … , M } j\in\{1,\ldots,M\} j ∈ { 1 , … , M }
p ( n ℓ K + k ( 0 ) , R 1 : v ∣ x t ℓ ( 0 ) , n ℓ K + 1 : ℓ K + k − 1 ( 0 ) , R 1 : v ) × q ~ ( n ℓ K + k ( 1 : M ) , R 1 : v ∣ x t ℓ ( a ℓ 0 : M ) , n ℓ K + 1 : ℓ K + k − 1 ( 0 : M ) , R 1 : v , n ℓ K + k ( 0 ) , R 1 : v ) = p ( n ℓ K + k ( j ) , R 1 : v ∣ x t ℓ ( a ℓ j ) , n ℓ K + 1 : ℓ K + k − 1 ( j ) , R 1 : v ) × q ~ ( n ℓ K + k ( − j ) , R 1 : v ∣ x t ℓ ( a ℓ 0 : M ) , n ℓ K + 1 : ℓ K + k − 1 ( 0 : M ) , R 1 : v , n ℓ K + k ( j ) , R 1 : v ) . (19) \begin{aligned}
&p\left(n_{\ell K+k}^{(0),R_{1:v}}\mid x_{t_\ell}^{(0)},n_{\ell K+1:\ell K+k-1}^{(0),R_{1:v}}\right)
\times\widetilde{q}\left(n_{\ell K+k}^{(1:M),R_{1:v}}\mid x_{t_\ell}^{(a_\ell^{0:M})},n_{\ell K+1:\ell K+k-1}^{(0:M),R_{1:v}},n_{\ell K+k}^{(0),R_{1:v}}\right)\\
&=p\left(n_{\ell K+k}^{(j),R_{1:v}}\mid x_{t_\ell}^{(a_\ell^j)},n_{\ell K+1:\ell K+k-1}^{(j),R_{1:v}}\right)
\times\widetilde{q}\left(n_{\ell K+k}^{(-j),R_{1:v}}\mid x_{t_\ell}^{(a_\ell^{0:M})},n_{\ell K+1:\ell K+k-1}^{(0:M),R_{1:v}},n_{\ell K+k}^{(j),R_{1:v}}\right). \tag*{(19)}
\end{aligned} p ( n ℓ K + k ( 0 ) , R 1 : v ∣ x t ℓ ( 0 ) , n ℓ K + 1 : ℓ K + k − 1 ( 0 ) , R 1 : v ) × q ( n ℓ K + k ( 1 : M ) , R 1 : v ∣ x t ℓ ( a ℓ 0 : M ) , n ℓ K + 1 : ℓ K + k − 1 ( 0 : M ) , R 1 : v , n ℓ K + k ( 0 ) , R 1 : v ) = p ( n ℓ K + k ( j ) , R 1 : v ∣ x t ℓ ( a ℓ j ) , n ℓ K + 1 : ℓ K + k − 1 ( j ) , R 1 : v ) × q ( n ℓ K + k ( − j ) , R 1 : v ∣ x t ℓ ( a ℓ 0 : M ) , n ℓ K + 1 : ℓ K + k − 1 ( 0 : M ) , R 1 : v , n ℓ K + k ( j ) , R 1 : v ) . (19) Appendix F.2 provides the proof of Lemma 4.
The remainder of the proof of Theorem 1 follows similarly to the proof of the validity of the standard particle Gibbs, as presented in Fearnhead and Sherlock [6 ] .
Proof of Theorem 1. The target posterior probability mass function of the full path X ( 0 : n ) τ X_{(0:n)\tau} X ( 0 : n ) τ
p ( x ( 0 : n ) τ ∣ y 1 : L , θ ) ∝ p 0 ( x 0 ∣ θ ) ∏ k = 1 n p ( n k R 1 : ν ∣ x 0 , n 1 : ( k − 1 ) R 1 : ν ) ∏ ℓ = 1 L g ( y ℓ ∣ x t ℓ ) = p 0 ( x 0 ∣ θ ) ∏ k = 1 n p ( n k R 1 : ν ∣ μ k R 1 : ν ) ∏ ℓ = 1 L g ( y ℓ ∣ x t ℓ ) . (20) \begin{aligned}
p(x_{(0:n)\tau}\mid y_{1:L},\theta)
&\propto p_0(x_0\mid\theta)\prod_{k=1}^{n}p\left(n_k^{R_{1:\nu}}\mid x_0,n_{1:(k-1)}^{R_{1:\nu}}\right)\prod_{\ell=1}^{L}g(y_\ell\mid x_{t_\ell})\\
&=p_0(x_0\mid\theta)\prod_{k=1}^{n}p\left(n_k^{R_{1:\nu}}\mid\mu_k^{R_{1:\nu}}\right)\prod_{\ell=1}^{L}g(y_\ell\mid x_{t_\ell}). \tag*{(20)}
\end{aligned} p ( x ( 0 : n ) τ ∣ y 1 : L , θ ) ∝ p 0 ( x 0 ∣ θ ) k = 1 ∏ n p ( n k R 1 : ν ∣ x 0 , n 1 : ( k − 1 ) R 1 : ν ) ℓ = 1 ∏ L g ( y ℓ ∣ x t ℓ ) = p 0 ( x 0 ∣ θ ) k = 1 ∏ n p ( n k R 1 : ν ∣ μ k R 1 : ν ) ℓ = 1 ∏ L g ( y ℓ ∣ x t ℓ ) . (20) where n 1 : 0 = ∅ n_{1:0}=\emptyset n 1 : 0 = ∅ x T ( i ) x_T^{(i)} x T ( i )
q ~ 0 ( x 0 ( 1 : M ) ∣ x 0 ( 0 ) ) × ∏ ℓ = 0 L − 1 ∏ k = 1 K q ~ ( n ℓ K + k ( 1 : M ) , R 1 : ν ∣ μ ℓ K + k ( 0 : M ) , R 1 : ν , n ℓ K + k ( 0 ) , R 1 : ν ) { ∏ ℓ = 1 L − 1 ∏ m = 1 M g ( y ℓ ∣ x t ℓ ( a ℓ m ) ) ∑ j = 0 M g ( y ℓ ∣ x t ℓ ( j ) ) } g ( y L ∣ x t L ( i ) ) ∑ j = 0 M g ( y L ∣ x t L ( j ) ) . (21) \widetilde{q}_0\left(x_0^{(1:M)}\mid x_0^{(0)}\right)\times\prod_{\ell=0}^{L-1}\prod_{k=1}^{K}\widetilde{q}\left(n_{\ell K+k}^{(1:M),R_{1:\nu}}\mid\mu_{\ell K+k}^{(0:M),R_{1:\nu}},n_{\ell K+k}^{(0),R_{1:\nu}}\right)
\left\{\prod_{\ell=1}^{L-1}\prod_{m=1}^{M}\frac{g\left(y_\ell\mid x_{t_\ell}^{(a_\ell^m)}\right)}{\sum_{j=0}^{M}g\left(y_\ell\mid x_{t_\ell}^{(j)}\right)}\right\}\frac{g\left(y_L\mid x_{t_L}^{(i)}\right)}{\sum_{j=0}^{M}g\left(y_L\mid x_{t_L}^{(j)}\right)}. \tag*{(21)} q 0 ( x 0 ( 1 : M ) ∣ x 0 ( 0 ) ) × ℓ = 0 ∏ L − 1 k = 1 ∏ K q ( n ℓ K + k ( 1 : M ) , R 1 : ν ∣ μ ℓ K + k ( 0 : M ) , R 1 : ν , n ℓ K + k ( 0 ) , R 1 : ν ) ⎩ ⎨ ⎧ ℓ = 1 ∏ L − 1 m = 1 ∏ M ∑ j = 0 M g ( y ℓ ∣ x t ℓ ( j ) ) g ( y ℓ ∣ x t ℓ ( a ℓ m ) ) ⎭ ⎬ ⎫ ∑ j = 0 M g ( y L ∣ x t L ( j ) ) g ( y L ∣ x t L ( i ) ) . (21) where a ℓ 0 = 0 a_\ell^0=0 a ℓ 0 = 0 ℓ = 0 , … , L \ell=0,\ldots,L ℓ = 0 , … , L
Based on (20) and (21) , the product of the posterior density of the current reference path and the probability of generating all the random variables within a single run of the conditional particle filter, followed by selecting the new path back from x T ( i ) x_T^{(i)} x T ( i )
The first term corresponds to the prior density of the initial state and the proposal used to generate the remaining initial particles,
p 0 ( x 0 ( 0 ) ∣ θ ) q ~ 0 ( x 0 ( 1 : M ) ∣ x 0 ( 0 ) ) . (22) p_0\left(x_0^{(0)}\mid\theta\right)\widetilde{q}_0\left(x_0^{(1:M)}\mid x_0^{(0)}\right). \tag*{(22)} p 0 ( x 0 ( 0 ) ∣ θ ) q 0 ( x 0 ( 1 : M ) ∣ x 0 ( 0 ) ) . (22) The second term accounts for the generation of reaction counts within each interval,
∏ ℓ = 0 L − 1 ∏ k = 1 K p ( n ℓ K + k ( 0 ) , R 1 : ν ∣ μ ℓ K + k ( 0 ) , R 1 : ν ) × q ~ ( n ℓ K + k ( 1 : M ) , R 1 : ν ∣ μ ℓ K + k ( 0 : M ) , R 1 : ν , n ℓ K + k ( 0 ) , R 1 : ν ) . (23) \prod_{\ell=0}^{L-1}\prod_{k=1}^{K}p\left(n_{\ell K+k}^{(0),R_{1:\nu}}\mid\mu_{\ell K+k}^{(0),R_{1:\nu}}\right)\times\widetilde{q}\left(n_{\ell K+k}^{(1:M),R_{1:\nu}}\mid\mu_{\ell K+k}^{(0:M),R_{1:\nu}},n_{\ell K+k}^{(0),R_{1:\nu}}\right). \tag*{(23)} ℓ = 0 ∏ L − 1 k = 1 ∏ K p ( n ℓ K + k ( 0 ) , R 1 : ν ∣ μ ℓ K + k ( 0 ) , R 1 : ν ) × q ( n ℓ K + k ( 1 : M ) , R 1 : ν ∣ μ ℓ K + k ( 0 : M ) , R 1 : ν , n ℓ K + k ( 0 ) , R 1 : ν ) . (23) Finally, the third term corresponds to the contribution of the observation likelihoods and the path resampling probabilities,
∏ ℓ = 1 L g ( y ℓ ∣ x t ℓ ( 0 ) ) [ ∏ ℓ = 1 L − 1 ∏ m = 1 M g ( y ℓ ∣ x t ℓ ( a ℓ m ) ) ∑ j = 0 M g ( y ℓ ∣ x t ℓ ( j ) ) ] × g ( y L ∣ x t L ( i ) ) ∑ j = 0 M g ( y L ∣ x t L ( j ) ) . (24) \prod_{\ell=1}^{L}g(y_\ell\mid x_{t_\ell}^{(0)})\left[\prod_{\ell=1}^{L-1}\prod_{m=1}^{M}\frac{g\left(y_\ell\mid x_{t_\ell}^{(a_\ell^m)}\right)}{\sum_{j=0}^{M}g\left(y_\ell\mid x_{t_\ell}^{(j)}\right)}\right]\times\frac{g\left(y_L\mid x_{t_L}^{(i)}\right)}{\sum_{j=0}^{M}g\left(y_L\mid x_{t_L}^{(j)}\right)}. \tag*{(24)} ℓ = 1 ∏ L g ( y ℓ ∣ x t ℓ ( 0 ) ) ℓ = 1 ∏ L − 1 m = 1 ∏ M ∑ j = 0 M g ( y ℓ ∣ x t ℓ ( j ) ) g ( y ℓ ∣ x t ℓ ( a ℓ m ) ) × ∑ j = 0 M g ( y L ∣ x t L ( j ) ) g ( y L ∣ x t L ( i ) ) . (24) From (5) we know that (22) is invariant to x 0 ( 0 ) ↔ x 0 ( b 0 i ) x_0^{(0)}\leftrightarrow x_0^{(b_0^i)} x 0 ( 0 ) ↔ x 0 ( b 0 i ) (23) is equal to
∏ ℓ = 0 L − 1 ∏ k = 1 K p ( n ℓ K + k ( b ℓ + 1 i ) , R 1 : ν ∣ μ ℓ K + k ( b ℓ + 1 i ) , R 1 : ν ) × q ~ ( n ℓ K + k ( − b ℓ + 1 i ) , R 1 : ν ∣ μ ℓ K + k ( 0 : M ) , R 1 : ν , n ℓ K + k ( b ℓ + 1 i ) , R 1 : ν ) . (25) \prod_{\ell=0}^{L-1}\prod_{k=1}^{K}p\left(n_{\ell K+k}^{(b_{\ell+1}^i),R_{1:\nu}}\mid\mu_{\ell K+k}^{(b_{\ell+1}^i),R_{1:\nu}}\right)\times\widetilde{q}\left(n_{\ell K+k}^{(-b_{\ell+1}^i),R_{1:\nu}}\mid\mu_{\ell K+k}^{(0:M),R_{1:\nu}},n_{\ell K+k}^{(b_{\ell+1}^i),R_{1:\nu}}\right). \tag*{(25)} ℓ = 0 ∏ L − 1 k = 1 ∏ K p ( n ℓ K + k ( b ℓ + 1 i ) , R 1 : ν ∣ μ ℓ K + k ( b ℓ + 1 i ) , R 1 : ν ) × q ( n ℓ K + k ( − b ℓ + 1 i ) , R 1 : ν ∣ μ ℓ K + k ( 0 : M ) , R 1 : ν , n ℓ K + k ( b ℓ + 1 i ) , R 1 : ν ) . (25) In applying Lemma 4, note that a ℓ + 1 b ℓ + 1 i = b ℓ i a_{\ell+1}^{b_{\ell+1}^i}=b_\ell^i a ℓ + 1 b ℓ + 1 i = b ℓ i a ℓ j = b ℓ i a_\ell^j=b_\ell^i a ℓ j = b ℓ i j = b ℓ + 1 i j=b_{\ell+1}^i j = b ℓ + 1 i ℓ = 1 , … , L − 1 \ell=1,\ldots,L-1 ℓ = 1 , … , L − 1 a ℓ 0 = 0 a_\ell^0=0 a ℓ 0 = 0 ℓ = 1 , … , L − 1 \ell=1,\ldots,L-1 ℓ = 1 , … , L − 1 (24) can be rewritten as
g ( y L ∣ x t L ( 0 ) ) g ( y L ∣ x t L ( i ) ) ∏ ℓ = 1 L − 1 ∏ m = 0 M g ( y ℓ ∣ x t ℓ ( a ℓ m ) ) ( ∏ ℓ = 1 L − 1 ∑ j = 0 M g ( y ℓ ∣ x t ℓ ( j ) ) ) M ∑ j = 0 M g ( y L ∣ x t L ( j ) ) \frac{g\left(y_L\mid x_{t_L}^{(0)}\right)g\left(y_L\mid x_{t_L}^{(i)}\right)\prod_{\ell=1}^{L-1}\prod_{m=0}^{M}g\left(y_\ell\mid x_{t_\ell}^{(a_\ell^m)}\right)}
{\left(\prod_{\ell=1}^{L-1}\sum_{j=0}^{M}g\left(y_\ell\mid x_{t_\ell}^{(j)}\right)\right)^{M}\sum_{j=0}^{M}g\left(y_L\mid x_{t_L}^{(j)}\right)} ( ∏ ℓ = 1 L − 1 ∑ j = 0 M g ( y ℓ ∣ x t ℓ ( j ) ) ) M ∑ j = 0 M g ( y L ∣ x t L ( j ) ) g ( y L ∣ x t L ( 0 ) ) g ( y L ∣ x t L ( i ) ) ∏ ℓ = 1 L − 1 ∏ m = 0 M g ( y ℓ ∣ x t ℓ ( a ℓ m ) ) which is invariant to x t L ( 0 ) ↔ x t L ( i ) x_{t_L}^{(0)}\leftrightarrow x_{t_L}^{(i)} x t L ( 0 ) ↔ x t L ( i ) i ∈ { 1 , … , M } i\in\{1,\ldots,M\} i ∈ { 1 , … , M } { x k τ ( 0 ) } k = 0 n \{x_{k\tau}^{(0)}\}_{k=0}^{n} { x k τ ( 0 ) } k = 0 n x T ( i ) x_T^{(i)} x T ( i ) x ( 0 : n ) τ ( 0 ) ↔ x ( 0 : n ) τ ( b i ) x_{(0:n)\tau}^{(0)}\leftrightarrow x_{(0:n)\tau}^{(b^i)} x ( 0 : n ) τ ( 0 ) ↔ x ( 0 : n ) τ ( b i ) b L i = i b_L^i=i b L i = i p ( x ( 0 : n ) τ ∣ y 0 : L , θ ) p(x_{(0:n)\tau}\mid y_{0:L},\theta) p ( x ( 0 : n ) τ ∣ y 0 : L , θ )
Proof of Lemma 4
Firstly, for each k ∈ { 1 , … , K } k \in\{1,\ldots,K\} k ∈ { 1 , … , K } x t ( 0 : M ) ℓ ( a ℓ 0 : M ) \boldsymbol{x}_{t}^{(0:M)}{}_{\ell}^{(a_{\ell}^{0:M})} x t ( 0 : M ) ℓ ( a ℓ 0 : M ) n k ( 0 : M ) , R 1 : v \boldsymbol{n}_{k}^{(0:M),R_{1:v}} n k ( 0 : M ) , R 1 : v μ ℓ K + k ( 0 : M ) , R 1 : v \boldsymbol{\mu}_{\ell K+k}^{(0:M),R_{1:v}} μ ℓ K + k ( 0 : M ) , R 1 : v j ∈ { 0 , … , M } j \in\{0,\ldots,M\} j ∈ { 0 , … , M } ( μ ℓ K + k ( 0 : M ) , R 1 : v , n ℓ K + k ( j ) , R 1 : v ) \left(\boldsymbol{\mu}_{\ell K+k}^{(0:M),R_{1:v}},\boldsymbol{n}_{\ell K+k}^{(j),R_{1:v}}\right) ( μ ℓ K + k ( 0 : M ) , R 1 : v , n ℓ K + k ( j ) , R 1 : v ) n ℓ K + k ( − j ) , R 1 : v \boldsymbol{n}_{\ell K+k}^{(-j),R_{1:v}} n ℓ K + k ( − j ) , R 1 : v ( x t ( a ℓ 0 : M ) , n ℓ K + 1 : ℓ K + k − 1 ( 0 : M ) , R 1 : v , n ℓ K + k ( j ) , R 1 : v ) \left(\boldsymbol{x}_{t}^{(a_{\ell}^{0:M})},\boldsymbol{n}_{\ell K+1:\ell K+k-1}^{(0:M),R_{1:v}},\boldsymbol{n}_{\ell K+k}^{(j),R_{1:v}}\right) ( x t ( a ℓ 0 : M ) , n ℓ K + 1 : ℓ K + k − 1 ( 0 : M ) , R 1 : v , n ℓ K + k ( j ) , R 1 : v ) ℓ = 0 \ell=0 ℓ = 0
p ( n k ( 0 ) , R 1 : v ∣ μ k ( 0 ) , R 1 : v ) q ( n k ( 1 : M ) , R 1 : v ∣ μ k ( 0 : M ) , R 1 : v , n k ( 0 ) , R 1 : v ) = p ( n k ( j ) , R 1 : v ∣ μ k ( j ) , R 1 : v ) q ( n k ( − j ) , R 1 : v ∣ μ k ( 0 : M ) , R 1 : v , n k ( j ) , R 1 : v ) . p\left(\boldsymbol{n}_{k}^{(0),R_{1:v}}\mid\boldsymbol{\mu}_{k}^{(0),R_{1:v}}\right)q\left(\boldsymbol{n}_{k}^{(1:M),R_{1:v}}\mid\boldsymbol{\mu}_{k}^{(0:M),R_{1:v}},\boldsymbol{n}_{k}^{(0),R_{1:v}}\right)=p\left(\boldsymbol{n}_{k}^{(j),R_{1:v}}\mid\boldsymbol{\mu}_{k}^{(j),R_{1:v}}\right)q\left(\boldsymbol{n}_{k}^{(-j),R_{1:v}}\mid\boldsymbol{\mu}_{k}^{(0:M),R_{1:v}},\boldsymbol{n}_{k}^{(j),R_{1:v}}\right). p ( n k ( 0 ) , R 1 : v ∣ μ k ( 0 ) , R 1 : v ) q ( n k ( 1 : M ) , R 1 : v ∣ μ k ( 0 : M ) , R 1 : v , n k ( 0 ) , R 1 : v ) = p ( n k ( j ) , R 1 : v ∣ μ k ( j ) , R 1 : v ) q ( n k ( − j ) , R 1 : v ∣ μ k ( 0 : M ) , R 1 : v , n k ( j ) , R 1 : v ) . Conditional on μ k ( 0 : M ) \boldsymbol{\mu}_{k}^{(0:M)} μ k ( 0 : M ) k k k
For a given reaction type, let q t h i n ( n ~ ( 0 ) ∣ n ( 0 ) ) , q R ( n ~ ( s i + 1 ) ∣ n ~ ( s i ) ) , q L ( n ~ ( s i − 1 ) ∣ n ~ ( s i ) ) q_{\mathrm{thin}}(\widetilde{n}^{(0)}\mid n^{(0)}), q_{\mathrm{R}}(\widetilde{n}^{(s_{i+1})}\mid\widetilde{n}^{(s_i)}), q_{\mathrm{L}}(\widetilde{n}^{(s_{i-1})}\mid\widetilde{n}^{(s_i)}) q thin ( n ( 0 ) ∣ n ( 0 ) ) , q R ( n ( s i + 1 ) ∣ n ( s i ) ) , q L ( n ( s i − 1 ) ∣ n ( s i ) ) q f i l l ( n ( j ) ∣ n ~ ( j ) ) q_{\mathrm{fill}}(n^{(j)}\mid\widetilde{n}^{(j)}) q fill ( n ( j ) ∣ n ( j ) ) μ ( j ) ≥ μ ( 0 ) \mu^{(j)}\geq\mu^{(0)} μ ( j ) ≥ μ ( 0 ) j = s r j=s_r j = s r r ≥ 0 r\geq0 r ≥ 0 0 ↔ j 0\leftrightarrow j 0 ↔ j
p ( n ( 0 ) ∣ μ ( 0 ) ) q t h i n ( n ~ ( 0 ) ∣ n ( 0 ) ) ∏ i = 0 h i − 1 q R ( n ~ ( s i + 1 ) ∣ n ~ ( s i ) ) ∏ i = − l o + 1 0 q L ( n ~ ( s i − 1 ) ∣ n ~ ( s i ) ) ∏ i = 1 M q f i l l ( n ( i ) ∣ n ~ ( i ) ) = p ( n ( j ) ∣ μ ( j ) ) q t h i n ( n ~ ( j ) ∣ n ( j ) ) ∏ i = r h i − 1 q R ( n ~ ( s i + 1 ) ∣ n ~ ( s i ) ) ∏ i = − l o + 1 r q L ( n ~ ( s i − 1 ) ∣ n ~ ( s i ) ) ∏ i = 0 i ≠ j M q f i l l ( n ( i ) ∣ n ~ ( i ) ) . \begin{aligned}
&p(n^{(0)}\mid\mu^{(0)})q_{\mathrm{thin}}(\widetilde{n}^{(0)}\mid n^{(0)})\prod_{i=0}^{h_i-1}q_{\mathrm{R}}(\widetilde{n}^{(s_{i+1})}\mid\widetilde{n}^{(s_i)})\prod_{i=-l_o+1}^{0}q_{\mathrm{L}}(\widetilde{n}^{(s_{i-1})}\mid\widetilde{n}^{(s_i)})\prod_{i=1}^{M}q_{\mathrm{fill}}(n^{(i)}\mid\widetilde{n}^{(i)})\\
&=p(n^{(j)}\mid\mu^{(j)})q_{\mathrm{thin}}(\widetilde{n}^{(j)}\mid n^{(j)})\prod_{i=r}^{h_i-1}q_{\mathrm{R}}(\widetilde{n}^{(s_{i+1})}\mid\widetilde{n}^{(s_i)})\prod_{i=-l_o+1}^{r}q_{\mathrm{L}}(\widetilde{n}^{(s_{i-1})}\mid\widetilde{n}^{(s_i)})\prod_{\substack{i=0\\i\ne j}}^{M}q_{\mathrm{fill}}(n^{(i)}\mid\widetilde{n}^{(i)}).
\end{aligned} p ( n ( 0 ) ∣ μ ( 0 ) ) q thin ( n ( 0 ) ∣ n ( 0 ) ) i = 0 ∏ h i − 1 q R ( n ( s i + 1 ) ∣ n ( s i ) ) i = − l o + 1 ∏ 0 q L ( n ( s i − 1 ) ∣ n ( s i ) ) i = 1 ∏ M q fill ( n ( i ) ∣ n ( i ) ) = p ( n ( j ) ∣ μ ( j ) ) q thin ( n ( j ) ∣ n ( j ) ) i = r ∏ h i − 1 q R ( n ( s i + 1 ) ∣ n ( s i ) ) i = − l o + 1 ∏ r q L ( n ( s i − 1 ) ∣ n ( s i ) ) i = 0 i = j ∏ M q fill ( n ( i ) ∣ n ( i ) ) . Integrating out n ~ ( 0 : M ) \widetilde{n}^{(0:M)} n ( 0 : M )
Both sides of the equation have the common terms ∏ i = r h i − 1 q R ( n ~ ( s i + 1 ) ∣ n ~ ( s i ) ) \prod_{i=r}^{h_i-1}q_{\mathrm{R}}(\widetilde{n}^{(s_{i+1})}\mid\widetilde{n}^{(s_i)}) ∏ i = r h i − 1 q R ( n ( s i + 1 ) ∣ n ( s i ) ) ∏ i = − l o + 1 0 q L ( n ~ ( s i − 1 ) ∣ n ~ ( s i ) ) \prod_{i=-l_o+1}^{0}q_{\mathrm{L}}(\widetilde{n}^{(s_{i-1})}\mid\widetilde{n}^{(s_i)}) ∏ i = − l o + 1 0 q L ( n ( s i − 1 ) ∣ n ( s i ) ) ∏ i = 1 , i ≠ j M q f i l l ( n ( i ) ∣ n ~ ( i ) ) \prod_{i=1,i\ne j}^{M}q_{\mathrm{fill}}(n^{(i)}\mid\widetilde{n}^{(i)}) ∏ i = 1 , i = j M q fill ( n ( i ) ∣ n ( i ) ) p ( n ∣ μ ) ≡ Pois ( n ; μ ) p(n\mid\mu)\equiv\operatorname{Pois}(n;\mu) p ( n ∣ μ ) ≡ Pois ( n ; μ )
Pois ( n ( 0 ) ; μ ( 0 ) ) q t h i n ( n ~ ( 0 ) ∣ n ( 0 ) ) { ∏ i = 0 r − 1 q R ( n ~ ( s i + 1 ) ∣ n ~ ( s i ) ) } q f i l l ( n ( j ) ∣ n ~ ( j ) ) = Pois ( n ( j ) ; μ ( j ) ) q t h i n ( n ~ ( j ) ∣ n ( j ) ) { ∏ i = 1 r q L ( n ~ ( s i − 1 ) ∣ n ~ ( s i ) ) } q f i l l ( n ( 0 ) ∣ n ~ ( 0 ) ) . (26) \begin{aligned}
&\operatorname{Pois}(n^{(0)};\mu^{(0)})q_{\mathrm{thin}}(\widetilde{n}^{(0)}\mid n^{(0)})\left\{\prod_{i=0}^{r-1}q_{\mathrm{R}}(\widetilde{n}^{(s_{i+1})}\mid\widetilde{n}^{(s_i)})\right\}q_{\mathrm{fill}}(n^{(j)}\mid\widetilde{n}^{(j)})\\
&=\operatorname{Pois}(n^{(j)};\mu^{(j)})q_{\mathrm{thin}}(\widetilde{n}^{(j)}\mid n^{(j)})\left\{\prod_{i=1}^{r}q_{\mathrm{L}}(\widetilde{n}^{(s_{i-1})}\mid\widetilde{n}^{(s_i)})\right\}q_{\mathrm{fill}}(n^{(0)}\mid\widetilde{n}^{(0)}). \tag*{(26)}
\end{aligned} Pois ( n ( 0 ) ; μ ( 0 ) ) q thin ( n ( 0 ) ∣ n ( 0 ) ) { i = 0 ∏ r − 1 q R ( n ( s i + 1 ) ∣ n ( s i ) ) } q fill ( n ( j ) ∣ n ( j ) ) = Pois ( n ( j ) ; μ ( j ) ) q thin ( n ( j ) ∣ n ( j ) ) { i = 1 ∏ r q L ( n ( s i − 1 ) ∣ n ( s i ) ) } q fill ( n ( 0 ) ∣ n ( 0 ) ) . (26) Now, defining μ ~ ( j ) : = ( 1 − δ ) μ ( j ) \widetilde{\mu}^{(j)}:=(1-\delta)\mu^{(j)} μ ( j ) := ( 1 − δ ) μ ( j )
q R ( n ~ ( s i ) ∣ n ~ ( s i − 1 ) ) = ( μ ~ ( s i ) − μ ~ ( s i − 1 ) ) n ~ ( s i ) − n ~ ( s i − 1 ) ( n ~ ( s i ) − n ~ ( s i − 1 ) ) ! exp { − ( μ ~ ( s i ) − μ ~ ( s i − 1 ) ) } . q_{\mathrm{R}}(\widetilde{n}^{(s_i)}\mid\widetilde{n}^{(s_{i-1})})=
\frac{\left(\widetilde{\mu}^{(s_i)}-\widetilde{\mu}^{(s_{i-1})}\right)^{\widetilde{n}^{(s_i)}-\widetilde{n}^{(s_{i-1})}}}{\left(\widetilde{n}^{(s_i)}-\widetilde{n}^{(s_{i-1})}\right)!}
\exp\left\{-\left(\widetilde{\mu}^{(s_i)}-\widetilde{\mu}^{(s_{i-1})}\right)\right\}. q R ( n ( s i ) ∣ n ( s i − 1 ) ) = ( n ( s i ) − n ( s i − 1 ) ) ! ( μ ( s i ) − μ ( s i − 1 ) ) n ( s i ) − n ( s i − 1 ) exp { − ( μ ( s i ) − μ ( s i − 1 ) ) } . and
q L ( n ~ ( s i − 1 ) ∣ n ~ ( s i ) ) = n ~ ( s i ) ! n ~ ( s i − 1 ) ! ( n ~ ( s i ) − n ~ ( s i − 1 ) ) ! ( μ ~ ( s i − 1 ) μ ~ ( s i ) ) n ~ ( s i − 1 ) ( μ ~ ( s i ) − μ ~ ( s i − 1 ) μ ~ ( s i ) ) n ~ ( s i ) − n ~ ( s i − 1 ) = q R ( n ~ ( s i ) ∣ n ~ ( s i − 1 ) ) n ~ ( s i ) ! n ~ ( s i − 1 ) ! ( μ ~ ( s i − 1 ) ) n ~ ( s i − 1 ) ( μ ~ ( s i ) ) n ~ ( s i ) exp ( μ ~ ( s i ) − μ ~ ( s i − 1 ) ) . \begin{aligned}
q_{\mathrm{L}}(\widetilde{n}^{(s_{i-1})}\mid\widetilde{n}^{(s_i)})
&=\frac{\widetilde{n}^{(s_i)}!}{\widetilde{n}^{(s_{i-1})}!\left(\widetilde{n}^{(s_i)}-\widetilde{n}^{(s_{i-1})}\right)!}
\left(\frac{\widetilde{\mu}^{(s_{i-1})}}{\widetilde{\mu}^{(s_i)}}\right)^{\widetilde{n}^{(s_{i-1})}}
\left(\frac{\widetilde{\mu}^{(s_i)}-\widetilde{\mu}^{(s_{i-1})}}{\widetilde{\mu}^{(s_i)}}\right)^{\widetilde{n}^{(s_i)}-\widetilde{n}^{(s_{i-1})}}\\
&=q_{\mathrm{R}}(\widetilde{n}^{(s_i)}\mid\widetilde{n}^{(s_{i-1})})\frac{\widetilde{n}^{(s_i)}!}{\widetilde{n}^{(s_{i-1})}!}\frac{\left(\widetilde{\mu}^{(s_{i-1})}\right)^{\widetilde{n}^{(s_{i-1})}}}{\left(\widetilde{\mu}^{(s_i)}\right)^{\widetilde{n}^{(s_i)}}}\exp\left(\widetilde{\mu}^{(s_i)}-\widetilde{\mu}^{(s_{i-1})}\right).
\end{aligned} q L ( n ( s i − 1 ) ∣ n ( s i ) ) = n ( s i − 1 ) ! ( n ( s i ) − n ( s i − 1 ) ) ! n ( s i ) ! ( μ ( s i ) μ ( s i − 1 ) ) n ( s i − 1 ) ( μ ( s i ) μ ( s i ) − μ ( s i − 1 ) ) n ( s i ) − n ( s i − 1 ) = q R ( n ( s i ) ∣ n ( s i − 1 ) ) n ( s i − 1 ) ! n ( s i ) ! ( μ ( s i ) ) n ( s i ) ( μ ( s i − 1 ) ) n ( s i − 1 ) exp ( μ ( s i ) − μ ( s i − 1 ) ) . So,
∏ i = 1 r q L ( n ~ ( s i − 1 ) ∣ n ~ ( s i ) ) = { ∏ i = 1 r q R ( n ~ ( s i ) ∣ n ~ ( s i − 1 ) ) } × n ~ ( s r ) ! n ~ ( s 0 ) ! ( μ ~ ( s 0 ) ) n ~ ( s 0 ) ( μ ~ ( s r ) ) n ~ ( s r ) exp ( μ ~ ( s r ) − μ ~ ( s 0 ) ) = { ∏ i = 1 r q R ( n ~ ( s i ) ∣ n ~ ( s i − 1 ) ) } × n ~ ( j ) ! n ~ ( 0 ) ! ( μ ~ ( 0 ) ) n ~ ( 0 ) ( μ ~ ( j ) ) n ~ ( j ) exp ( μ ~ ( j ) − μ ~ ( 0 ) ) = { ∏ i = 1 r q R ( n ~ ( s i ) ∣ n ~ ( s i − 1 ) ) } × n ~ ( j ) ! n ~ ( 0 ) ! ( μ ( 0 ) ) n ~ ( 0 ) ( μ ( j ) ) n ~ ( j ) ( 1 − δ ) n ~ ( 0 ) − n ~ ( j ) × exp { ( 1 − δ ) μ ( j ) − ( 1 − δ ) μ ( 0 ) } . \begin{aligned}
\prod_{i=1}^{r}q_{\mathrm{L}}(\widetilde{n}^{(s_{i-1})}\mid\widetilde{n}^{(s_i)})
&=\left\{\prod_{i=1}^{r}q_{\mathrm{R}}(\widetilde{n}^{(s_i)}\mid\widetilde{n}^{(s_{i-1})})\right\}
\times\frac{\widetilde{n}^{(s_r)}!}{\widetilde{n}^{(s_0)}!}
\frac{\left(\widetilde{\mu}^{(s_0)}\right)^{\widetilde{n}^{(s_0)}}}{\left(\widetilde{\mu}^{(s_r)}\right)^{\widetilde{n}^{(s_r)}}}
\exp\left(\widetilde{\mu}^{(s_r)}-\widetilde{\mu}^{(s_0)}\right)\\
&=\left\{\prod_{i=1}^{r}q_{\mathrm{R}}(\widetilde{n}^{(s_i)}\mid\widetilde{n}^{(s_{i-1})})\right\}
\times\frac{\widetilde{n}^{(j)}!}{\widetilde{n}^{(0)}!}
\frac{\left(\widetilde{\mu}^{(0)}\right)^{\widetilde{n}^{(0)}}}{\left(\widetilde{\mu}^{(j)}\right)^{\widetilde{n}^{(j)}}}
\exp\left(\widetilde{\mu}^{(j)}-\widetilde{\mu}^{(0)}\right)\\
&=\left\{\prod_{i=1}^{r}q_{\mathrm{R}}(\widetilde{n}^{(s_i)}\mid\widetilde{n}^{(s_{i-1})})\right\}
\times\frac{\widetilde{n}^{(j)}!}{\widetilde{n}^{(0)}!}
\frac{\left(\mu^{(0)}\right)^{\widetilde{n}^{(0)}}}{\left(\mu^{(j)}\right)^{\widetilde{n}^{(j)}}}
(1-\delta)^{\widetilde{n}^{(0)}-\widetilde{n}^{(j)}}\\
&\qquad\times\exp\left\{(1-\delta)\mu^{(j)}-(1-\delta)\mu^{(0)}\right\}.
\end{aligned} i = 1 ∏ r q L ( n ( s i − 1 ) ∣ n ( s i ) ) = { i = 1 ∏ r q R ( n ( s i ) ∣ n ( s i − 1 ) ) } × n ( s 0 ) ! n ( s r ) ! ( μ ( s r ) ) n ( s r ) ( μ ( s 0 ) ) n ( s 0 ) exp ( μ ( s r ) − μ ( s 0 ) ) = { i = 1 ∏ r q R ( n ( s i ) ∣ n ( s i − 1 ) ) } × n ( 0 ) ! n ( j ) ! ( μ ( j ) ) n ( j ) ( μ ( 0 ) ) n ( 0 ) exp ( μ ( j ) − μ ( 0 ) ) = { i = 1 ∏ r q R ( n ( s i ) ∣ n ( s i − 1 ) ) } × n ( 0 ) ! n ( j ) ! ( μ ( j ) ) n ( j ) ( μ ( 0 ) ) n ( 0 ) ( 1 − δ ) n ( 0 ) − n ( j ) × exp { ( 1 − δ ) μ ( j ) − ( 1 − δ ) μ ( 0 ) } . Let
b : = ( μ ( 0 ) ) n ( 0 ) ( μ ( j ) ) n ( j ) δ n ( 0 ) − n ~ ( 0 ) + n ( j ) − n ~ ( j ) exp { − μ ( 0 ) − μ ( j ) } / [ ( n ( 0 ) − n ~ ( 0 ) ) ! ( n ( j ) − n ~ ( j ) ) ! ] . b := \left(\mu^{(0)}\right)^{n^{(0)}}\left(\mu^{(j)}\right)^{n^{(j)}}\delta^{n^{(0)}-\widetilde{n}^{(0)}+n^{(j)}-\widetilde{n}^{(j)}}\exp\left\{-\mu^{(0)}-\mu^{(j)}\right\}/\left[\left(n^{(0)}-\widetilde{n}^{(0)}\right)!\left(n^{(j)}-\widetilde{n}^{(j)}\right)!\right]. b := ( μ ( 0 ) ) n ( 0 ) ( μ ( j ) ) n ( j ) δ n ( 0 ) − n ( 0 ) + n ( j ) − n ( j ) exp { − μ ( 0 ) − μ ( j ) } / [ ( n ( 0 ) − n ( 0 ) ) ! ( n ( j ) − n ( j ) ) ! ] . The multiplier of ∏ i = 1 r q L ( n ~ ( s i − 1 ) ∣ n ~ ( s i ) ) \prod_{i=1}^{r}q_L\left(\widetilde{n}^{(s_{i-1})}\mid\widetilde{n}^{(s_i)}\right) ∏ i = 1 r q L ( n ( s i − 1 ) ∣ n ( s i ) ) (26) is
( μ ( j ) ) n ( j ) n ( j ) ! exp { − μ ( j ) } n ( j ) ! n ~ ( j ) ! ( n ( j ) − n ~ ( j ) ) ! ( 1 − δ ) n ~ ( j ) δ n ( j ) − n ~ ( j ) ( μ ( 0 ) δ ) n ( 0 ) − n ~ ( 0 ) ( n ( 0 ) − n ~ ( 0 ) ) ! exp { − δ μ ( 0 ) } = b × 1 n ~ ( j ) ! ( 1 − δ ) n ~ ( j ) ( μ ( 0 ) ) − n ~ ( 0 ) exp { ( 1 − δ ) μ ( 0 ) } . \begin{aligned}
\frac{\left(\mu^{(j)}\right)^{n^{(j)}}}{n^{(j)}!}\exp\left\{-\mu^{(j)}\right\}\frac{n^{(j)}!}{\widetilde{n}^{(j)}!\left(n^{(j)}-\widetilde{n}^{(j)}\right)!}(1-\delta)^{\widetilde{n}^{(j)}}\delta^{n^{(j)}-\widetilde{n}^{(j)}}\frac{\left(\mu^{(0)}\delta\right)^{n^{(0)}-\widetilde{n}^{(0)}}}{\left(n^{(0)}-\widetilde{n}^{(0)}\right)!}\exp\left\{-\delta\mu^{(0)}\right\} \\
&= b\times\frac{1}{\widetilde{n}^{(j)}!}(1-\delta)^{\widetilde{n}^{(j)}}\left(\mu^{(0)}\right)^{-\widetilde{n}^{(0)}}\exp\left\{(1-\delta)\mu^{(0)}\right\}.
\end{aligned} n ( j ) ! ( μ ( j ) ) n ( j ) exp { − μ ( j ) } n ( j ) ! ( n ( j ) − n ( j ) ) ! n ( j ) ! ( 1 − δ ) n ( j ) δ n ( j ) − n ( j ) ( n ( 0 ) − n ( 0 ) ) ! ( μ ( 0 ) δ ) n ( 0 ) − n ( 0 ) exp { − δ μ ( 0 ) } = b × n ( j ) ! 1 ( 1 − δ ) n ( j ) ( μ ( 0 ) ) − n ( 0 ) exp { ( 1 − δ ) μ ( 0 ) } . That is, the right hand side of (26) is
b × ∏ i = 1 r q L ( n ~ ( s i − 1 ) ∣ n ~ ( s i ) ) × 1 n ~ ( j ) ! ( 1 − δ ) n ~ ( j ) ( μ ( 0 ) ) − n ~ ( 0 ) exp { ( 1 − δ ) μ ( 0 ) } . b\times\prod_{i=1}^{r}q_L\left(\widetilde{n}^{(s_{i-1})}\mid\widetilde{n}^{(s_i)}\right)\times\frac{1}{\widetilde{n}^{(j)}!}(1-\delta)^{\widetilde{n}^{(j)}}\left(\mu^{(0)}\right)^{-\widetilde{n}^{(0)}}\exp\left\{(1-\delta)\mu^{(0)}\right\}. b × i = 1 ∏ r q L ( n ( s i − 1 ) ∣ n ( s i ) ) × n ( j ) ! 1 ( 1 − δ ) n ( j ) ( μ ( 0 ) ) − n ( 0 ) exp { ( 1 − δ ) μ ( 0 ) } . However, by (27) , this is
b × ∏ i = 1 r q R ( n ~ ( s i ) ∣ n ~ ( s i − 1 ) ) × 1 n ~ ( 0 ) ! ( 1 − δ ) n ~ ( 0 ) ( μ ( j ) ) − n ~ ( j ) exp { ( 1 − δ ) μ ( j ) } . b\times\prod_{i=1}^{r}q_R\left(\widetilde{n}^{(s_i)}\mid\widetilde{n}^{(s_{i-1})}\right)\times\frac{1}{\widetilde{n}^{(0)}!}(1-\delta)^{\widetilde{n}^{(0)}}\left(\mu^{(j)}\right)^{-\widetilde{n}^{(j)}}\exp\left\{(1-\delta)\mu^{(j)}\right\}. b × i = 1 ∏ r q R ( n ( s i ) ∣ n ( s i − 1 ) ) × n ( 0 ) ! 1 ( 1 − δ ) n ( 0 ) ( μ ( j ) ) − n ( j ) exp { ( 1 − δ ) μ ( j ) } . which is the left hand side of (26) , completing the proof.
Justification for Backward Sampling
We prove the validity of xPGBS in this section. The proof is similar to the theoretical justification of Ancestor Sampling in Lindsten et al. [11 ] . Recall from § 3.3 that for Backwards sampling, the state at observation time t 1 t_1 t 1 x 1 = ( x 0 , x τ , … , x t ℓ ) ⊤ x_1=\left(x_0,x_\tau,\ldots,x_{t_\ell}\right)^\top x 1 = ( x 0 , x τ , … , x t ℓ ) ⊤ t ℓ t_\ell t ℓ x ℓ = ( x t ℓ − 1 + τ , x t ℓ − 1 + 2 τ , … , x t ℓ ) ⊤ x_\ell=\left(x_{t_{\ell-1}+\tau},x_{t_{\ell-1}+2\tau},\ldots,x_{t_\ell}\right)^\top x ℓ = ( x t ℓ − 1 + τ , x t ℓ − 1 + 2 τ , … , x t ℓ ) ⊤ ℓ = 2 , … , L \ell=2,\ldots,L ℓ = 2 , … , L
Following Lindsten et al. [11 ] , let us define a sequence of unnormalized densities, for ℓ = 1 , … , L \ell=1,\ldots,L ℓ = 1 , … , L
γ θ , ℓ ( x 1 : ℓ ) = p ( x 1 : ℓ , y 1 : ℓ ∣ θ ) = p 0 ( x 0 ∣ θ ) ∏ s = 0 ℓ − 1 ∏ k = 1 K p ( n s K + k ∣ μ s K + k ) ∏ s = 1 ℓ g ( y s ∣ x t s , θ ) . (27) \begin{aligned}
\gamma_{\theta,\ell}\left(x_{1:\ell}\right)=p\left(x_{1:\ell},y_{1:\ell}\mid\theta\right) \\
&=p_0\left(x_0\mid\theta\right)\prod_{s=0}^{\ell-1}\prod_{k=1}^{K}p\left(n_{sK+k}\mid\mu_{sK+k}\right)\prod_{s=1}^{\ell}g\left(y_s\mid x_{t_s},\theta\right). \tag*{(27)}
\end{aligned} γ θ , ℓ ( x 1 : ℓ ) = p ( x 1 : ℓ , y 1 : ℓ ∣ θ ) = p 0 ( x 0 ∣ θ ) s = 0 ∏ ℓ − 1 k = 1 ∏ K p ( n sK + k ∣ μ sK + k ) s = 1 ∏ ℓ g ( y s ∣ x t s , θ ) . (27) where μ s K + k \mu_{sK+k} μ sK + k n s K + k n_{sK+k} n sK + k ℓ = 0 \ell=0 ℓ = 0
γ ˉ θ , 0 ( x 0 ) : = γ θ , 0 ( x 0 ) = p 0 ( x 0 ∣ θ ) . \bar{\gamma}_{\theta,0}\left(x_0\right):=\gamma_{\theta,0}\left(x_0\right)=p_0\left(x_0\mid\theta\right). γ ˉ θ , 0 ( x 0 ) := γ θ , 0 ( x 0 ) = p 0 ( x 0 ∣ θ ) . For ℓ ≥ 1 \ell\ge1 ℓ ≥ 1
γ ˉ θ , ℓ ( x 1 : ℓ ) = γ θ , ℓ ( x 1 : ℓ ) Z θ , ℓ = p ( x 1 : ℓ ∣ y 1 : ℓ , θ ) , \bar{\gamma}_{\theta,\ell}\left(x_{1:\ell}\right)=\frac{\gamma_{\theta,\ell}\left(x_{1:\ell}\right)}{Z_{\theta,\ell}}=p\left(x_{1:\ell}\mid y_{1:\ell},\theta\right), γ ˉ θ , ℓ ( x 1 : ℓ ) = Z θ , ℓ γ θ , ℓ ( x 1 : ℓ ) = p ( x 1 : ℓ ∣ y 1 : ℓ , θ ) , where Z θ , ℓ : = ∫ γ θ , ℓ ( x 0 , x 1 : ℓ ) d x 1 : ℓ Z_{\theta,\ell}:=\int\gamma_{\theta,\ell}\left(x_0,x_{1:\ell}\right)\,dx_{1:\ell} Z θ , ℓ := ∫ γ θ , ℓ ( x 0 , x 1 : ℓ ) d x 1 : ℓ [11 ] , we define a probability function ϕ θ \phi_\theta ϕ θ
ϕ θ ( x 1 : L ( 0 : M ) , a 1 : L − 1 ( 0 : M ) , b L ) = ϕ θ ( x 1 : L ( b 1 : L ) , b 1 : L ) ϕ θ ( x 1 : L ( − b 1 : L ) , a 1 : L − 1 ( − b 2 : L ) ∣ x 1 : L ( b 1 : L ) , b 1 : L ) . \phi_\theta\left(x_{1:L}^{(0:M)},a_{1:L-1}^{(0:M)},b_L\right)=\phi_\theta\left(x_{1:L}^{(b_{1:L})},b_{1:L}\right)\phi_\theta\left(x_{1:L}^{(-b_{1:L})},a_{1:L-1}^{(-b_{2:L})}\mid x_{1:L}^{(b_{1:L})},b_{1:L}\right). ϕ θ ( x 1 : L ( 0 : M ) , a 1 : L − 1 ( 0 : M ) , b L ) = ϕ θ ( x 1 : L ( b 1 : L ) , b 1 : L ) ϕ θ ( x 1 : L ( − b 1 : L ) , a 1 : L − 1 ( − b 2 : L ) ∣ x 1 : L ( b 1 : L ) , b 1 : L ) . where the first term on the right hand side is given by
ϕ θ ( x 1 : L ( b 1 : L ) , b 1 : L ) = γ ˉ θ , ℓ ( x 1 : L ( b 1 : L ) ) ( M + 1 ) L ; \phi_\theta\left(x_{1:L}^{(b_{1:L})},b_{1:L}\right)=\frac{\bar{\gamma}_{\theta,\ell}\left(x_{1:L}^{(b_{1:L})}\right)}{(M+1)^L}; ϕ θ ( x 1 : L ( b 1 : L ) , b 1 : L ) = ( M + 1 ) L γ ˉ θ , ℓ ( x 1 : L ( b 1 : L ) ) ; that is, we first sample the path x 1 : L ( b 1 : L ) x_{1:L}^{(b_{1:L})} x 1 : L ( b 1 : L ) γ ˉ θ , ℓ ( x 1 : L ( b 1 : L ) ) \bar{\gamma}_{\theta,\ell}\left(x_{1:L}^{(b_{1:L})}\right) γ ˉ θ , ℓ ( x 1 : L ( b 1 : L ) ) b 1 : L b_{1:L} b 1 : L { 0 , 1 , … , M } \{0,1,\ldots,M\} { 0 , 1 , … , M }
Similar to the proof in Lindsten et al. [11 ] , we first show that the the xPGBS algorithm is equivalent to the following steps for sampling from the full joint probability function ϕ θ \phi_\theta ϕ θ
Procedure 1 (Instrumental reformulation of xPGBS). Given x 1 : L ′ ( b 1 : L ′ ) x_{1:L}^{\prime(b'_{1:L})} x 1 : L ′ ( b 1 : L ′ ) b 1 : L ′ ∈ { 0 , … , M } L b'_{1:L} \in\{0,\ldots,M\}^{L} b 1 : L ′ ∈ { 0 , … , M } L
Run the conditional particle filter forwards by drawing:
{ x 1 ( − b 1 ′ ) } ∼ ϕ θ ( ⋅ ∣ x 1 : L ′ ( b 1 : L ′ ) , b 1 : L ′ ) , \left\{x_1^{(-b'_1)}\right\}\sim\phi_\theta\left(\,\cdot\mid x_{1:L}^{\prime(b'_{1:L})},b'_{1:L}\right), { x 1 ( − b 1 ′ ) } ∼ ϕ θ ( ⋅ ∣ x 1 : L ′ ( b 1 : L ′ ) , b 1 : L ′ ) , and for ℓ = 2 , … , L \ell=2,\ldots,L ℓ = 2 , … , L
{ x ℓ ( − b ℓ ′ ) , a ℓ − 1 − b ℓ ′ } ∼ ϕ θ ( ⋅ ∣ x 1 : ℓ − 1 ′ ( − b 1 : ℓ − 1 ′ ) , a 1 : ℓ − 2 ( 0 : M ) , x 1 : L ′ ( b 1 : L ′ ) , b ℓ − 1 : L ′ ) . \left\{x_\ell^{(-b'_\ell)},a_{\ell-1}^{-b'_\ell}\right\}\sim\phi_\theta\left(\,\cdot\mid x_{1:\ell-1}^{\prime(-b'_{1:\ell-1})},a_{1:\ell-2}^{(0:M)},x_{1:L}^{\prime(b'_{1:L})},b'_{\ell-1:L}\right). { x ℓ ( − b ℓ ′ ) , a ℓ − 1 − b ℓ ′ } ∼ ϕ θ ( ⋅ ∣ x 1 : ℓ − 1 ′ ( − b 1 : ℓ − 1 ′ ) , a 1 : ℓ − 2 ( 0 : M ) , x 1 : L ′ ( b 1 : L ′ ) , b ℓ − 1 : L ′ ) . Draw a new reference path via backward sampling:
b L ∼ ϕ θ ( ⋅ ∣ x 1 : L ′ ( − b 1 : L ′ ) , a 1 : L − 1 ( 0 : M ) , x 1 : L ′ ( b 1 : L ′ ) ) , b_L\sim\phi_\theta\left(\,\cdot\mid x_{1:L}^{\prime(-b'_{1:L})},a_{1:L-1}^{(0:M)},x_{1:L}^{\prime(b'_{1:L})}\right), b L ∼ ϕ θ ( ⋅ ∣ x 1 : L ′ ( − b 1 : L ′ ) , a 1 : L − 1 ( 0 : M ) , x 1 : L ′ ( b 1 : L ′ ) ) , and for ℓ = L − 1 , … , 1 \ell=L-1,\ldots,1 ℓ = L − 1 , … , 1
b ℓ ≡ a ℓ b ℓ + 1 ∼ ϕ θ ( ⋅ ∣ x 1 : ℓ ′ ( − b 1 : ℓ ′ ) , a 1 : ℓ − 1 ( 0 : M ) , x 1 : ℓ ′ ( b 1 : ℓ ′ ) , x ℓ + 1 : L ( b ℓ + 1 : L ) , b ℓ + 1 : L ) . (28) b_\ell\equiv a_\ell^{b_{\ell+1}}\sim\phi_\theta\left(\,\cdot\mid x_{1:\ell}^{\prime(-b'_{1:\ell})},a_{1:\ell-1}^{(0:M)},x_{1:\ell}^{\prime(b'_{1:\ell})},x_{\ell+1:L}^{(b_{\ell+1:L})},b_{\ell+1:L}\right). \tag*{(28)} b ℓ ≡ a ℓ b ℓ + 1 ∼ ϕ θ ( ⋅ ∣ x 1 : ℓ ′ ( − b 1 : ℓ ′ ) , a 1 : ℓ − 1 ( 0 : M ) , x 1 : ℓ ′ ( b 1 : ℓ ′ ) , x ℓ + 1 : L ( b ℓ + 1 : L ) , b ℓ + 1 : L ) . (28) Lemma 5. Procedure 1 leaves the joint distribution
ϕ θ ( x 1 : L ( 0 : M ) , a 1 : L − 1 ( 0 : M ) , b L ) \phi_\theta\left(x_{1:L}^{(0:M)},a_{1:L-1}^{(0:M)},b_L\right) ϕ θ ( x 1 : L ( 0 : M ) , a 1 : L − 1 ( 0 : M ) , b L ) invariant.
Proof. By the chain rule for conditional distributions, it the straightforward to verify that step 1 samples { x 1 : L ( − b 1 : L ) , a 1 : L − 1 ( − b 2 : L ) } \{x_{1:L}^{(-b_{1:L})},a_{1:L-1}^{(-b_{2:L})}\} { x 1 : L ( − b 1 : L ) , a 1 : L − 1 ( − b 2 : L ) } ϕ θ \phi_\theta ϕ θ
{ x 1 : L ( − b 1 : L ) , a 1 : L − 1 ( − b 2 : L ) } ∼ ϕ θ ( ⋅ ∣ x 1 : L ′ ( b 1 : L ′ ) , b 1 : L ′ ) . \left\{x_{1:L}^{(-b_{1:L})},a_{1:L-1}^{(-b_{2:L})}\right\}\sim\phi_\theta\left(\,\cdot\mid x_{1:L}^{\prime(b'_{1:L})},b'_{1:L}\right). { x 1 : L ( − b 1 : L ) , a 1 : L − 1 ( − b 2 : L ) } ∼ ϕ θ ( ⋅ ∣ x 1 : L ′ ( b 1 : L ′ ) , b 1 : L ′ ) . In addition, (28) in step 2 is equivalent to sampling b ℓ b_\ell b ℓ ϕ θ \phi_\theta ϕ θ
b ℓ ≡ a ℓ b ℓ + 1 ∼ ϕ θ ( ⋅ ∣ x 1 : ℓ ′ ( − b 1 : ℓ ′ ) , a 1 : ℓ − 1 ( 0 : M ) , x 1 : ℓ ′ ( b 1 : ℓ ′ ) , x ℓ + 1 : L ( b ℓ + 1 : L ) , b ℓ + 1 : L , x ℓ + 1 : L ( − b ℓ + 1 : L ) , a ℓ : L − 1 − b ℓ + 1 : L ) . b_\ell\equiv a_\ell^{b_{\ell+1}}\sim\phi_\theta\left(\,\cdot\mid x_{1:\ell}^{\prime(-b'_{1:\ell})},a_{1:\ell-1}^{(0:M)},x_{1:\ell}^{\prime(b'_{1:\ell})},x_{\ell+1:L}^{(b_{\ell+1:L})},b_{\ell+1:L},x_{\ell+1:L}^{(-b_{\ell+1:L})},a_{\ell:L-1}^{-b_{\ell+1:L}}\right). b ℓ ≡ a ℓ b ℓ + 1 ∼ ϕ θ ( ⋅ ∣ x 1 : ℓ ′ ( − b 1 : ℓ ′ ) , a 1 : ℓ − 1 ( 0 : M ) , x 1 : ℓ ′ ( b 1 : ℓ ′ ) , x ℓ + 1 : L ( b ℓ + 1 : L ) , b ℓ + 1 : L , x ℓ + 1 : L ( − b ℓ + 1 : L ) , a ℓ : L − 1 − b ℓ + 1 : L ) . because b ℓ b_\ell b ℓ { x ℓ + 1 : L ( b ℓ + 1 : L ) , b ℓ + 1 : L } \{x_{\ell+1:L}^{(b_{\ell+1:L})},b_{\ell+1:L}\} { x ℓ + 1 : L ( b ℓ + 1 : L ) , b ℓ + 1 : L } b ℓ b_\ell b ℓ b 1 : L b_{1:L} b 1 : L ϕ θ \phi_\theta ϕ θ
Lemma 6. The algorithm of xPGBS is equivalent to Procedure 1, conditionally on x 0 ′ ( b 0 ) = x 0 ( 0 ) x_0^{\prime(b_0)}=x_0^{(0)} x 0 ′ ( b 0 ) = x 0 ( 0 ) x 1 : L ′ ( b 1 : L ′ ) = x 1 : L ( 0 ) x_{1:L}^{\prime(b'_{1:L})}=x_{1:L}^{(0)} x 1 : L ′ ( b 1 : L ′ ) = x 1 : L ( 0 ) b 0 : L ′ = ( 0 , … , 0 ) b_{0:L}'=(0,\ldots,0) b 0 : L ′ = ( 0 , … , 0 )
Proof. According to the algorithm of xPG, we have
ϕ θ ( x 1 : L ( − b 1 : L ) , a 1 : L − 1 − b 2 : L ∣ x 1 : L ( b 1 : L ) , b 1 : L ) = q ~ 0 ( x 0 ( − b 1 ) ∣ x 0 ( b 1 ) ) ∏ ℓ = 0 L − 1 ∏ k = 1 K q ( n ℓ K + k − b ℓ + 1 ∣ μ ℓ K + k ( 0 : M ) , n ℓ K + k b ℓ + 1 ) \phi_\theta\left(x_{1:L}^{(-b_{1:L})},a_{1:L-1}^{-b_{2:L}}\mid x_{1:L}^{(b_{1:L})},b_{1:L}\right)=\widetilde q_0\left(x_0^{(-b_1)}\mid x_0^{(b_1)}\right)\prod_{\ell=0}^{L-1}\prod_{k=1}^{K}q\left(n_{\ell K+k}^{-b_{\ell+1}}\mid\mu_{\ell K+k}^{(0:M)},n_{\ell K+k}^{b_{\ell+1}}\right) ϕ θ ( x 1 : L ( − b 1 : L ) , a 1 : L − 1 − b 2 : L ∣ x 1 : L ( b 1 : L ) , b 1 : L ) = q 0 ( x 0 ( − b 1 ) ∣ x 0 ( b 1 ) ) ℓ = 0 ∏ L − 1 k = 1 ∏ K q ( n ℓ K + k − b ℓ + 1 ∣ μ ℓ K + k ( 0 : M ) , n ℓ K + k b ℓ + 1 ) × { ∏ ℓ = 1 L − 1 ∏ m = 0 m ≠ b ℓ + 1 M g ( y ℓ ∣ x ℓ ( a ℓ m ) ) ∑ j = 0 M q ( y ℓ ∣ x ℓ ( j ) ) } . {}\times\left\{\prod_{\ell=1}^{L-1}\prod_{\substack{m=0\\m\ne b_{\ell+1}}}^{M}\frac{g\left(y_\ell\mid x_\ell^{(a_\ell^m)}\right)}{\sum_{j=0}^{M}q\left(y_\ell\mid x_\ell^{(j)}\right)}\right\}. × ⎩ ⎨ ⎧ ℓ = 1 ∏ L − 1 m = 0 m = b ℓ + 1 ∏ M ∑ j = 0 M q ( y ℓ ∣ x ℓ ( j ) ) g ( y ℓ ∣ x ℓ ( a ℓ m ) ) ⎭ ⎬ ⎫ . Therefore, we have
ϕ θ ( x 1 ( − b 1 ) ∣ x 1 : L ( b 1 : L ) , b 1 : L ) = q ~ 0 ( x 0 ( − b 1 ) ∣ x 0 ( b 1 ) ) ∏ k = 1 K q ( n k − b 1 ∣ μ k ( 0 : M ) , n k b 1 ) . \phi_\theta\left(x_1^{(-b_1)}\mid x_{1:L}^{(b_{1:L})},b_{1:L}\right)=\widetilde q_0\left(x_0^{(-b_1)}\mid x_0^{(b_1)}\right)\prod_{k=1}^{K}q\left(n_k^{-b_1}\mid\mu_k^{(0:M)},n_k^{b_1}\right). ϕ θ ( x 1 ( − b 1 ) ∣ x 1 : L ( b 1 : L ) , b 1 : L ) = q 0 ( x 0 ( − b 1 ) ∣ x 0 ( b 1 ) ) k = 1 ∏ K q ( n k − b 1 ∣ μ k ( 0 : M ) , n k b 1 ) . and
ϕ θ ( x ℓ ( − b ℓ ) , a ℓ − 1 − b ℓ ∣ x 1 : ℓ − 1 ( − b 1 : ℓ − 1 ) , a 1 : ℓ − 2 − b 2 : ℓ − 1 , x 1 : L ( b 1 : L ) , b 1 : L ) = ∏ k = 1 K q ( n ( ℓ − 1 ) K + k − b ℓ ∣ μ ( ℓ − 1 ) K + k ( 0 : M ) , n ( ℓ − 1 ) K + k b ℓ ) (29) \phi_\theta\left(x_\ell^{(-b_\ell)},a_{\ell-1}^{-b_\ell}\mid x_{1:\ell-1}^{(-b_{1:\ell-1})},a_{1:\ell-2}^{-b_{2:\ell-1}},x_{1:L}^{(b_{1:L})},b_{1:L}\right)=\prod_{k=1}^{K}q\left(n_{(\ell-1)K+k}^{-b_\ell}\mid\mu_{(\ell-1)K+k}^{(0:M)},n_{(\ell-1)K+k}^{b_\ell}\right) \tag*{(29)} ϕ θ ( x ℓ ( − b ℓ ) , a ℓ − 1 − b ℓ ∣ x 1 : ℓ − 1 ( − b 1 : ℓ − 1 ) , a 1 : ℓ − 2 − b 2 : ℓ − 1 , x 1 : L ( b 1 : L ) , b 1 : L ) = k = 1 ∏ K q ( n ( ℓ − 1 ) K + k − b ℓ ∣ μ ( ℓ − 1 ) K + k ( 0 : M ) , n ( ℓ − 1 ) K + k b ℓ ) (29) × ∏ m = 0 m ≠ b ℓ M g ( y ℓ − 1 ∣ x t ℓ − 1 ( a ℓ − 1 m ) ) ∑ j = 0 M g ( y ℓ − 1 ∣ x t ℓ − 1 ( j ) ) . \times\prod_{\substack{m=0\\m\ne b_{\ell}}}^{M}
\frac{g\left(y_{\ell-1}\mid x_{t_{\ell-1}}^{(a_{\ell-1}^{m})}\right)}
{\sum_{j=0}^{M}g\left(y_{\ell-1}\mid x_{t_{\ell-1}}^{(j)}\right)}. × m = 0 m = b ℓ ∏ M ∑ j = 0 M g ( y ℓ − 1 ∣ x t ℓ − 1 ( j ) ) g ( y ℓ − 1 ∣ x t ℓ − 1 ( a ℓ − 1 m ) ) . With b ℓ = 0 b_{\ell}=0 b ℓ = 0 ℓ = 1 , … , L \ell=1,\ldots,L ℓ = 1 , … , L
For the backward sampling step, for ℓ = L − 1 , … , 1 \ell=L-1,\ldots,1 ℓ = L − 1 , … , 1 a ℓ b ℓ + 1 = b ℓ a_{\ell}^{b_{\ell+1}}=b_{\ell} a ℓ b ℓ + 1 = b ℓ ℓ = 2 , … , L \ell=2,\ldots,L ℓ = 2 , … , L b 1 : ℓ − 1 = a 1 : ℓ − 1 b 2 : ℓ b_{1:\ell-1}=a_{1:\ell-1}^{b_{2:\ell}} b 1 : ℓ − 1 = a 1 : ℓ − 1 b 2 : ℓ a 1 : ℓ − 1 − b 2 : ℓ a_{1:\ell-1}^{-b_{2:\ell}} a 1 : ℓ − 1 − b 2 : ℓ a 1 − b 2 a_{1}^{-b_{2}} a 1 − b 2 a 2 − b 3 , … , a ℓ − 1 − b ℓ a_{2}^{-b_{3}},\ldots,a_{\ell-1}^{-b_{\ell}} a 2 − b 3 , … , a ℓ − 1 − b ℓ
ϕ θ ( a ℓ b ℓ + 1 ∣ x 1 : ℓ ( 0 : M ) , a 1 : ℓ − 1 ( 0 : M ) , x ℓ + 1 : L ( b ℓ + 1 : L ) , b ℓ + 1 : L ) ∝ ϕ θ ( x 1 : ℓ ( 0 : M ) , a 1 : ℓ − 1 ( 0 : M ) , x ℓ + 1 : L ( b ℓ + 1 : L ) , b ℓ : L ) = ϕ θ ( x 1 : L ( b 1 : L ) , b 1 : L ) ϕ θ ( x 1 : ℓ ( − b 1 : ℓ ) , a 1 : ℓ − 1 − b 2 : ℓ ∣ x 1 : L ( b 1 : L ) , b 1 : L ) = γ ˉ θ , L ( x 1 : L ( b 1 : L ) ) ( M + 1 ) L + 1 q ~ θ ( x 0 ( − b 1 ) ∣ x 0 ( b 1 ) ) ∏ s = 0 ℓ − 1 ∏ k = 1 K q ( n s K + k − b s + 1 ∣ μ s K + k ( 0 : M ) , n s K + k b s + 1 ) { ∏ s = 1 ℓ − 1 ∏ m = 0 m ≠ b s + 1 M g ( y s ∣ x t s ( a s m ) ) ∑ j = 0 M g ( y s ∣ x t s ( j ) ) } = γ ˉ θ , L ( x 1 : L ( b 1 : L ) ) γ θ , ℓ ( x 1 : ℓ ( b 1 : ℓ ) ) ( M + 1 ) L + 1 γ θ , ℓ ( x 1 : ℓ ( b 1 : ℓ ) ) q ~ θ ( x 0 ( − b 1 ) ∣ x 0 ( b 1 ) ) ∏ s = 0 ℓ − 1 ∏ k = 1 K q ( n s K + k − b s + 1 ∣ μ s K + k ( 0 : M ) , n s K + k b s + 1 ) × { ∏ s = 1 ℓ − 1 ∏ m = 0 m ≠ b s + 1 M g ( y s ∣ x t s ( a s m ) ) ∑ j = 0 M g ( y s ∣ x t s ( j ) ) } . (30) \begin{aligned}
\phi_{\theta}\left(a_{\ell}^{b_{\ell+1}}\mid x_{1:\ell}^{(0:M)},a_{1:\ell-1}^{(0:M)},x_{\ell+1:L}^{(b_{\ell+1:L})},b_{\ell+1:L}\right)
\propto\phi_{\theta}\left(x_{1:\ell}^{(0:M)},a_{1:\ell-1}^{(0:M)},x_{\ell+1:L}^{(b_{\ell+1:L})},b_{\ell:L}\right) \\
&=\phi_{\theta}\left(x_{1:L}^{(b_{1:L})},b_{1:L}\right)\phi_{\theta}\left(x_{1:\ell}^{(-b_{1:\ell})},a_{1:\ell-1}^{-b_{2:\ell}}\mid x_{1:L}^{(b_{1:L})},b_{1:L}\right) \\
&=\frac{\bar{\gamma}_{\theta,L}\left(x_{1:L}^{(b_{1:L})}\right)}{(M+1)^{L+1}}\widetilde{q}_{\theta}\left(x_{0}^{(-b_{1})}\mid x_{0}^{(b_{1})}\right)\prod_{s=0}^{\ell-1}\prod_{k=1}^{K}q\left(n_{sK+k}^{-b_{s+1}}\mid\mu_{sK+k}^{(0:M)},n_{sK+k}^{b_{s+1}}\right)
\left\{\prod_{s=1}^{\ell-1}\prod_{\substack{m=0\\m\ne b_{s+1}}}^{M}\frac{g\left(y_s\mid x_{t_s}^{(a_s^m)}\right)}{\sum_{j=0}^{M}g\left(y_s\mid x_{t_s}^{(j)}\right)}\right\} \\
&=\frac{\bar{\gamma}_{\theta,L}\left(x_{1:L}^{(b_{1:L})}\right)\gamma_{\theta,\ell}\left(x_{1:\ell}^{(b_{1:\ell})}\right)}
{(M+1)^{L+1}\gamma_{\theta,\ell}\left(x_{1:\ell}^{(b_{1:\ell})}\right)}
\widetilde{q}_{\theta}\left(x_{0}^{(-b_{1})}\mid x_{0}^{(b_{1})}\right)
\prod_{s=0}^{\ell-1}\prod_{k=1}^{K}q\left(n_{sK+k}^{-b_{s+1}}\mid\mu_{sK+k}^{(0:M)},n_{sK+k}^{b_{s+1}}\right) \\
&\times\left\{\prod_{s=1}^{\ell-1}\prod_{\substack{m=0\\m\ne b_{s+1}}}^{M}\frac{g\left(y_s\mid x_{t_s}^{(a_s^m)}\right)}{\sum_{j=0}^{M}g\left(y_s\mid x_{t_s}^{(j)}\right)}\right\}. \tag*{(30)}
\end{aligned} ϕ θ ( a ℓ b ℓ + 1 ∣ x 1 : ℓ ( 0 : M ) , a 1 : ℓ − 1 ( 0 : M ) , x ℓ + 1 : L ( b ℓ + 1 : L ) , b ℓ + 1 : L ) ∝ ϕ θ ( x 1 : ℓ ( 0 : M ) , a 1 : ℓ − 1 ( 0 : M ) , x ℓ + 1 : L ( b ℓ + 1 : L ) , b ℓ : L ) = ϕ θ ( x 1 : L ( b 1 : L ) , b 1 : L ) ϕ θ ( x 1 : ℓ ( − b 1 : ℓ ) , a 1 : ℓ − 1 − b 2 : ℓ ∣ x 1 : L ( b 1 : L ) , b 1 : L ) = ( M + 1 ) L + 1 γ ˉ θ , L ( x 1 : L ( b 1 : L ) ) q θ ( x 0 ( − b 1 ) ∣ x 0 ( b 1 ) ) s = 0 ∏ ℓ − 1 k = 1 ∏ K q ( n sK + k − b s + 1 ∣ μ sK + k ( 0 : M ) , n sK + k b s + 1 ) ⎩ ⎨ ⎧ s = 1 ∏ ℓ − 1 m = 0 m = b s + 1 ∏ M ∑ j = 0 M g ( y s ∣ x t s ( j ) ) g ( y s ∣ x t s ( a s m ) ) ⎭ ⎬ ⎫ = ( M + 1 ) L + 1 γ θ , ℓ ( x 1 : ℓ ( b 1 : ℓ ) ) γ ˉ θ , L ( x 1 : L ( b 1 : L ) ) γ θ , ℓ ( x 1 : ℓ ( b 1 : ℓ ) ) q θ ( x 0 ( − b 1 ) ∣ x 0 ( b 1 ) ) s = 0 ∏ ℓ − 1 k = 1 ∏ K q ( n sK + k − b s + 1 ∣ μ sK + k ( 0 : M ) , n sK + k b s + 1 ) × ⎩ ⎨ ⎧ s = 1 ∏ ℓ − 1 m = 0 m = b s + 1 ∏ M ∑ j = 0 M g ( y s ∣ x t s ( j ) ) g ( y s ∣ x t s ( a s m ) ) ⎭ ⎬ ⎫ . (30) Using (27) to expand γ θ , ℓ ( x 1 : ℓ ( b 1 : ℓ ) ) \gamma_{\theta,\ell}\left(x_{1:\ell}^{(b_{1:\ell})}\right) γ θ , ℓ ( x 1 : ℓ ( b 1 : ℓ ) ) (30) , we find that the expression above is proportional to the ratio
γ ˉ θ , L ( x 1 : L ( b 1 : L ) ) γ θ , ℓ ( x 1 : ℓ ( b 1 : ℓ ) ) \frac{\bar{\gamma}_{\theta,L}\left(x_{1:L}^{(b_{1:L})}\right)}
{\gamma_{\theta,\ell}\left(x_{1:\ell}^{(b_{1:\ell})}\right)} γ θ , ℓ ( x 1 : ℓ ( b 1 : ℓ ) ) γ ˉ θ , L ( x 1 : L ( b 1 : L ) ) multiplied by the product of three terms.
The first term accounts for both the prior density of the initial state and the proposal distribution for the remaining initial states,
q 0 ( x 0 ( 0 : M ) ) ≡ p 0 ( x 0 ( b 1 ) ) q ~ 0 ( x 0 ( − b 1 ) ) . (31) q_0\left(x_0^{(0:M)}\right)\equiv p_0\left(x_0^{(b_1)}\right)\widetilde{q}_0\left(x_0^{(-b_1)}\right). \tag*{(31)} q 0 ( x 0 ( 0 : M ) ) ≡ p 0 ( x 0 ( b 1 ) ) q 0 ( x 0 ( − b 1 ) ) . (31) The second term combines the marginal distribution of the reaction counts with the proposal distribution for the remaining reaction counts at each time step up to t = ℓ K t=\ell K t = ℓ K
∏ s = 0 ℓ − 1 ∏ k = 1 K p ( n s K + k b s + 1 ∣ μ s K + k b s + 1 ) q ~ ( n s K + k − b s + 1 ∣ μ s K + k ( 0 : M ) , n s K + k b s + 1 ) . (32) \prod_{s=0}^{\ell-1}\prod_{k=1}^{K}p\left(n_{sK+k}^{b_{s+1}}\mid\mu_{sK+k}^{b_{s+1}}\right)\widetilde{q}\left(n_{sK+k}^{-b_{s+1}}\mid\mu_{sK+k}^{(0:M)},n_{sK+k}^{b_{s+1}}\right). \tag*{(32)} s = 0 ∏ ℓ − 1 k = 1 ∏ K p ( n sK + k b s + 1 ∣ μ sK + k b s + 1 ) q ( n sK + k − b s + 1 ∣ μ sK + k ( 0 : M ) , n sK + k b s + 1 ) . (32) Finally, the third term captures the contributions from the observation likelihoods and path resampling probabilities up to t = ℓ K t=\ell K t = ℓ K
g ( y ℓ ∣ x t ℓ ( b ℓ ) ) ∏ s = 1 ℓ − 1 { g ( y s ∣ x t s ( b s ) ) ∏ m = 0 m ≠ b s + 1 M g ( y s ∣ x t s ( a s m ) ) ∑ j = 0 M g ( y s ∣ x t s ( j ) ) } = g ( y ℓ ∣ x t ℓ ( b ℓ ) ) ∏ s = 1 ℓ − 1 { ( ∑ j = 0 M g ( y s ∣ x t s ( j ) ) ) − M ∏ m = 0 M g ( y s ∣ x t s ( a s m ) ) } . (33) \begin{aligned}
g\left(y_{\ell}\mid x_{t_{\ell}}^{(b_{\ell})}\right)\prod_{s=1}^{\ell-1}\left\{g\left(y_s\mid x_{t_s}^{(b_s)}\right)\prod_{\substack{m=0\\m\ne b_{s+1}}}^{M}\frac{g\left(y_s\mid x_{t_s}^{(a_s^m)}\right)}{\sum_{j=0}^{M}g\left(y_s\mid x_{t_s}^{(j)}\right)}\right\} \\
&=g\left(y_{\ell}\mid x_{t_{\ell}}^{(b_{\ell})}\right)\prod_{s=1}^{\ell-1}\left\{\left(\sum_{j=0}^{M}g\left(y_s\mid x_{t_s}^{(j)}\right)\right)^{-M}\prod_{m=0}^{M}g\left(y_s\mid x_{t_s}^{(a_s^m)}\right)\right\}. \tag*{(33)}
\end{aligned} g ( y ℓ ∣ x t ℓ ( b ℓ ) ) s = 1 ∏ ℓ − 1 ⎩ ⎨ ⎧ g ( y s ∣ x t s ( b s ) ) m = 0 m = b s + 1 ∏ M ∑ j = 0 M g ( y s ∣ x t s ( j ) ) g ( y s ∣ x t s ( a s m ) ) ⎭ ⎬ ⎫ = g ( y ℓ ∣ x t ℓ ( b ℓ ) ) s = 1 ∏ ℓ − 1 ⎩ ⎨ ⎧ ( j = 0 ∑ M g ( y s ∣ x t s ( j ) ) ) − M m = 0 ∏ M g ( y s ∣ x t s ( a s m ) ) ⎭ ⎬ ⎫ . (33) since a s b s + 1 = b s a_s^{b_{s+1}}=b_s a s b s + 1 = b s σ ( 0 : M ) \sigma(0:M) σ ( 0 : M ) b ℓ b_{\ell} b ℓ
Sampling from (28) is equivalent to selecting b ℓ b_\ell b ℓ b ℓ ( 0 : M ) b_\ell^{(0:M)} b ℓ ( 0 : M ) b ℓ ( i ) b_\ell^{(i)} b ℓ ( i ) b ℓ = b ℓ ( i ) b_\ell=b_\ell^{(i)} b ℓ = b ℓ ( i ) b ℓ ( i ) b_\ell^{(i)} b ℓ ( i )
g ( y ℓ ∣ x t ℓ ( b 1 : ℓ ( i ) ) ) γ ˉ θ , L ( x 1 : ℓ ( b 1 : ℓ ( i ) ) , x ℓ + 1 : L ( b ℓ + 1 : L ) ) γ θ , ℓ ( x 1 : ℓ ( b 1 : ℓ ( i ) ) ) (34) g\left(y_\ell\mid x_{t_\ell}^{(b_{1:\ell}^{(i)})}\right)
\frac{\bar{\gamma}_{\theta,L}\left(x_{1:\ell}^{(b_{1:\ell}^{(i)})},x_{\ell+1:L}^{(b_{\ell+1:L})}\right)}
{\gamma_{\theta,\ell}\left(x_{1:\ell}^{(b_{1:\ell}^{(i)})}\right)} \tag*{(34)} g ( y ℓ ∣ x t ℓ ( b 1 : ℓ ( i ) ) ) γ θ , ℓ ( x 1 : ℓ ( b 1 : ℓ ( i ) ) ) γ ˉ θ , L ( x 1 : ℓ ( b 1 : ℓ ( i ) ) , x ℓ + 1 : L ( b ℓ + 1 : L ) ) (34) where b 1 : ℓ ( i ) b_{1:\ell}^{(i)} b 1 : ℓ ( i ) b ℓ = i b_\ell=i b ℓ = i
Fixing b s = b ~ s b_s=\tilde{b}_s b s = b ~ s s = ℓ + 1 , … , L s=\ell+1,\ldots,L s = ℓ + 1 , … , L b ℓ ( i ) = i b_\ell^{(i)}=i b ℓ ( i ) = i w ℓ i = g ( y ℓ ∣ x t ℓ ( i ) ) w_\ell^i=g(y_\ell\mid x_{t_\ell}^{(i)}) w ℓ i = g ( y ℓ ∣ x t ℓ ( i ) )
g ( y ℓ ∣ x t ℓ ( i ) ) p ( x ~ t ℓ + τ ( b ℓ + 1 ) ∣ x t ℓ ( i ) , θ ) = w ℓ i ⋅ p ( x t ℓ + τ ( b ~ ℓ + 1 ) ∣ x t ℓ ( i ) , θ ) , g\left(y_\ell\mid x_{t_\ell}^{(i)}\right)p\left(\tilde{x}_{t_{\ell+\tau}}^{(b_{\ell+1})}\mid x_{t_\ell}^{(i)},\theta\right)
=w_\ell^i\cdot p\left(x_{t_{\ell+\tau}}^{(\tilde{b}_{\ell+1})}\mid x_{t_\ell}^{(i)},\theta\right), g ( y ℓ ∣ x t ℓ ( i ) ) p ( x ~ t ℓ + τ ( b ℓ + 1 ) ∣ x t ℓ ( i ) , θ ) = w ℓ i ⋅ p ( x t ℓ + τ ( b ~ ℓ + 1 ) ∣ x t ℓ ( i ) , θ ) , which corresponds to the backward sampling step described in Secton 3.3.
Justification of the tuning recommendations
Rationale for the assumptions
Assumption 1 states that, whether or not a particle has the reference path in its ancestry does not affect its average weight just prior to the current resampling step. The mixing of the Markov chain { X t } t ∈ [ 0 , T ] \{X_t\}_{t\in[0,T]} { X t } t ∈ [ 0 , T ] s ≪ t s\ll t s ≪ t t t t s s s t t t δ \delta δ
The intuition behind Assumption 2 is as follows. After a large number of MCMC iterations, the full reference path would achieve its stationary distribution. Since the observations are fixed, it is therefore natural to assume that the average of the particle weights also has its own stationary distribution which does not change substantially with the number of proposed particles M M M M M M
Assumption 2 also states that the number of particles from/not from the reference path at the previous observation time (after resampling) is independent of the current weight of the reference path and the average weight of all other particles at the current observation time. Again, we do not expect the assumption to hold precisely, but rather that any dependence is weak enough to be ignored. Suppose the latent process { X t ℓ } ℓ = 1 L \{X_{t_\ell}\}_{\ell=1}^L { X t ℓ } ℓ = 1 L { Y ℓ } ℓ = 1 L \{Y_\ell\}_{\ell=1}^L { Y ℓ } ℓ = 1 L { ( X t ℓ , Y ℓ ) } ℓ = 1 L \{(X_{t_\ell},Y_\ell)\}_{\ell=1}^L {( X t ℓ , Y ℓ ) } ℓ = 1 L { ( X t ℓ ( 1 : M ) , Y ℓ ) } ℓ = 1 L \{(X_{t_\ell}^{(1:M)},Y_\ell)\}_{\ell=1}^L {( X t ℓ ( 1 : M ) , Y ℓ ) } ℓ = 1 L W ℓ ( 0 ) = g ( Y ℓ ∣ X t ℓ ( 0 ) , θ ) W_\ell^{(0)}=g(Y_\ell\mid X_{t_\ell}^{(0)},\theta) W ℓ ( 0 ) = g ( Y ℓ ∣ X t ℓ ( 0 ) , θ ) W ‾ ℓ = 1 M ∑ i = 1 M g ( Y ℓ ∣ X t ℓ ( i ) , θ ) \overline{W}_\ell=\frac{1}{M}\sum_{i=1}^M g(Y_\ell\mid X_{t_\ell}^{(i)},\theta) W ℓ = M 1 ∑ i = 1 M g ( Y ℓ ∣ X t ℓ ( i ) , θ ) { W ℓ ( 0 ) } ℓ = 1 L \{W_\ell^{(0)}\}_{\ell=1}^L { W ℓ ( 0 ) } ℓ = 1 L { W ‾ ℓ } ℓ = 1 L \{\overline{W}_\ell\}_{\ell=1}^L { W ℓ } ℓ = 1 L y 1 : L y_{1:L} y 1 : L
Proof of Lemma 2
At the first observation time t 1 t_1 t 1
M 1 G ∼ Bin ( M , p 1 ( G ) ) , M_1^G\sim\operatorname{Bin}(M,p_1(\mathcal{G})), M 1 G ∼ Bin ( M , p 1 ( G )) , where
p 1 ( G ) = M W ‾ 1 M W ‾ 1 + W 1 ( 0 ) = P 1 ( M ) . p_1(\mathcal{G})=\frac{M\overline{W}_1}{M\overline{W}_1+W_1^{(0)}}=P_1^{(M)}. p 1 ( G ) = M W 1 + W 1 ( 0 ) M W 1 = P 1 ( M ) . At observation time t ℓ t_\ell t ℓ ℓ = 2 , … , L \ell=2,\ldots,L ℓ = 2 , … , L M ℓ − 1 G M_{\ell-1}^G M ℓ − 1 G p ℓ ( G ∣ M ℓ − 1 G ) p_\ell(\mathcal{G}\mid M_{\ell-1}^G) p ℓ ( G ∣ M ℓ − 1 G ) M + 1 M+1 M + 1 t ℓ t_\ell t ℓ M ℓ − 1 G M_{\ell-1}^G M ℓ − 1 G
M ℓ G ∼ Bin ( M , p ℓ ( G ∣ M ℓ − 1 G ) ) , M_\ell^G\sim\operatorname{Bin}\left(M,p_\ell(\mathcal{G}\mid M_{\ell-1}^G)\right), M ℓ G ∼ Bin ( M , p ℓ ( G ∣ M ℓ − 1 G ) ) , where p ℓ ( G ∣ M ℓ − 1 G ) = M ℓ − 1 G W ‾ ℓ M W ‾ ℓ + W ℓ ( 0 ) = M ℓ − 1 G W ‾ ℓ M W ‾ ℓ + W ℓ ( 0 ) p_{\ell}(G \mid M_{\ell-1}^{G}) = \frac{M_{\ell-1}^{G}\overline{W}_{\ell}}{M\overline{W}_{\ell} + W_{\ell}^{(0)}} = \frac{M_{\ell-1}^{G}\overline{W}_{\ell}}{M\overline{W}_{\ell} + W_{\ell}^{(0)}} p ℓ ( G ∣ M ℓ − 1 G ) = M W ℓ + W ℓ ( 0 ) M ℓ − 1 G W ℓ = M W ℓ + W ℓ ( 0 ) M ℓ − 1 G W ℓ
where the last equality follows from Assumption 1. We prove the result by induction on ℓ \ell ℓ ℓ = 2 \ell= 2 ℓ = 2
E M 1 G [ p 2 ( G ) ] = E M 1 G [ p 2 ( G ∣ M 1 G ) ] = E [ M 1 G ] W ‾ 2 M W ‾ 2 + W 2 ( 0 ) = P 1 ( M ) M W ‾ 2 M W ‾ 2 + W 2 ( 0 ) = P 1 ( M ) P 2 ( M ) . \mathbb{E}_{M_{1}^{G}}[p_{2}(G)] = \mathbb{E}_{M_{1}^{G}}\left[p_{2}(G \mid M_{1}^{G})\right] = \frac{\mathbb{E}[M_{1}^{G}]\overline{W}_{2}}{M\overline{W}_{2} + W_{2}^{(0)}} = \frac{P_{1}^{(M)}M\overline{W}_{2}}{M\overline{W}_{2} + W_{2}^{(0)}} = P_{1}^{(M)}P_{2}^{(M)}. E M 1 G [ p 2 ( G )] = E M 1 G [ p 2 ( G ∣ M 1 G ) ] = M W 2 + W 2 ( 0 ) E [ M 1 G ] W 2 = M W 2 + W 2 ( 0 ) P 1 ( M ) M W 2 = P 1 ( M ) P 2 ( M ) . Suppose now that it holds for some 2 ≤ ℓ ≤ L − 1 2 \le\ell\le L - 1 2 ≤ ℓ ≤ L − 1 ℓ + 1 \ell+ 1 ℓ + 1 E M 1 : ℓ G [ p ℓ + 1 ( G ) ] \mathbb{E}_{M_{1:\ell}^{G}}[p_{\ell+1}(G)] E M 1 : ℓ G [ p ℓ + 1 ( G )]
E M ℓ G [ E M 1 : ( ℓ − 1 ) G [ p ℓ + 1 ( G ∣ M ℓ G ) ] ] = E M ℓ G [ E M 1 : ( ℓ − 1 ) G [ M ℓ G W ‾ ℓ + 1 M W ‾ ℓ + 1 + W ℓ + 1 ( 0 ) ] ] = M E M 1 : ( ℓ − 1 ) G [ p ℓ ( G ) ] W ‾ ℓ + 1 M W ‾ ℓ + 1 + W ℓ + 1 ( 0 ) = { ∏ k = 1 ℓ P k ( M ) } P ℓ + 1 ( M ) . \begin{aligned}
\mathbb{E}_{M_{\ell}^{G}}\left[\mathbb{E}_{M_{1:(\ell-1)}^{G}}\left[p_{\ell+1}(G \mid M_{\ell}^{G})\right]\right] \\
&= \mathbb{E}_{M_{\ell}^{G}}\left[\mathbb{E}_{M_{1:(\ell-1)}^{G}}\left[M_{\ell}^{G}\frac{\overline{W}_{\ell+1}}{M\overline{W}_{\ell+1} + W_{\ell+1}^{(0)}}\right]\right] \\
&= M\mathbb{E}_{M_{1:(\ell-1)}^{G}}\left[p_{\ell}(G)\right]\frac{\overline{W}_{\ell+1}}{M\overline{W}_{\ell+1} + W_{\ell+1}^{(0)}} \\
&= \left\{\prod_{k=1}^{\ell}P_{k}^{(M)}\right\}P_{\ell+1}^{(M)}.
\end{aligned} E M ℓ G [ E M 1 : ( ℓ − 1 ) G [ p ℓ + 1 ( G ∣ M ℓ G ) ] ] = E M ℓ G [ E M 1 : ( ℓ − 1 ) G [ M ℓ G M W ℓ + 1 + W ℓ + 1 ( 0 ) W ℓ + 1 ] ] = M E M 1 : ( ℓ − 1 ) G [ p ℓ ( G ) ] M W ℓ + 1 + W ℓ + 1 ( 0 ) W ℓ + 1 = { k = 1 ∏ ℓ P k ( M ) } P ℓ + 1 ( M ) . Therefore, by induction, the claim holds for all ℓ ∈ { 2 , … , L } \ell\in\{2,\ldots,L\} ℓ ∈ { 2 , … , L }
Proof of Lemma 3
To simplify notation, we abbreviate P ( L , M ) P^{(L,M)} P ( L , M ) E ( L , M ) E^{(L,M)} E ( L , M ) P ( M ) P^{(M)} P ( M ) E ( M ) E^{(M)} E ( M )
Proposition 1. For x ≥ 0 x \ge0 x ≥ 0 y ≥ 0 y \ge0 y ≥ 0
∣ e − x − e − y ∣ ≤ ∣ x − y ∣ . \left|e^{-x} - e^{-y}\right| \le|x-y|. e − x − e − y ≤ ∣ x − y ∣. Proof. Let f ( x ) = e − x f(x)=e^{-x} f ( x ) = e − x
∣ f ′ ( x ) ∣ = e − x ≤ 1 , x ≥ 0. |f'(x)| = e^{-x} \le1,\qquad x \ge0. ∣ f ′ ( x ) ∣ = e − x ≤ 1 , x ≥ 0. By the mean value theorem,
f ( x ) − f ( y ) = ( x − y ) f ′ ( ξ ) f(x)-f(y)=(x-y)f'(\xi) f ( x ) − f ( y ) = ( x − y ) f ′ ( ξ ) for some ξ \xi ξ x x x
∣ f ( x ) − f ( y ) ∣ = ∣ x − y ∣ ∣ f ′ ( ξ ) ∣ ≤ ∣ x − y ∣ , |f(x)-f(y)|=|x-y||f'(\xi)|\le|x-y|, ∣ f ( x ) − f ( y ) ∣ = ∣ x − y ∣∣ f ′ ( ξ ) ∣ ≤ ∣ x − y ∣ , since ξ ≥ 0 \xi\ge0 ξ ≥ 0
Proposition 2. For x ≥ 0 x \ge0 x ≥ 0
1 − x ≤ e − x ≤ 1 1 + x . 1-x \le e^{-x} \le\frac{1}{1+x}. 1 − x ≤ e − x ≤ 1 + x 1 . Proof. For the first inequality, apply Taylor’s theorem to f ( x ) = e − x f(x)=e^{-x} f ( x ) = e − x
f ( x ) = f ( 0 ) + x f ′ ( 0 ) + x 2 2 f ′ ′ ( ξ ) f(x)=f(0)+xf'(0)+\frac{x^{2}}{2}f''(\xi) f ( x ) = f ( 0 ) + x f ′ ( 0 ) + 2 x 2 f ′′ ( ξ ) for some ξ ∈ ( 0 , x ) \xi\in(0,x) ξ ∈ ( 0 , x ) f ( 0 ) = 1 f(0)=1 f ( 0 ) = 1 f ′ ( 0 ) = − 1 f'(0)=-1 f ′ ( 0 ) = − 1 f ′ ′ ( ξ ) = e − ξ ≥ 0 f''(\xi)=e^{-\xi}\ge0 f ′′ ( ξ ) = e − ξ ≥ 0
e − x = 1 − x + x 2 2 e − ξ ≥ 1 − x . e^{-x}=1-x+\frac{x^{2}}{2}e^{-\xi}\ge1-x. e − x = 1 − x + 2 x 2 e − ξ ≥ 1 − x . For the second inequality, apply the same argument to g ( x ) = e x g(x)=e^{x} g ( x ) = e x
e x = 1 + x + x 2 2 e η e^{x}=1+x+\frac{x^{2}}{2}e^{\eta} e x = 1 + x + 2 x 2 e η for some η ∈ ( 0 , x ) \eta\in(0,x) η ∈ ( 0 , x )
e x ≥ 1 + x . e^x \ge1+x. e x ≥ 1 + x . Taking reciprocals yields
e − x ≤ 1 1 + x . e^{-x} \le\frac{1}{1+x}. e − x ≤ 1 + x 1 . This proves the result. □
Lemma 7.
E [ ∣ exp ( − 1 M ∑ ℓ = 1 L R ℓ ) − exp ( − 1 M ∑ ℓ = 1 L E ℓ ) ∣ ] ≤ 1 M Var ( ∑ ℓ = 1 L R ℓ ) 1 / 2 . \mathbb{E}\left[\left|\exp\left(-\frac{1}{M}\sum_{\ell=1}^{L}R_\ell\right)-\exp\left(-\frac{1}{M}\sum_{\ell=1}^{L}E_\ell\right)\right|\right]\leq\frac{1}{M}\operatorname{Var}\left(\sum_{\ell=1}^{L}R_\ell\right)^{1/2}. E [ exp ( − M 1 ℓ = 1 ∑ L R ℓ ) − exp ( − M 1 ℓ = 1 ∑ L E ℓ ) ] ≤ M 1 Var ( ℓ = 1 ∑ L R ℓ ) 1/2 . Proof. From Proposition 1,
∣ exp ( − 1 M ∑ ℓ = 1 L R ℓ ) − exp ( − 1 M ∑ ℓ = 1 L E ℓ ) ∣ ≤ 1 M ∣ ∑ ℓ = 1 L R ℓ − ∑ ℓ = 1 L E ℓ ∣ . \left|\exp\left(-\frac{1}{M}\sum_{\ell=1}^{L}R_\ell\right)-\exp\left(-\frac{1}{M}\sum_{\ell=1}^{L}E_\ell\right)\right|\leq\frac{1}{M}\left|\sum_{\ell=1}^{L}R_\ell-\sum_{\ell=1}^{L}E_\ell\right|. exp ( − M 1 ℓ = 1 ∑ L R ℓ ) − exp ( − M 1 ℓ = 1 ∑ L E ℓ ) ≤ M 1 ℓ = 1 ∑ L R ℓ − ℓ = 1 ∑ L E ℓ . So, by Jensen’s inequality,
E [ ∣ exp ( − 1 M ∑ ℓ = 1 L R ℓ ) − exp ( − 1 M ∑ ℓ = 1 L E ℓ ) ∣ ] ≤ 1 M E [ ∣ ∑ ℓ = 1 L R ℓ − ∑ ℓ = 1 L E ℓ ∣ ] ≤ 1 M E [ ( ∑ ℓ = 1 L R ℓ − ∑ ℓ = 1 L E ℓ ) 2 ] 1 / 2 . \begin{aligned}
\mathbb{E}\left[\left|\exp\left(-\frac{1}{M}\sum_{\ell=1}^{L}R_\ell\right)-\exp\left(-\frac{1}{M}\sum_{\ell=1}^{L}E_\ell\right)\right|\right]\leq\frac{1}{M}\mathbb{E}\left[\left|\sum_{\ell=1}^{L}R_\ell-\sum_{\ell=1}^{L}E_\ell\right|\right] \\
&\leq\frac{1}{M}\mathbb{E}\left[\left(\sum_{\ell=1}^{L}R_\ell-\sum_{\ell=1}^{L}E_\ell\right)^2\right]^{1/2}.
\end{aligned} E [ exp ( − M 1 ℓ = 1 ∑ L R ℓ ) − exp ( − M 1 ℓ = 1 ∑ L E ℓ ) ] ≤ M 1 E [ ℓ = 1 ∑ L R ℓ − ℓ = 1 ∑ L E ℓ ] ≤ M 1 E ( ℓ = 1 ∑ L R ℓ − ℓ = 1 ∑ L E ℓ ) 2 1/2 . Since E ℓ = E [ R ℓ ] E_\ell=\mathbb{E}[R_\ell] E ℓ = E [ R ℓ ]
1 M Var ( ∑ ℓ = 1 L R ℓ ) 1 / 2 , \frac{1}{M}\operatorname{Var}\left(\sum_{\ell=1}^{L}R_\ell\right)^{1/2}, M 1 Var ( ℓ = 1 ∑ L R ℓ ) 1/2 , which proves the result. □
Lemma 8.
∣ P ( M ) − exp ( − 1 M ∑ ℓ = 1 L R ℓ ) ∣ ≤ 1 M 2 ∑ ℓ = 1 L { R ℓ } 2 \left|P^{(M)}-\exp\left(-\frac{1}{M}\sum_{\ell=1}^{L}R_\ell\right)\right|\leq\frac{1}{M^2}\sum_{\ell=1}^{L}\{R_\ell\}^2 P ( M ) − exp ( − M 1 ℓ = 1 ∑ L R ℓ ) ≤ M 2 1 ℓ = 1 ∑ L { R ℓ } 2 Proof. Write
P ( M ) = ∏ ℓ = 1 L a ℓ , where a ℓ : = P ℓ ( M ) and 0 ≤ a ℓ ≤ 1. P^{(M)}=\prod_{\ell=1}^{L}a_\ell,\qquad\text{where }a_\ell:=P_\ell^{(M)}\text{ and }0\leq a_\ell\leq1. P ( M ) = ℓ = 1 ∏ L a ℓ , where a ℓ := P ℓ ( M ) and 0 ≤ a ℓ ≤ 1. Write
exp ( − 1 M ∑ ℓ = 1 L R ℓ ) = ∏ ℓ = 1 L b ℓ , where b ℓ : = exp ( − R ℓ M ) and 0 ≤ b ℓ ≤ 1. \exp\left(-\frac{1}{M}\sum_{\ell=1}^{L}R_\ell\right)=\prod_{\ell=1}^{L}b_\ell,\qquad\text{where }b_\ell:=\exp\left(-\frac{R_\ell}{M}\right)\text{ and }0\leq b_\ell\leq1. exp ( − M 1 ℓ = 1 ∑ L R ℓ ) = ℓ = 1 ∏ L b ℓ , where b ℓ := exp ( − M R ℓ ) and 0 ≤ b ℓ ≤ 1. Then, using the convention that the product over an empty set is 1,
∏ ℓ = 1 L a ℓ − ∏ ℓ = 1 L b ℓ = ∑ j = 1 L ( ∏ ℓ = 1 j − 1 a ℓ ) ( a j − b j ) ( ∏ ℓ = j + 1 L b ℓ ) . \prod_{\ell=1}^{L}a_\ell-\prod_{\ell=1}^{L}b_\ell=\sum_{j=1}^{L}\left(\prod_{\ell=1}^{j-1}a_\ell\right)(a_j-b_j)\left(\prod_{\ell=j+1}^{L}b_\ell\right). ℓ = 1 ∏ L a ℓ − ℓ = 1 ∏ L b ℓ = j = 1 ∑ L ( ℓ = 1 ∏ j − 1 a ℓ ) ( a j − b j ) ℓ = j + 1 ∏ L b ℓ . So
∣ ∏ ℓ = 1 L a ℓ − ∏ ℓ = 1 L b ℓ ∣ ≤ ∑ j = 1 L ( ∏ ℓ = 1 j − 1 a ℓ ) ( a j − b j ) ( ∏ ℓ = j + 1 L b ℓ ) ≤ ∑ j = 1 L ∣ a j − b j ∣ , \left|\prod_{\ell=1}^{L}a_\ell-\prod_{\ell=1}^{L}b_\ell\right|\leq\sum_{j=1}^{L}\left(\prod_{\ell=1}^{j-1}a_\ell\right)(a_j-b_j)\left(\prod_{\ell=j+1}^{L}b_\ell\right)\leq\sum_{j=1}^{L}|a_j-b_j|, ℓ = 1 ∏ L a ℓ − ℓ = 1 ∏ L b ℓ ≤ j = 1 ∑ L ( ℓ = 1 ∏ j − 1 a ℓ ) ( a j − b j ) ℓ = j + 1 ∏ L b ℓ ≤ j = 1 ∑ L ∣ a j − b j ∣ , because all a j a_j a j b j b_j b j
From Proposition 2,
∣ 1 1 + R ℓ / M − exp ( − 1 M R ℓ ) ∣ ≤ ∣ 1 1 + R ℓ / M − ( 1 − 1 M R ℓ ) ∣ = R ℓ 2 / M 2 1 + R ℓ / M ≤ 1 M 2 R ℓ 2 , \left|\frac{1}{1+R_{\ell}/M}-\exp\left(-\frac{1}{M}R_{\ell}\right)\right|
\leq\left|\frac{1}{1+R_{\ell}/M}-\left(1-\frac{1}{M}R_{\ell}\right)\right|
= \frac{R_{\ell}^{2}/M^{2}}{1+R_{\ell}/M}
\leq\frac{1}{M^{2}}R_{\ell}^{2}, 1 + R ℓ / M 1 − exp ( − M 1 R ℓ ) ≤ 1 + R ℓ / M 1 − ( 1 − M 1 R ℓ ) = 1 + R ℓ / M R ℓ 2 / M 2 ≤ M 2 1 R ℓ 2 , Therefore, we have
∣ P ( M ) − exp ( − 1 M ∑ ℓ = 1 L R ℓ ) ∣ ≤ ∑ ℓ = 1 L ∣ P ℓ ( M ) − exp ( − R ℓ M ) ∣ ≤ 1 M 2 ∑ ℓ = 1 L { R ℓ } 2 . \left|P^{(M)}-\exp\left(-\frac{1}{M}\sum_{\ell=1}^{L}R_{\ell}\right)\right|
\leq\sum_{\ell=1}^{L}\left|P_{\ell}^{(M)}-\exp\left(-\frac{R_{\ell}}{M}\right)\right|
\leq\frac{1}{M^{2}}\sum_{\ell=1}^{L}\{R_{\ell}\}^{2}. P ( M ) − exp ( − M 1 ℓ = 1 ∑ L R ℓ ) ≤ ℓ = 1 ∑ L P ℓ ( M ) − exp ( − M R ℓ ) ≤ M 2 1 ℓ = 1 ∑ L { R ℓ } 2 . □ \square □
We now combine Lemmas 7 and 8 to prove Lemma 3. First define R ( M ) : = exp ( − 1 M ∑ ℓ = 1 L R ℓ ) R^{(M)}:=\exp\left(-\frac{1}{M}\sum_{\ell=1}^{L}R_{\ell}\right) R ( M ) := exp ( − M 1 ∑ ℓ = 1 L R ℓ )
E [ ∣ P ( M ) − E ( M ) ∣ ] ≤ E [ ∣ P ( M ) − R ( M ) ∣ ] + E [ ∣ R ( M ) − E ( M ) ∣ ] ≤ 1 M 2 E [ ∑ ℓ = 0 L R ℓ 2 ] + 1 M Var [ ∑ ℓ = 0 L R ℓ ] 1 2 . \begin{aligned}
\mathbb{E}\left[\left|P^{(M)}-E^{(M)}\right|\right]
\leq\mathbb{E}\left[\left|P^{(M)}-R^{(M)}\right|\right]
+\mathbb{E}\left[\left|R^{(M)}-E^{(M)}\right|\right] \\
&\leq\frac{1}{M^{2}}\mathbb{E}\left[\sum_{\ell=0}^{L}R_{\ell}^{2}\right]
+\frac{1}{M}\operatorname{Var}\left[\sum_{\ell=0}^{L}R_{\ell}\right]^{\frac{1}{2}}.
\end{aligned} E [ P ( M ) − E ( M ) ] ≤ E [ P ( M ) − R ( M ) ] + E [ R ( M ) − E ( M ) ] ≤ M 2 1 E [ ℓ = 0 ∑ L R ℓ 2 ] + M 1 Var [ ℓ = 0 ∑ L R ℓ ] 2 1 . Inferred numbers infectious over time
Figure 6. Plots of the numbers infectious over time in four regions in the SEIR model after running xPG for 10 6 10^{6} 1 0 6 β = 0.6 \beta= 0.6 β = 0.6 η = 0.6 \eta= 0.6 η = 0.6 γ = 0.2 \gamma= 0.2 γ = 0.2
Figure 7. Plots of the numbers infectious over time in four regions in the SEIR model after running xPGBS for 10 6 10^{6} 1 0 6 β = 0.6 \beta= 0.6 β = 0.6 η = 0.6 \eta= 0.6 η = 0.6 γ = 0.2 \gamma= 0.2 γ = 0.2
Posterior distributions of model parameters
Figure 8. Density estimates of the posterior distributions of β \beta β η \eta η γ \gamma γ M = 30 M = 30 M = 30 δ = 0.02 \delta= 0.02 δ = 0.02 M = 10 M = 10 M = 10 δ = 0.02 \delta= 0.02 δ = 0.02 2 × 10 6 2 \times10^{6} 2 × 1 0 6
Tuning diagnostics
As is illustrated in § 4.3 , Figure 9 a and Figure 9 b show that the value of M M M Figure 9 c and Figure 9 d, together with Figure 9 e and Figure 9 f, indicate that the suggested value of M M M M M M
Figure 9. Contour plots of sample efficiency for xPG applied to simulated data from four coupled SIR models (left column) and four coupled SEIR models (right column), using 10 6 10^{6} 1 0 6 t 1 t_{1} t 1 exp ( − 1 ) \exp(-1) exp ( − 1 )
Despite the overestimation of the optimal value of M M M Figure 10 and Figure 11 .
Figure 10. Sample efficiency of the worst mixing state over all regions and time points (the number of susceptible individuals) versus the acceptance rate at the first observation time for xPG applied to the four-coupled SIR model, based on 10 6 10^6 1 0 6
Figure 11. Sample efficiency of the worst mixing state over all regions and time points versus the acceptance rate at the first observation time for xPG applied to the four-coupled SEIR model, based on 10 6 10^6 1 0 6
Supplementary inference results
Figure 12. Contour plots of sample efficiency for xPG (left column) and xPGBS (right column) applied to simulated data from four coupled SEIR models, using 10 6 10^{6} 1 0 6
