# Bayesian Computation Methods for Inference in Stochastic Kinetic Models

Authors: Eugenia Koblents; Inés P. Mariño; Joaquín Míguez
Journal: Complexity
 (2019)

---

## Abstract

In this paper we investigate Monte Carlo methods for the approximation of the posterior probability distributions in stochastic kinetic models (SKMs). SKMs are multivariate Markov jump processes that model the interactions among species in biological systems according to a set of usually unknown parameters. The tracking of the species populations together with the estimation of the interaction parameters is a Bayesian inference problem for which Markov chain Monte Carlo (MCMC) methods have been a typical computational tool. Specifically, the particle MCMC (pMCMC) method has been shown to be effective, while computationally demanding method applicable to this problem. Recently, it has been shown that an alternative approach to Bayesian computation, namely, the class of adaptive importance samplers, may be more efficient than classical MCMC-like schemes, at least for certain applications. For example, the nonlinear population Monte Carlo (NPMC) algorithm has yielded promising results with a low dimensional SKM (the classical predator-prey model). In this paper we explore the application of both pMCMC and NPMC to analyze complex autoregulatory feedback networks modelled by SKMs. We demonstrate numerically how the populations of the relevant species in the network can be tracked and their interaction rates estimated, even in scenarios with partial observations. NPMC schemes attain an appealing trade-off between accuracy and computational cost that can make them advantageous in many practical applications.


---

## Body

## 1. Introduction

Stochastic kinetic models (SKMs) are multivariate systems that model interactions among species in ecological, biological, and chemical problems, according to a set of unknown rate parameters [1]. The aim of this paper is to study computational methods for the approximation of the posterior distribution of the rate parameters and the populations of all species, provided a set of discrete, noisy observations is available. This inference problem has been traditionally addressed in a Bayesian framework using Markov chain Monte Carlo (MCMC) schemes [1–4]. In [5], a particle MCMC (pMCMC) method [6] has been applied to Lotka-Volterra and prokaryotic autoregulatory models. The pMCMC technique relies on a sequential Monte Carlo (SMC) approximation of the posterior distribution of the populations to compute the Metropolis-Hastings (MH) acceptance ratio.However, MCMC methods in general and pMCMC in particular suffer from a number of problems. The convergence of the Markov chain is hard to assess and the final set of samples presents correlations which can greatly reduce its efficiency. Besides, MCMC methods do not (easily) allow for parallel implementation and turn out to be computationally intensive, even for relatively simple models. A classical technique to reduce the complexity of the existing MCMC methods when applied to SKMs is the so-called diffusion approximation of the underlying stochastic jump process [7]. However, the diffusion approximation breaks down in low-concentration scenarios [5]. Other techniques, including linear noise approximations [8–11] and Gaussian-process estimators of the likelihood [12], have also been proposed in the past few years to reduce the computational cost. The parameters of the MCMC proposal are also hard to choose and they determine the performance of the algorithm.An appealing alternative to MCMC methods is the class of the population Monte Carlo (PMC) algorithms [13]. PMC is an adaptive importance sampling (IS) scheme [14] that yields a discrete approximation of a target probability distribution. The PMC algorithm has important advantages with respect to MCMC techniques. It provides independent samples and asymptotically unbiased estimates at all iterations, which avoids the need of a convergence period. Additionally, PMC may be easily parallelized.On the other hand, the main weakness of IS and PMC is their low efficiency in high dimensional problems, due to the well-known degeneracy problem [15]. The recently proposed nonlinear PMC (NPMC) scheme [16, 17] mitigates this difficulty by computing nonlinear transformations of the importance weights (IWs), in order to smooth their variations and avoid degeneracy. In [17] a rigorous proof of convergence of the nonlinear importance sampling scheme as the number of Monte Carlo samples increases is presented. Similarly to the pMCMC method in [5], the NPMC algorithm resorts to an SMC approximation of the posterior distribution of the populations to compute the IWs.In this paper we apply the NPMC method to the estimation of both the parameters and the unobserved populations in a complex SKM modelling an autoregulatory feedback network with several species and noisy observations. We present numerical results to compare the performance of the state-of-the-art pMCMC and the proposed NPMC in two scenarios of different dimension and with two different observation models. We show that the NPMC method outperforms the pMCMC method for the same computational cost.The rest of the paper is organized as follows. In Section2 we present an introduction to the basics of SKMs and the usual solutions to this Bayesian inference problem. In Sections 3 and 4 we describe the pMCMC and NPMC methods, respectively, when applied to the approximation of posterior distributions in SKMs. In Section 5 we numerically compare the performance of pMCMC and NPMC schemes when applied to an autoregulatory feedback network, an SKM that can be interpreted as a generalization of the paradigmatic Lotka-Volterra model. Finally, Section 6 is devoted to the conclusions.


## 2. Bayesian Inference for Stochastic Kinetic Models

### 2.1. Stochastic Kinetic Models

A SKM is a multivariate continuous-time jump process modelling the interactions among different species or molecules in ecological, biological, and chemical systems [1].Consider an ecological interaction network that describes the time evolution of the population ofV species x
1
,
…
,
x
V related by means of K rate equations r
1
,
…
,
r
K
(1)
r
1
:
p
11
x
1
+
p
12
x
2
+
…
+
p
1
V
x
V
→
c
1
q
11
x
1
+
q
12
x
2
+
…
+
q
1
V
x
V
,
r
2
:
p
21
x
1
+
p
22
x
2
+
…
+
p
2
V
x
V
→
c
2
q
21
x
1
+
q
22
x
2
+
…
+
q
2
V
x
V
,
⋮
⋮
r
K
:
p
K
1
x
1
+
p
K
2
x
2
+
…
+
p
K
V
x
V
→
c
K
q
K
1
x
1
+
q
K
2
x
2
+
…
+
q
K
V
x
V
,where c
k
>
0, k
=
1
,
…
,
K are the random constant rate parameters and p
k
v and q
k
v, k
=
1
,
…
,
K, v
=
1
,
…
,
V, denote the coefficients of the population of the v-th species, before and after the k-th interaction, respectively. We will refer to p
k
v and q
k
v as reactant and product coefficients, respectively, by analogy with terminology commonly used in biochemical reactions. A matrix P of size K
×
V contains the reactant coefficients p
k
v and, similarly, Q contains the product coefficients q
k
v. The stoichiometry matrix of size V
×
K is defined as S
=
(
Q
-
P
)
⊤. The vector c
=
[
c
1
,
…
,
c
K
]
⊤ contains the rate parameters.Letx
v
(
t
), v
=
1
,
…
,
V, denote the nonnegative, integer population of species x
v at time t, and let x
(
t
)
=
[
x
1
(
t
)
,
…
,
x
V
(
t
)
]
⊤ denote the state of the system at this time instant. Let x
n
=
[
x
1
,
n
,
…
,
x
V
,
n
]
⊤ denote the state of the system at discrete time instants t
=
n
Δ, n
=
1
,
…
,
N, i.e., x
v
,
n
=
x
v
(
n
Δ
) where Δ denotes a time-discretization period. We denote by x the V
N
×
1 vector containing the population of all species at N consecutive discrete time instants, i.e., x
=
[
x
1
⊤
,
…
,
x
N
⊤
]
⊤.Thek-th rate equation describes an interaction that takes place stochastically according to its instantaneous rate or hazard function (2)
h
k
t
=
c
k
∏
v
=
1
V
x
v
t
p
k
v
,
k
=
1
,
…
,
K
, where the product of binomial coefficients represents the number of combinations in which the k-th interaction can occur, as a function of the population of the reactant species with population x
v, v
=
1
,
…
,
V. We additionally define the vector h
(
t
)
=
[
h
1
(
t
)
,
…
,
h
K
(
t
)
]
⊤. The waiting time to the next interaction is exponentially distributed with parameter h
0
(
t
)
=
∑
k
=
1
K
h
k
(
t
), and the probability of each type of interaction is given by h
k
(
t
)
/
h
0
(
t
)
,
k
=
1
,
…
,
K.See, e.g., [1] for a detailed description and discussion of the class of SKMs.


### 2.2. Bayesian Inference for SKMs

We consider the log-transformed rate parametersθ
=
[
θ
1
,
…
,
θ
K
]
⊤, where θ
k
=
log
⁡
(
c
k
), k
=
1
,
…
,
K, with prior pdf p
(
θ
). The prior pdf (for simplicity of notation, in this section we use p to denote the pdfs in the model. We write conditional pdfs as p
y
∣
x, and joint densities as p
θ
=
p
θ
1
,
…
,
θ
K
. This is an argument-wise notation, hence p
(
θ
1
) denotes the distribution of θ
1
, possibly different from p
(
θ
2
)
) of the initial population vector x
0 is denoted by p
(
x
0
). We assume that a linear combination of the populations of a subset of species is observed at discrete time instants corrupted by Gaussian noise, i.e.,(3)
y
n
=
M
x
n
+
w
n
,
n
=
1
,
…
,
N
,where M is the observation matrix with dimensions D
×
V and w
n
~
N
D
(
w
n
;
0
,
Λ
) is a multivariate Gaussian noise component with zero mean and covariance matrix Λ. We denote the complete observation vector with dimensionDN × 1 as y
=
[
y
1
⊤
,
…
,
y
N
⊤
]
⊤.The dynamical behavior of an arbitrary SKM may be described in terms of the following set of equations (that describe a state-space model with unknown parameters):(4)
θ
~
p
θ
p
a
r
a
m
e
t
e
r
s
p
r
i
o
r
,
x
0
~
p
x
0
p
o
p
u
l
a
t
i
o
n
s
p
r
i
o
r
,
x
n
~
p
x
n
∣
x
n
-
1
,
θ
t
r
a
n
s
i
t
i
o
n
e
q
u
a
t
i
o
n
,
y
n
~
p
y
n
∣
x
n
o
b
s
e
r
v
a
t
i
o
n
e
q
u
a
t
i
o
n
,where p
(
x
n
∣
x
n
-
1
,
θ
) and p
(
y
n
∣
x
n
) denote the transition pdf and the likelihood function, respectively. The Gillespie algorithm [18] (see Appendix A for a detailed description) allows to perform exact forward simulations of arbitrary SKMs, drawing samples from the transition densities p
(
x
n
∣
x
n
-
1
,
θ
), n
=
1
,
…
,
N, given a set of log-rate parameters θ and an initial population x
0. While the assumption Gaussian observational noise may not be adequate for some applications of SKMs, it has been employed by other authors in the past [5, 11] and, in particular, it has been used for computer experiments involving the prokaryotic autoregulation model [5] that we study numerically in Section 5.In this paper, we aim at obtaining a Monte Carlo approximation of the full joint posterior distribution of the log-rate parametersθ and the populations x, with density(5)
p
θ
,
x
∣
y
∝
p
y
∣
x
p
x
∣
x
0
,
θ
p
x
0
p
θ
,given the prior distributions p
(
θ
) and p
(
x
0
), the transition pdf p
(
x
∣
x
0
,
θ
)
=
∏
n
=
1
N
p
(
x
n
∣
x
n
-
1
,
θ
), and the likelihood function p
(
y
∣
x
)
=
∏
n
=
1
N
p
(
y
n
∣
x
n
) constructed from (3). We note that this state-space model can represent the SKM reaction equations of Section 2.1, yet it is abstract enough to make it applicable to other, relatively similar, models. For example, the assumption of mass action kinetics is not explicitly needed for this representation; hence other reaction rates could be incorporated into the state-space model (and the inference algorithms to be derived from it).We are also interested in computing approximations of the posterior marginals of the rate parametersp
(
θ
∣
y
)
=
∫
p
(
θ
,
x
∣
y
)
d
x and the species populations p
(
x
∣
y
)
=
∫
p
(
θ
,
x
∣
y
)
d
θ as well as their moments (e.g., the posterior mean), which are of the form (6)
E
p
θ
∣
y
f
θ
=
∫
f
θ
p
θ
∣
y
d
θ
,
E
p
x
∣
y
f
x
=
∫
f
x
p
x
∣
y
d
x
,
r
e
s
p
e
c
t
i
v
e
l
y
, where f is a real, integrable function.Bayesian inference based on exact stochastic simulations fromp
(
x
n
∣
x
n
-
1
,
θ
) generated via the Gillespie algorithm often becomes practically intractable even for models of modest complexity [7]. Thus, it is very common to resort to a continuous approximation of the underlying stochastic process, which is known as the diffusion approximation [1]. This approximation is known to be poor in low-concentration scenarios and thus should be avoided for models involving species with a very low population. Alternatively, in [3] the authors propose a solution based on a moment closure approximation of the stochastic process. Other techniques for the approximate modelling of the underlying jump process that have become relatively popular in the past few years include Gaussian processes [12] and the linear noise approximation [8–10].This inference problem has been traditionally addressed using MCMC methods, and IS based schemes have been avoided due to their inefficiency in high dimensional spaces [1]. In [2] various MCMC algorithms are evaluated in data-poor scenarios. In [5] a likelihood-free pMCMC scheme [6] is applied to this problem. This method is, to the best of our knowledge, the most powerful, yet computationally expensive, method provided so far for this kind of applications.In [16] a NPMC scheme is proposed for the approximation of the marginal posterior pdf p
(
θ
∣
y
). The performance of the NPMC method is tested in a simple SKM known as predator-prey model [20], providing excellent results with a low computational cost.In this paper we compare the performances of the pMCMC and the NPMC methods in the approximation of the full joint posterior pdfp
(
θ
,
x
∣
y
) in (5), which allows performing Bayesian inference for the rate parameters θ and the full sample path x, including unobserved components.


## 2.1. Stochastic Kinetic Models

A SKM is a multivariate continuous-time jump process modelling the interactions among different species or molecules in ecological, biological, and chemical systems [1].Consider an ecological interaction network that describes the time evolution of the population ofV species x
1
,
…
,
x
V related by means of K rate equations r
1
,
…
,
r
K
(1)
r
1
:
p
11
x
1
+
p
12
x
2
+
…
+
p
1
V
x
V
→
c
1
q
11
x
1
+
q
12
x
2
+
…
+
q
1
V
x
V
,
r
2
:
p
21
x
1
+
p
22
x
2
+
…
+
p
2
V
x
V
→
c
2
q
21
x
1
+
q
22
x
2
+
…
+
q
2
V
x
V
,
⋮
⋮
r
K
:
p
K
1
x
1
+
p
K
2
x
2
+
…
+
p
K
V
x
V
→
c
K
q
K
1
x
1
+
q
K
2
x
2
+
…
+
q
K
V
x
V
,where c
k
>
0, k
=
1
,
…
,
K are the random constant rate parameters and p
k
v and q
k
v, k
=
1
,
…
,
K, v
=
1
,
…
,
V, denote the coefficients of the population of the v-th species, before and after the k-th interaction, respectively. We will refer to p
k
v and q
k
v as reactant and product coefficients, respectively, by analogy with terminology commonly used in biochemical reactions. A matrix P of size K
×
V contains the reactant coefficients p
k
v and, similarly, Q contains the product coefficients q
k
v. The stoichiometry matrix of size V
×
K is defined as S
=
(
Q
-
P
)
⊤. The vector c
=
[
c
1
,
…
,
c
K
]
⊤ contains the rate parameters.Letx
v
(
t
), v
=
1
,
…
,
V, denote the nonnegative, integer population of species x
v at time t, and let x
(
t
)
=
[
x
1
(
t
)
,
…
,
x
V
(
t
)
]
⊤ denote the state of the system at this time instant. Let x
n
=
[
x
1
,
n
,
…
,
x
V
,
n
]
⊤ denote the state of the system at discrete time instants t
=
n
Δ, n
=
1
,
…
,
N, i.e., x
v
,
n
=
x
v
(
n
Δ
) where Δ denotes a time-discretization period. We denote by x the V
N
×
1 vector containing the population of all species at N consecutive discrete time instants, i.e., x
=
[
x
1
⊤
,
…
,
x
N
⊤
]
⊤.Thek-th rate equation describes an interaction that takes place stochastically according to its instantaneous rate or hazard function (2)
h
k
t
=
c
k
∏
v
=
1
V
x
v
t
p
k
v
,
k
=
1
,
…
,
K
, where the product of binomial coefficients represents the number of combinations in which the k-th interaction can occur, as a function of the population of the reactant species with population x
v, v
=
1
,
…
,
V. We additionally define the vector h
(
t
)
=
[
h
1
(
t
)
,
…
,
h
K
(
t
)
]
⊤. The waiting time to the next interaction is exponentially distributed with parameter h
0
(
t
)
=
∑
k
=
1
K
h
k
(
t
), and the probability of each type of interaction is given by h
k
(
t
)
/
h
0
(
t
)
,
k
=
1
,
…
,
K.See, e.g., [1] for a detailed description and discussion of the class of SKMs.


## 2.2. Bayesian Inference for SKMs

We consider the log-transformed rate parametersθ
=
[
θ
1
,
…
,
θ
K
]
⊤, where θ
k
=
log
⁡
(
c
k
), k
=
1
,
…
,
K, with prior pdf p
(
θ
). The prior pdf (for simplicity of notation, in this section we use p to denote the pdfs in the model. We write conditional pdfs as p
y
∣
x, and joint densities as p
θ
=
p
θ
1
,
…
,
θ
K
. This is an argument-wise notation, hence p
(
θ
1
) denotes the distribution of θ
1
, possibly different from p
(
θ
2
)
) of the initial population vector x
0 is denoted by p
(
x
0
). We assume that a linear combination of the populations of a subset of species is observed at discrete time instants corrupted by Gaussian noise, i.e.,(3)
y
n
=
M
x
n
+
w
n
,
n
=
1
,
…
,
N
,where M is the observation matrix with dimensions D
×
V and w
n
~
N
D
(
w
n
;
0
,
Λ
) is a multivariate Gaussian noise component with zero mean and covariance matrix Λ. We denote the complete observation vector with dimensionDN × 1 as y
=
[
y
1
⊤
,
…
,
y
N
⊤
]
⊤.The dynamical behavior of an arbitrary SKM may be described in terms of the following set of equations (that describe a state-space model with unknown parameters):(4)
θ
~
p
θ
p
a
r
a
m
e
t
e
r
s
p
r
i
o
r
,
x
0
~
p
x
0
p
o
p
u
l
a
t
i
o
n
s
p
r
i
o
r
,
x
n
~
p
x
n
∣
x
n
-
1
,
θ
t
r
a
n
s
i
t
i
o
n
e
q
u
a
t
i
o
n
,
y
n
~
p
y
n
∣
x
n
o
b
s
e
r
v
a
t
i
o
n
e
q
u
a
t
i
o
n
,where p
(
x
n
∣
x
n
-
1
,
θ
) and p
(
y
n
∣
x
n
) denote the transition pdf and the likelihood function, respectively. The Gillespie algorithm [18] (see Appendix A for a detailed description) allows to perform exact forward simulations of arbitrary SKMs, drawing samples from the transition densities p
(
x
n
∣
x
n
-
1
,
θ
), n
=
1
,
…
,
N, given a set of log-rate parameters θ and an initial population x
0. While the assumption Gaussian observational noise may not be adequate for some applications of SKMs, it has been employed by other authors in the past [5, 11] and, in particular, it has been used for computer experiments involving the prokaryotic autoregulation model [5] that we study numerically in Section 5.In this paper, we aim at obtaining a Monte Carlo approximation of the full joint posterior distribution of the log-rate parametersθ and the populations x, with density(5)
p
θ
,
x
∣
y
∝
p
y
∣
x
p
x
∣
x
0
,
θ
p
x
0
p
θ
,given the prior distributions p
(
θ
) and p
(
x
0
), the transition pdf p
(
x
∣
x
0
,
θ
)
=
∏
n
=
1
N
p
(
x
n
∣
x
n
-
1
,
θ
), and the likelihood function p
(
y
∣
x
)
=
∏
n
=
1
N
p
(
y
n
∣
x
n
) constructed from (3). We note that this state-space model can represent the SKM reaction equations of Section 2.1, yet it is abstract enough to make it applicable to other, relatively similar, models. For example, the assumption of mass action kinetics is not explicitly needed for this representation; hence other reaction rates could be incorporated into the state-space model (and the inference algorithms to be derived from it).We are also interested in computing approximations of the posterior marginals of the rate parametersp
(
θ
∣
y
)
=
∫
p
(
θ
,
x
∣
y
)
d
x and the species populations p
(
x
∣
y
)
=
∫
p
(
θ
,
x
∣
y
)
d
θ as well as their moments (e.g., the posterior mean), which are of the form (6)
E
p
θ
∣
y
f
θ
=
∫
f
θ
p
θ
∣
y
d
θ
,
E
p
x
∣
y
f
x
=
∫
f
x
p
x
∣
y
d
x
,
r
e
s
p
e
c
t
i
v
e
l
y
, where f is a real, integrable function.Bayesian inference based on exact stochastic simulations fromp
(
x
n
∣
x
n
-
1
,
θ
) generated via the Gillespie algorithm often becomes practically intractable even for models of modest complexity [7]. Thus, it is very common to resort to a continuous approximation of the underlying stochastic process, which is known as the diffusion approximation [1]. This approximation is known to be poor in low-concentration scenarios and thus should be avoided for models involving species with a very low population. Alternatively, in [3] the authors propose a solution based on a moment closure approximation of the stochastic process. Other techniques for the approximate modelling of the underlying jump process that have become relatively popular in the past few years include Gaussian processes [12] and the linear noise approximation [8–10].This inference problem has been traditionally addressed using MCMC methods, and IS based schemes have been avoided due to their inefficiency in high dimensional spaces [1]. In [2] various MCMC algorithms are evaluated in data-poor scenarios. In [5] a likelihood-free pMCMC scheme [6] is applied to this problem. This method is, to the best of our knowledge, the most powerful, yet computationally expensive, method provided so far for this kind of applications.In [16] a NPMC scheme is proposed for the approximation of the marginal posterior pdf p
(
θ
∣
y
). The performance of the NPMC method is tested in a simple SKM known as predator-prey model [20], providing excellent results with a low computational cost.In this paper we compare the performances of the pMCMC and the NPMC methods in the approximation of the full joint posterior pdfp
(
θ
,
x
∣
y
) in (5), which allows performing Bayesian inference for the rate parameters θ and the full sample path x, including unobserved components.


## 3. Particle MCMC for SKM’s

The particle marginal Metropolis-Hastings (PMMH) algorithm is a pMCMC method originally proposed in [6] for Monte Carlo sampling from the full posterior distribution p
(
θ
,
x
∣
y
). The PMMH scheme suggests a proposal mechanism of the form q
(
θ
⋆
∣
θ
)
p
^
J
(
x
⋆
∣
y
,
θ
⋆
). To be specific, a new candidate in the parameter space, θ
⋆, is drawn from an arbitrary proposal distribution q
(
θ
⋆
∣
θ
), while the new candidate in the variable space, x
⋆, is generated using an approximation of the posterior marginal p
(
x
⋆
∣
y
,
θ
⋆
) constructed by means of an SMC algorithm (i.e., a particle filter) with J particles and denoted p
^
J
(
x
⋆
∣
y
,
θ
⋆
). The probability of accepting the proposed pair (
θ
⋆
,
x
⋆
) is (7)
m
i
n
1
,
p
^
J
y
∣
θ
⋆
p
θ
⋆
p
^
J
y
∣
θ
p
θ
×
q
θ
∣
θ
⋆
q
θ
⋆
∣
θ
, where p
^
J
(
y
∣
θ
⋆
) is an unbiased approximation of the marginal likelihood of θ
⋆ (i.e., p
(
y
∣
θ
⋆
)), computed, again, by way of a particle filter with J particles. The PMMH algorithm is reproduced in Algorithm 1, and the SMC approximations of p
(
y
∣
θ
∗
) and p
(
x
∗
∣
y
,
θ
∗
) are described in Appendix B. Full details can be found in [6]. Note that the forward simulation of the stochastic process in the particle filter may be performed exactly with the Gillespie algorithm or using a diffusion approximation.Algorithm 1:Particle MCMC algorithm targeting p
(
θ
,
x
∣
y
) [6].
Initialization  (i
=
0):
1. Sample θ
(
0
)
~
p
(
θ
) and
2. run a SMC scheme targeting p
(
x
∣
y
,
θ
(
0
)
). Draw x
(
0
)
~
p
^
J
(
x
∣
y
,
θ
(
0
)
) from the SMC approximation and let p
^
J
(
y
∣
θ
(
0
)
)
denote the marginal likelihood estimate.
Iteration  (i
=
1
,
…
,
I):
1. Sample θ
⋆
~
q
(
θ
⋆
∣
θ
(
i
-
1
)
) and
2. run a SMC scheme targeting p
(
x
∣
y
,
θ
⋆
). Draw x
⋆
~
p
^
J
(
x
∣
y
,
θ
⋆
), let p
^
J
(
y
∣
θ
⋆
) denote the marginal likelihood estimate,
and
3. with probability
m
i
n
1
,
p
^
J
(
y
∣
θ
⋆
)
p
(
θ
⋆
)
p
^
J
(
y
∣
θ
(
i
-
1
)
)
p
(
θ
(
i
-
1
)
)
×
q
(
θ
(
i
-
1
)
∣
θ
⋆
)
q
(
θ
⋆
∣
θ
(
i
-
1
)
)
accept the move setting θ
(
i
)
=
θ
⋆, x
(
i
)
=
x
⋆ and p
^
J
(
y
∣
θ
(
i
)
)
=
p
^
J
(
y
∣
θ
⋆
). Otherwise store the current values θ
(
i
)
=
θ
(
i
-
1
),
x
(
i
)
=
x
(
i
-
1
) and p
^
J
(
y
∣
θ
(
i
)
)
=
p
^
J
(
y
∣
θ
(
i
-
1
)
).In [5] the proposal is selected as a Gaussian random walk q
(
θ
⋆
∣
θ
)
=
N
K
(
θ
⋆
;
θ
,
γ
2
I
), ignoring the correlation structure of the target distribution. The variance γ
2 has to be tuned and partly determines the performance of the algorithm (see [21] for some practical advice on the choice of this parameter). Low values of γ result in a poor exploration of the space of θ, while high γ values yield low acceptance rates. In both situations the resulting chain presents highly correlated samples and slow mixing properties. To reduce the correlation among samples it is common tothin the obtained chain, keeping a subset of equally spaced samples and discarding the rest. After removing the initial burn-in samples and thinning the output, we obtain a Markov chain {
θ
(
i
)
,
x
(
i
)
}
i
=
1
M with M correlated samples. Then, we may construct a sample approximation of the marginal posterior distributions of the parameters θ and the populations x, as (8)
p
^
M
d
θ
∣
y
=
1
M
∑
i
=
1
M
δ
θ
i
d
θ
,
p
^
M
d
x
∣
y
=
1
M
∑
i
=
1
M
δ
x
i
d
x
, respectively, where δ
θ
(
i
) and δ
x
(
i
) denote the unit delta measure centered at θ
(
i
) and x
(
i
), respectively. The approximation of the full joint posterior has the form (9)
p
^
M
d
θ
,
d
x
∣
y
=
1
M
∑
i
=
1
M
δ
θ
i
,
x
i
d
θ
,
d
x
.


## 4. Nonlinear PMC for SKM’s

The PMC method [13] is an iterative IS scheme that generates a sequence of proposal pdf’s q
l
(
·
), l
=
1
,
…
,
L, that approximate a target pdf π along the iterations. The NPMC scheme is proposed in [16] and it introduces nonlinearly transformed IWs (TIWs) in order to mitigate the numerical problems caused by degeneracy in the proposal update scheme.


### 4.1. NPMC Targetingp
(
θ
∣
y
)

We first consider as a target density the marginal posterior pdf of the parametersθ given the observation vector y, i.e., π
(
θ
)
=
p
(
θ
∣
y
). As in [16], we construct the proposal pdf q
l
(
θ
), l
=
2
,
…
,
L, as a Gaussian approximation of the target pdf obtained at the previous iteration l
-
1, whose mean and covariance parameters are selected to match the moments of the previous sample set. The NPMC algorithm is displayed in Algorithm 2. Details and some simple convergence results can be found in [16]. See [17] for a rigorous analysis.Algorithm 2:Nonlinear PMC targeting π
(
θ
)
=
p
(
θ
∣
y
).
Iteration (l
=
1
,
…
,
L):
1. Draw a set of M samples {
θ
l
(
i
)
}
i
=
1
M from the proposal density q
l
(
θ
):
(i) at iteration l
=
1, let q
1
(
θ
)
=
p
(
θ
).
(ii) at iterations l
=
2
,
…
,
L the proposal q
l
(
θ
) is the Gaussian approximation of p
(
θ
∣
y
) obtained at iteration l
-
1.
2. For i
=
1
,
…
,
M, run a SMC scheme with J particles targeting p
(
x
∣
y
,
θ
l
(
i
)
) and compute the marginal likelihood
estimate p
^
l
J
(
y
∣
θ
l
(
i
)
).
3. For i
=
1
,
…
,
M, compute the unnormalized IWs
w
l
(
i
)
∗
∝
p
^
l
J
(
y
∣
θ
l
(
i
)
)
p
(
θ
l
(
i
)
)
q
l
(
θ
l
(
i
)
)
.
4. For i
=
1
,
…
,
M, compute normalized TIWs, w
-
l
(
i
), by clipping the original IWs as
w
-
l
i
∗
=
min
⁡
w
l
i
∗
,
T
l
M
T
,
w
-
l
(
i
)
=
w
-
l
i
∗
∑
j
=
1
M
w
-
l
(
j
)
∗,
where the threshold value T
l
M
T denotes the M
T-th highest unnormalized IW w
l
(
i
)
∗, with 1
<
M
T
<
M.
5. Resample to obtain an unweighted set {
θ
~
l
(
i
)
}
i
=
1
M: for i
,
j
=
1
,
…
,
M, let θ
~
l
(
i
)
=
θ
l
(
j
) with probability w
-
l
(
j
).
6. Construct a Gaussian approximation q
l
+
1
(
θ
)
=
N
(
θ
;
μ
l
,
Σ
l
) of the posterior p
(
θ
∣
y
), where the mean vector and
covariance matrix are computed as
μ
l
=
1
M
∑
i
=
1
M
θ
~
l
i
a
n
d
Σ
l
=
1
M
∑
i
=
1
M
(
θ
~
l
(
i
)
-
μ
l
)
θ
~
l
i
-
μ
l
⊤
.
(
∗
)Similar to the pMCMC algorithm, in the NPMC implementation the densitiesp
(
x
∣
y
,
θ
) and p
(
y
∣
θ
) required in steps 2 and 3 are replaced by their SMC approximations, which are given in Appendix B. The NPMC method may also use either exact or approximate samples of the stochastic process, depending on the computational capabilities.For theclipping procedure performed in step 4 we consider, at each iteration l, a permutation i
1
,
…
,
i
M of the indices in {
1
,
…
,
M
} such that w
l
(
i
1
)
∗
≥
…
≥
w
l
(
i
M
)
∗ and choose aclipping parameter M
T
<
M. We select a threshold value T
l
M
=
w
l
(
i
M
T
)
∗ andclip the largest IWs, w
l
(
i
k
)
∗
≥
T
l
M, k
=
1
,
…
,
M
T
-
1. This transformation leads to M
T flat TIWs in the region of interest of θ, allowing for a robust update of the proposal. The performance of the algorithm shows little sensitivity to the selection of theclipping parameter M
T [16]. For simplicity, step 5 performs multinomial resampling.At each iteration of the NPMC algorithm we may construct a discrete approximation of the posterior probability distribution described by the pdfp
(
θ
∣
y
), based on the set of samples and TIWs, as (10)
p
^
l
M
d
θ
∣
y
=
∑
i
=
1
M
w
¯
l
i
δ
θ
l
i
d
θ
.The choice of a Gaussian approximation of the proposal q
l
+
1
(
θ
) in step 6 is arbitrary. We have adopted it here for simplicity. Other families of pdfs can be used without modifying the rest of the algorithm.


### 4.2. NPMC Targetingp
(
θ
,
x
∣
y
)

The NPMC method proposed in [16] may be readily applied to the approximation of the full joint posterior p
(
θ
,
x
∣
y
), in an manner equivalent to the pMCMC algorithm. We consider a sampling mechanism of the form q
(
θ
)
p
^
J
(
x
∣
y
,
θ
), where samples θ
(
i
) are again generated from the latest proposal q
(
θ
) and x
(
i
) are drawn from the SMC approximation p
^
J
(
x
∣
y
,
θ
(
i
)
) obtained via particle filtering (the iteration index has been omitted for simplicity). Then, the standard, unnormalized IW associated with the pair (
θ
(
i
)
,
x
(
i
)
) is computed as (11)
w
i
∗
=
p
^
J
θ
i
,
x
i
∣
y
q
θ
i
p
^
J
x
i
∣
y
,
θ
i
∝
p
^
J
x
i
,
y
∣
θ
i
p
θ
i
q
θ
i
p
^
J
x
i
∣
y
,
θ
i
∝
p
^
J
y
∣
θ
i
p
θ
i
q
θ
i and is independent of x. This reveals that, when samples x
l
(
i
) are drawn from p
^
J
(
d
x
∣
y
,
θ
), the algorithm yields a discrete approximation of the posterior distribution of the unobserved populations x constructed as (12)
p
^
l
M
d
x
∣
y
=
∑
i
=
1
M
w
¯
l
i
δ
x
l
i
d
x
.Even though the proposed NPMC and the pMCMC require very similar computations for each pair of samples of{
θ
,
x
} and thus have an equivalent computational cost, the NPMC has a set of important advantages with respect to its MCMC counterpart. PMC methods provide independent sets of samples at all iterations and do not require a burn-in period. On the other hand, the nonlinearity applied in the NPMC mitigates weight degeneracy, which is the main problem arising in conventional IS based methods, dramatically increasing its efficiency in high dimensional problems. As a consequence, we claim that the total number of samples, M
L (and thus, the computational complexity), required by the NPMC can be significantly lower than that of pMCMC, I. Additionally, contrary to pMCMC, which requires a careful choice of the proposal tuning parameter, the proposed method does not require the accurate fitting of any parameters.The NPMC method processes a set ofM i.i.d. samples at each iteration, requiring a low number of iterations (around 10 for the type of problems addressed here) for convergence to the target distribution. The bulk of the computational cost of NPMC, as well as of pMCMC, is the SMC approximation of the likelihood function. In the pMCMC algorithm the samples θ
(
i
) are processed sequentially (one after the other), and this process cannot be parallelized. On the contrary, at each iteration l of the NPMC method, the process of drawing M samples θ
l
(
i
) and computing the associated IWs w
l
(
i
)
∗ can be performed independently for each sample i. Thus, steps 1, 2, and 3 can be easily parallelized, reducing the total execution time up to that of a single sample θ
l
(
i
). On the other hand, steps 4, 5, and 6 require the complete set of samples and weights {
θ
l
(
i
)
,
w
l
(
i
)
∗
}
i
=
1
M and must be performed in a centralized manner. However, these computations have a negligible cost in comparison with a single likelihood approximation. Thus, the parallelization of the NPMC method can allow for a reduction in execution time up to a factor ≈
1
/
M. Note that we refer here to the parallelization of the PMC method and not of the SMC filter used to approximate the likelihood function.An extensive numerical comparison of pMCMC versus NPMC for a complex SKM that can be interpreted as a generalization of the paradigmatic Lotka-Volterra model is presented in Section5.


## 4.1. NPMC Targetingp
(
θ
∣
y
)

We first consider as a target density the marginal posterior pdf of the parametersθ given the observation vector y, i.e., π
(
θ
)
=
p
(
θ
∣
y
). As in [16], we construct the proposal pdf q
l
(
θ
), l
=
2
,
…
,
L, as a Gaussian approximation of the target pdf obtained at the previous iteration l
-
1, whose mean and covariance parameters are selected to match the moments of the previous sample set. The NPMC algorithm is displayed in Algorithm 2. Details and some simple convergence results can be found in [16]. See [17] for a rigorous analysis.Algorithm 2:Nonlinear PMC targeting π
(
θ
)
=
p
(
θ
∣
y
).
Iteration (l
=
1
,
…
,
L):
1. Draw a set of M samples {
θ
l
(
i
)
}
i
=
1
M from the proposal density q
l
(
θ
):
(i) at iteration l
=
1, let q
1
(
θ
)
=
p
(
θ
).
(ii) at iterations l
=
2
,
…
,
L the proposal q
l
(
θ
) is the Gaussian approximation of p
(
θ
∣
y
) obtained at iteration l
-
1.
2. For i
=
1
,
…
,
M, run a SMC scheme with J particles targeting p
(
x
∣
y
,
θ
l
(
i
)
) and compute the marginal likelihood
estimate p
^
l
J
(
y
∣
θ
l
(
i
)
).
3. For i
=
1
,
…
,
M, compute the unnormalized IWs
w
l
(
i
)
∗
∝
p
^
l
J
(
y
∣
θ
l
(
i
)
)
p
(
θ
l
(
i
)
)
q
l
(
θ
l
(
i
)
)
.
4. For i
=
1
,
…
,
M, compute normalized TIWs, w
-
l
(
i
), by clipping the original IWs as
w
-
l
i
∗
=
min
⁡
w
l
i
∗
,
T
l
M
T
,
w
-
l
(
i
)
=
w
-
l
i
∗
∑
j
=
1
M
w
-
l
(
j
)
∗,
where the threshold value T
l
M
T denotes the M
T-th highest unnormalized IW w
l
(
i
)
∗, with 1
<
M
T
<
M.
5. Resample to obtain an unweighted set {
θ
~
l
(
i
)
}
i
=
1
M: for i
,
j
=
1
,
…
,
M, let θ
~
l
(
i
)
=
θ
l
(
j
) with probability w
-
l
(
j
).
6. Construct a Gaussian approximation q
l
+
1
(
θ
)
=
N
(
θ
;
μ
l
,
Σ
l
) of the posterior p
(
θ
∣
y
), where the mean vector and
covariance matrix are computed as
μ
l
=
1
M
∑
i
=
1
M
θ
~
l
i
a
n
d
Σ
l
=
1
M
∑
i
=
1
M
(
θ
~
l
(
i
)
-
μ
l
)
θ
~
l
i
-
μ
l
⊤
.
(
∗
)Similar to the pMCMC algorithm, in the NPMC implementation the densitiesp
(
x
∣
y
,
θ
) and p
(
y
∣
θ
) required in steps 2 and 3 are replaced by their SMC approximations, which are given in Appendix B. The NPMC method may also use either exact or approximate samples of the stochastic process, depending on the computational capabilities.For theclipping procedure performed in step 4 we consider, at each iteration l, a permutation i
1
,
…
,
i
M of the indices in {
1
,
…
,
M
} such that w
l
(
i
1
)
∗
≥
…
≥
w
l
(
i
M
)
∗ and choose aclipping parameter M
T
<
M. We select a threshold value T
l
M
=
w
l
(
i
M
T
)
∗ andclip the largest IWs, w
l
(
i
k
)
∗
≥
T
l
M, k
=
1
,
…
,
M
T
-
1. This transformation leads to M
T flat TIWs in the region of interest of θ, allowing for a robust update of the proposal. The performance of the algorithm shows little sensitivity to the selection of theclipping parameter M
T [16]. For simplicity, step 5 performs multinomial resampling.At each iteration of the NPMC algorithm we may construct a discrete approximation of the posterior probability distribution described by the pdfp
(
θ
∣
y
), based on the set of samples and TIWs, as (10)
p
^
l
M
d
θ
∣
y
=
∑
i
=
1
M
w
¯
l
i
δ
θ
l
i
d
θ
.The choice of a Gaussian approximation of the proposal q
l
+
1
(
θ
) in step 6 is arbitrary. We have adopted it here for simplicity. Other families of pdfs can be used without modifying the rest of the algorithm.


## 4.2. NPMC Targetingp
(
θ
,
x
∣
y
)

The NPMC method proposed in [16] may be readily applied to the approximation of the full joint posterior p
(
θ
,
x
∣
y
), in an manner equivalent to the pMCMC algorithm. We consider a sampling mechanism of the form q
(
θ
)
p
^
J
(
x
∣
y
,
θ
), where samples θ
(
i
) are again generated from the latest proposal q
(
θ
) and x
(
i
) are drawn from the SMC approximation p
^
J
(
x
∣
y
,
θ
(
i
)
) obtained via particle filtering (the iteration index has been omitted for simplicity). Then, the standard, unnormalized IW associated with the pair (
θ
(
i
)
,
x
(
i
)
) is computed as (11)
w
i
∗
=
p
^
J
θ
i
,
x
i
∣
y
q
θ
i
p
^
J
x
i
∣
y
,
θ
i
∝
p
^
J
x
i
,
y
∣
θ
i
p
θ
i
q
θ
i
p
^
J
x
i
∣
y
,
θ
i
∝
p
^
J
y
∣
θ
i
p
θ
i
q
θ
i and is independent of x. This reveals that, when samples x
l
(
i
) are drawn from p
^
J
(
d
x
∣
y
,
θ
), the algorithm yields a discrete approximation of the posterior distribution of the unobserved populations x constructed as (12)
p
^
l
M
d
x
∣
y
=
∑
i
=
1
M
w
¯
l
i
δ
x
l
i
d
x
.Even though the proposed NPMC and the pMCMC require very similar computations for each pair of samples of{
θ
,
x
} and thus have an equivalent computational cost, the NPMC has a set of important advantages with respect to its MCMC counterpart. PMC methods provide independent sets of samples at all iterations and do not require a burn-in period. On the other hand, the nonlinearity applied in the NPMC mitigates weight degeneracy, which is the main problem arising in conventional IS based methods, dramatically increasing its efficiency in high dimensional problems. As a consequence, we claim that the total number of samples, M
L (and thus, the computational complexity), required by the NPMC can be significantly lower than that of pMCMC, I. Additionally, contrary to pMCMC, which requires a careful choice of the proposal tuning parameter, the proposed method does not require the accurate fitting of any parameters.The NPMC method processes a set ofM i.i.d. samples at each iteration, requiring a low number of iterations (around 10 for the type of problems addressed here) for convergence to the target distribution. The bulk of the computational cost of NPMC, as well as of pMCMC, is the SMC approximation of the likelihood function. In the pMCMC algorithm the samples θ
(
i
) are processed sequentially (one after the other), and this process cannot be parallelized. On the contrary, at each iteration l of the NPMC method, the process of drawing M samples θ
l
(
i
) and computing the associated IWs w
l
(
i
)
∗ can be performed independently for each sample i. Thus, steps 1, 2, and 3 can be easily parallelized, reducing the total execution time up to that of a single sample θ
l
(
i
). On the other hand, steps 4, 5, and 6 require the complete set of samples and weights {
θ
l
(
i
)
,
w
l
(
i
)
∗
}
i
=
1
M and must be performed in a centralized manner. However, these computations have a negligible cost in comparison with a single likelihood approximation. Thus, the parallelization of the NPMC method can allow for a reduction in execution time up to a factor ≈
1
/
M. Note that we refer here to the parallelization of the PMC method and not of the SMC filter used to approximate the likelihood function.An extensive numerical comparison of pMCMC versus NPMC for a complex SKM that can be interpreted as a generalization of the paradigmatic Lotka-Volterra model is presented in Section5.


## 5. Example: An Autoregulatory Feedback Network

In this section, we compare the performance of the pMCMC and the NPMC methods when applied to the problem of approximating the posterior distributions of the log-rate parametersp
(
θ
∣
y
) and the populations p
(
x
∣
y
) in a complex SKM given some observed data y. This model is significantly more involved than the standard predator-prey system. It has been introduced in [7] and further analyzed in [1, 5] for biochemical processes. In particular, it has also been considered as a model for representing the mechanisms for autoregulation in prokaryotes. The model is minimal in terms of the level of details included and offers a simplistic view of the mechanisms involved in gene autoregulation. However, it contains many of the interesting features of general autoregulatory feedback networks and it does provide sufficient detail to capture the network dynamics.


### 5.1. Prokaryotic Autoregulatory Model

The prokaryotic autoregulatory model is a SKM that involvesV
=
5 chemical species and K
=
8 reaction equations, r
1
,
…
,
r
K, given by [7] (13)
r
1
:
x
D
N
A
+
x
P
2
→
c
1
x
D
N
A
·
P
2
,
r
2
:
x
D
N
A
·
P
2
→
c
2
x
D
N
A
+
x
P
2
,
r
3
:
x
D
N
A
→
c
3
x
D
N
A
+
x
R
N
A
,
r
4
:
x
R
N
A
→
c
4
x
R
N
A
+
x
P
,
r
5
:
2
x
P
→
c
5
x
P
2
,
r
6
:
x
P
2
→
c
6
2
x
P
,
r
7
:
x
R
N
A
→
c
7
0
,
r
8
:
x
P
→
c
8
0
.We construct theV-dimensional vector containing the population of all species at time instant t as x
(
t
)
=
[
x
R
N
A
(
t
)
,
x
P
(
t
)
,
x
P
2
(
t
)
,
x
D
N
A
·
P
2
(
t
)
,
x
D
N
A
(
t
)
]
⊤. Thus, we obtain a stoichiometry matrix of the form (14)
S
=
0
0
1
0
0
0
-
1
0
0
0
0
1
-
2
2
0
-
1
-
1
1
0
0
1
-
1
0
0
1
-
1
0
0
0
0
0
0
-
1
1
0
0
0
0
0
0 and the hazard vector is given by(15)
h
t
=
c
1
x
D
N
A
x
P
2
,
c
2
x
D
N
A
·
P
2
,
c
3
x
D
N
A
,
c
4
x
R
N
A
,
c
5
x
P
x
P
-
1
2
,
c
6
x
P
2
,
c
7
x
R
N
A
,
c
8
x
P
⊤
,where the time dependance of the population of each species is omitted for notational simplicity.This model involves a conservation law given by the relationx
D
N
A
·
P
2
+
x
D
N
A
=
C, where C is the number of copies of this gene in the genome. We could use this relation to remove x
D
N
A
·
P
2 from the model, replacing any occurrences of the latter in the hazard function with C
-
x
D
N
A, but in this paper we abide by the notation in equation (15). Further details of this model can be found in [1].


### 5.2. Simulation Setup

We have selected most of the simulation parameters following [5]. The true vector of rate parameters which we aim to estimate has been set to (16)
c
=
0.1,0.7,0.35,0.2,0.1,0.9,0.3,0.1
⊤
, which yields log-transformed rate parameters (17)
θ
=
-
2.30,0.36,1.05,1.61,2.30,0.10,1.20,2.30
⊤
. The initial populations and the conservation constant have been set to x
0
=
[
x
1
(
0
)
,
…
,
x
V
(
0
)
]
⊤
=
[
8,8
,
8,5
,
5
]
⊤ and C
=
10, respectively. In all the simulations in this paper we have performed exact sampling from the stochastic model with the Gillespie algorithm to obtain the likelihood approximation via particle filtering. The number of particles for the SMC approximation p
^
J
(
x
∣
θ
,
y
) has been set to J
=
100 for all the simulations (both the pMCMC and the NPMC algorithms converge (as I
→
∞ and M
→
∞, respectively) even for fixed J[6, 17]. We have tested numerically that for J
>
100 performance hardly improves, while running times obviously become larger).Independent uniform priorsU
(
θ
k
;
-
7,2
) have been taken for each θ
k
=
log
⁡
(
c
k
). Opposite to [5], the initial populations x
0 are assumed unknown for the inference algorithm, which employs independent Poisson priors p
(
x
v
(
0
)
)
=
P
(
x
v
(
0
)
;
λ
v
), with the λ
v parameters set to the true initial populations, that is, λ
v
=
x
v
(
0
), v
=
1
,
…
,
V. Following [5], the conservation constant C is assumed to be known in the simulations.We consider two different observation scenarios. In the complete observation (CO) scenario we assume that all speciesx
v, v
=
1
,
…
,
V, are observed at regular time intervals of length Δ
=
1 and corrupted by Gaussian noise with variance σ
2
=
4 (assumed to be known). Thus, the observation matrix is of the form M
=
I
V and the observations are given by (18)
y
n
=
x
n
+
w
n
,
n
=
1
,
…
,
N
. In the CO case the complete vector of observations y
=
[
y
1
⊤
,
…
,
y
N
⊤
]
⊤ has dimension V
N
×
1.In the partial observation scenario (PO) only a linear combination of the proteinsx
P
+
2
x
P
2 is observed, also contaminated by Gaussian noise, i.e., the observation matrix is given by M
=
[
0,1
,
2,0
,
0
] (with dimension 1
×
V) and the observations are generated as (19)
y
n
=
x
2
,
n
+
2
x
3
,
n
+
w
n
,
w
h
e
r
e
w
n
~
N
1
w
n
;
0
,
σ
2
. In the PO case, a vector of scalar observations with dimension N
×
1 is constructed as y
=
[
y
1
,
…
,
y
N
]
⊤.We remark that the assumption of the observation noise terms (eitherw
n or w
n) displaying a Gaussian distribution is not needed to apply neither the pMCMC method nor the NPMC algorithm. Both techniques can be applied as long as the joint likelihood p
(
y
n
∣
x
n
,
θ
) of the state and the parameters can be evaluated, numerically and point-wise (either exactly or approximately; see [22]). In particular, it is possible to plug a Gaussian-process estimate of the likelihood [12] into the NPMC algorithm while running the rest of the algorithm in the same way as described in Section 4 (and still rely on the theoretical guarantees introduced in [22]).


### 5.3. Performance Evaluation

To evaluate the performance of the pMCMC and the NPMC methods we compute, in all the simulation runs, the mean square error (MSE) attained by the sample set that approximates the marginal posterior ofθ, generated by both schemes.For the pMCMC method, we compute the MSE of each parameterθ
k based on the M-size final output (after removing the burn-in period and thinning), as (20)
M
S
E
k
=
1
M
∑
i
=
1
M
θ
k
i
-
θ
k
2
,
k
∈
1
,
…
,
K
.For the NPMC method, we compute the MSE associated with each parameterθ
k, k
=
1
,
…
,
K, based on the unweighted sample set at the l-th iteration {
θ
~
l
i
}
i
=
1
M
,
l
=
1
,
…
,
L, as (21)
M
S
E
l
,
k
=
1
M
∑
i
=
1
M
θ
~
l
,
k
i
-
θ
k
2
=
μ
l
,
k
-
θ
k
2
+
σ
l
,
k
2
, where μ
l
,
k is the k-th component of the mean vector μ
l and the variance term σ
l
,
k
2 is the (
k
,
k
) component of matrix Σ
l.We have chosen the MSE as a figure of merit because the observations are generated synthetically and, therefore, we have access to ground truth values for the rate parameters that can be used to compute MSEs. A direct comparison in terms of the posterior distributions (e.g., using the Kullback-Leibler divergence or the total variation distance) is not possible because the posterior distribution ofθ conditional on the observations y is intractable.However, the MSE cannot be computed in real problems, where the true parametersθ
k are unknown. To monitor the stability and the efficiency of the two sampling schemes based on the generated sample alone, we resort to the so-called normalized effective sample size (NESS), which is often defined differently for MCMC and IS schemes [23].In the MCMC literature, the NESS gives the relative size of an i.i.d. (independent and identically distributed) sample with the same variance as the current sample and thus indicates the loss in efficiency due to the use of a Markov chain [23]. For pMCMC we compute the NESS from the final chain (after removing the burn-in period and thinning) as (22)
M
n
e
f
f
=
1
1
+
2
∑
j
=
1
∞
ρ
^
j
, where ρ
^
j
=
corr
⁡
(
θ
(
0
)
,
θ
(
j
)
) is the average autocorrelation function (ACF) at lag j. For the computation of the NESS, we truncate j when ρ
^
(
j
)
<
0.1.For IS methods, the NESS may be interpreted as the relative size of a sample generated from the target distribution with the same variance as the current sample. Even when high values of the NESS do not guarantee a low approximation error, the NESS is often used as an indicator of the numerical stability of the algorithm [24]. It cannot be evaluated exactly but we may compute an approximation of the NESS at each iteration of the NPMC scheme based on the set of TIWs as (23)
M
l
n
e
f
f
=
1
M
∑
i
=
1
M
w
¯
l
i
2
,
l
=
1
,
…
,
L
.In the next subsection we present, together with other numerical results, some comparisons involving the NESS for the PMCMC and PMC methods. The comparison should be taken with some caution because, although the interpretation of the NESS is the same for both techniques, the estimators have a different form and they are derived in a different way. Nevertheless, we believe that the NESS is a good numerical indicator of degeneracy for the two schemes, because it directly reflects poor mixing (high correlation) of the Markov chain and concentration of the weight in importance samplers, and we plot it with that purpose. See [25] for more details on the NESS and related statistics.


### 5.4. Simulation Results

We consider two simulation scenarios in which a different number of parameters is estimated.


#### 5.4.1. Estimation of a Single Rate Parameterθ
1

In this section we present numerical results regarding the approximation of the posterior distributionp
(
θ
1
,
x
∣
θ
∖
1
,
y
) of a single rate parameter θ
1
=
log
⁡
c
1 and the populations x, when the rest of parameters θ
∖
1
=
[
θ
2
,
…
,
θ
K
]
⊤, are assumed known.We compare the pMCMC and the NPMC methods in this simple scenario in order to show that the two algorithms perform almost equivalently in low dimensional problems. This is in a clear contrast with the results presented in Section5.4.2, which show that the NPMC method can be more efficient in higher dimensional settings.We have performedP
=
100 independent simulation runs of the pMCMC and the NPMC schemes in the CO and the PO scenarios, with independent population and observation vectors in each simulation p
=
1
,
…
,
P. Both in the CO and the PO cases, the same true population trajectories x
(
p
) were used, but the observations in the CO scenario, y
C
O
(
p
), and in the PO scenario, y
P
O
(
p
), differ. The number of observation times has been set to N
=
100 and exactly the same data has been used for NPMC and pMCMC.Following [5], as a proposal pdf q
(
θ
⋆
∣
θ
) in the pMCMC scheme we consider a Gaussian random walk update with variance γ
2. We have performed simulations with different γ values and the best results were obtained with γ
2
=
1, which are presented here. A total number of I
=
1
0
4 iterations have been run in each simulation. A final sample of size M
=
1
0
3 has been obtained from each Markov chain by discarding a burn-in period of 1
0
3 samples and thinning the output by a factor of 9.In the NPMC scheme, the number of iterations has been set toL
=
10, the number of samples per iteration is M
=
1
0
3 and theclipping parameter is M
T
=
100. In this way, the computational effort of the two methods is approximately the same, as they both generate I
=
M
L
=
1
0
4 samples in the space of θ.In Figure1 the final MSE obtained by the pMCMC (left) and the NPMC (right) algorithms for each simulation run is depicted versus the final NESS, in the CO and the PO scenarios. Note that the NESS is computed differently for pMCMC and NPMC. It can be observed that both algorithms perform similarly in this case, with an equivalent computational cost. Both algorithms attain on average lower MSE values in the CO scenario, as expected. However, the NESS also takes lower values in the CO case, which indicates a worse mixing of the Markov chains in the pMCMC algorithm and also higher degeneracy in the NPMC algorithm.Figure 1
Performance of the pMCMC (left) and the NPMC (right) methods for the estimation of a unique rate parameter θ
1: MSE (in logarithmic scale) obtained from the final output versus the NESS for each simulation run in the CO and the PO scenario. The big circles (CO) and squares (PO) represent simulation runs with a final mean MSE close to the global average.In Figure2 the evolution of the NESS (left) and the MSE (right) along the iterations of the NPMC algorithm is represented, for the CO and the PO scenarios. It can be observed that both measures attain a steady value by the 5-th iteration, both in the CO and the PO case, which suggests that actually less iterations are sufficient for this problem. Again, we observe that in the CO scenario both the NESS and the MSE reach lower values.Figure 2
Evolution along the iterations of the NPMC algorithm of the average NESS (left) and MSE (right) in the CO and PO scenarios, estimating a single parameter θ
1.Figure3 (left) plots the average ACF of the final pMCMC sample, after removing the burn-in period and thinning the Markov chain by a factor of 9. Particularly high correlations are present in the CO case, leading to a poor NESS. Related to the ACF, the average sample acceptance probability in the pMCMC scheme in the PO scenario is 0.091, while in the CO scenario it is only 0.0034, which means that 910 samples are accepted out of I
=
1
0
4 in the CO case and only 34 in the CO case.Figure 3
Left. Average ACF based on the final sample of size M
=
1
0
3 of the pMCMC scheme in the CO and the PO scenarios, averaged over P
=
100 simulation runs.Right. Marginal posterior pdf estimates p
^
(
θ
1
,
∣
θ
∖
1
,
y
) of an average simulation run, for pMCMC and NPMC in the CO and PO scenarios. The true value θ
1 is also shown.In Figure3 (right) the final pdf estimates p
^
(
θ
1
∣
θ
∖
1
,
y
) of the average simulation runs represented as big circles and crosses in Figure 1 are represented in the CO and the PO scenario, for the pMCMC and the NPMC schemes. For the pMCMC method we have built a Gaussian approximation of the posterior density p
(
θ
1
∣
θ
∖
1
,
y
) based on the final MCMC sample {
θ
1
(
i
)
}
i
=
1
M. For the NPMC method, this approximation corresponds to the proposal pdf for the next iteration L
+
1, i.e., p
^
(
θ
1
∣
θ
∖
1
,
y
)
=
q
L
+
1
(
θ
1
)
=
N
(
θ
1
;
μ
L
,
1
,
σ
L
,
1
2
), where the mean and variance terms μ
L
,
1 and σ
L
,
1
2 are computed as in ∗
. It can be observed in Figure 3 (right) that very similar results are obtained by both algorithms in this scenario. The final MSE values obtained by the pMCMC and the NPMC methods, averaged over P
=
100 simulation runs, are shown in Table 1, together with the MSE corresponding to the prior distribution.Table 1
Final mean and standard deviation (std) values of the MSE forθ
1 in the CO and PO scenarios, for pMCMC and NPMC. The prior values are included for comparison.
mean MSE
std MSE
Prior
6.789
0
PO
pMCMC
0.215
0.171
NPMC
0.195
0.170
CO
pMCMC
0.027
0.026
NPMC
0.022
0.016Figure4 depicts the posterior mean of the populations, x
^
=
E
p
^
(
x
∣
y
)
[
x
], obtained with pMCMC (left) as x
^
=
1
/
M
∑
i
=
1
M
x
(
i
) and with NPMC (right) as x
^
=
∑
i
=
1
M
w
¯
L
(
i
)
x
L
(
i
) in the PO scenario. The results correspond to the particular simulation runs (different for pMCMC and NPMC) identified with big squares in Figure 1 and whose posterior approximations, p
^
(
θ
1
∣
θ
∖
1
,
y
), are shown in Figure 3 (right). It can be observed that, in the PO scenario, the trend of the population of all the species is reasonably identified, even though only a linear combination of the proteins is observed. In the CO scenario the populations of all species are accurately estimated and are not shown for conciseness. Note that the populations of all species are very low, which suggests that the diffusion approximation may perform poorly in this scenario.Figure 4
Posterior mean,x
^
=
E
p
^
(
x
∣
y
)
[
x
], of the populations obtained in a single simulation run of pMCMC (left) and NPMC (right) in the PO scenario (only a linear combination of the proteins is observed, corrupted by noise).The results presented in this section reveal a very similar performance of the two methods in this simple scenario, which suggests that the parameters of pMCMC have been properly tuned. Also in terms of computational complexity pMCMC and NPMC perform very similarly. The execution time per1
0
3 samples (one NPMC iteration and 1
0
3 pMCMC iterations) for the pMCMC scheme is 312 seconds, while for NPMC it is 325 seconds, both in the CO and in the PO cases, on a 3-GHz Intel Core 2 Duo CPU, with 4 GB of RAM. The stochastic forward simulation of the prokaryotic model with the Gillespie algorithm has been implemented in C, and the rest of the code in Matlab R2014b. Out of the execution time per iteration of the NPMC method, only 0.06 second correspond to batch steps 4, 5, and 6, and the rest of operations can be easily parallelized, as previously discussed.However, the pMCMC method does not allow for an easy parallelization, provides a set of highly correlated samples (especially in the CO scenario) and requires the setting of the proposal varianceγ
2, the burn-in and the thinning parameters, which may not be straightforward, and determines the performance of the algorithm. On the contrary, the NPMC scheme provides uncorrelated sets of samples at each iteration and does not require the precise fitting of any parameters. Additionally, the computer simulations suggest that the convergence of the NPMC algorithm may be assessed observing the evolution of the NESS, which usually reaches a steady value simultaneously with the MSE.


#### 5.4.2. Estimation of All the Parametersθ
k
,
k
=
1
,
…
,
K

In this section we present simulation results to evaluate the performance of the pMCMC and the NPMC schemes in the approximation of the posterior distribution of the rate parameters and the populations of all species,p
(
θ
,
x
∣
y
), assuming that all the rate parameters are unknown, again in the CO and the PO scenarios.In this case,N
=
200 observation times are assumed for all the simulations. Again, P
=
100 independent simulation runs of each algorithm have been performed. The NPMC scheme has been run for L
=
15 iterations, with M
=
1
0
3 samples per iteration andclipping parameter M
T
=
100. The pMCMC scheme has been run with I
=
15
×
1
0
3 iterations in each simulation run, a burn-in period of 1
0
3 iterations, and thinning the output by a factor of 14. With this setup the computational effort is approximately the same in the two schemes. The variance of the random walk proposal has again been set to γ
2
=
1, since this value seems to yield the best performance.In Figure5 the MSE (in logarithmic scale), averaged over the parameters θ
k, attained by the pMCMC (left) and the NPMC (right) algorithms, is represented versus the NESS, in the CO and PO scenarios. Simulation runs which attained a final MSE close to the global average value are indicated with big circles (CO) and squares (PO) on both plots. It can be observed that the pMCMC method performs similarly in both scenarios, in terms of MSE and NESS, yielding poor results in both cases. On the contrary, the NPMC method provides significantly better MSE results in the CO scenario, where a larger amount of information is available. The NPMC method does not present degradation due to the high degeneracy occurring in the CO scenario.Figure 5
Performance of the pMCMC (left) and the NPMC (right) methods for the estimation of the whole set of rate parameters θ: MSE (in logarithmic scale) versus the final NESS, for each simulation run in the CO and the PO scenario. The big circles (CO) and squares (PO) represent simulation runs with a final mean MSE close to the global average.Figure6 depicts the evolution along the iterations of the NESS (left) and the MSE (right) averaged over P
=
100 independent simulation runs for the NPMC algorithm. Both indices converge to a steady value in a low number of iterations also in this complex scenario. As expected, a significantly higher final MSE is attained in the extremely data-poor PO scenario.Figure 6
Evolution along the NPMC iterations of the average NESS (left) and MSE (right) in the CO and the PO scenario.In Figure7 (left) the average ACF attained by the pMCMC in the CO and the PO cases is represented. Even after thinning the output, the sample correlation is extremely high in both scenarios, which leads to a very low NESS. The acceptance rate is also very low and very long chains are required to obtain reasonable results. In the PO scenario ≈
44 samples are accepted on average in a simulation run of I
=
15
×
1
0
3 samples (acceptance rate 0.0029). In the CO case, only ≈
23 samples are accepted on average (rate 0.0015).Figure 7
Left. Autocorrelations based on the final sample of size 1
0
3 of the pMCMC scheme in the CO and the PO scenarios, averaged over P
=
100 simulation runs.Right. Markov chain provided by the pMCMC method in the PO scenario, corresponding to the average simulation run depicted with a big square in Figure 5 (left).Figure7 (right) depicts the final Markov chain provided by the pMCMC method (after removing the burn-in period and thinning the output) in the average simulation run represented with a big square in Figure 5 (left). It can be observed that the mixing of the chain is very poor, with a total number of accepted samples of 46 (close to the average). Many other simulations, in both the PO and the CO scenarios, provide even lower number of accepted samples and thus very inconsistent results.Figure8 depicts the final Gaussian approximations of the marginal posteriors p
(
θ
k
∣
y
), k
=
1
,
…
,
8, obtained by the pMCMC and the NPMC methods, in the CO and PO scenarios, for the average simulation runs represented as big circles and squares in Figure 5. We can observe that the NPMC method provides a significantly better approximation of the log-rate parameters in the CO scenario, where a larger amount of data is available, which is also clear from Figure 5 (right). However, the pMCMC on average performs similarly in both scenarios, due to the low efficiency of the pMCMC sampling scheme when the dimension of the problem (either K or N) increases.Figure 8
Marginal posterior pdf approximations of each parameterp
^
(
θ
k
∣
y
), k
=
1
,
…
,
K, attained in an average simulation run by the pMCMC and the NPMC, in the CO and in the PO case.In Table2 the MSE of each parameter θ
k averaged over P
=
100 independent simulation runs is shown, as obtained with the pMCMC and the NPMC schemes, for the CO and the PO experiments. In the CO case, NPMC provides homogeneous results for all parameters. On the contrary, in the PO case, some of the parameters (specially θ
5 and θ
6) are significantly poorly estimated, presenting a final MSE close to the initial value (which corresponds to the prior knowledge). The pMCMC scheme presents significantly higher MSE values than NPMC in both observation scenarios and for all parameters θ
k.Table 2
Final MSE for the parametersθ
k
,
k
=
1
,
…
,
K in the CO and PO experiments, averaged over the simulation runs. The last two columns correspond to the mean and standard deviation (std) values of the global MSE (averaged over the parameters). The prior values are included for comparison.
θ
1
θ
2
θ
3
θ
4
θ
5
θ
6
θ
7
θ
8
mean MSE
std MSE
Prior
6.789
11.344
8.853
7.543
6.789
12.484
8.430
6.789
8.628
0
PO
pMCMC
3.412
3.319
5.543
3.200
7.059
8.929
6.799
4.371
5.329
2.926
NPMC
1.246
1.011
2.214
1.490
4.073
7.015
2.311
1.856
2.652
1.020
CO
pMCMC
2.899
2.958
1.676
1.572
1.604
1.547
1.573
1.468
1.912
1.476
NPMC
0.305
0.302
0.162
0.167
0.280
0.280
0.156
0.168
0.228
0.091Figure9 depicts the population posterior mean x
^
=
E
p
(
x
∣
y
)
[
x
] corresponding to the average simulation runs of the pMCMC and the NPMC methods in the PO scenario, represented as big squares in Figure 5. Again, the NPMC method provides more accurate estimates of the unobserved populations than the pMCMC method, especially for x
R
N
A. In the CO scenario both methods provide good approximations of the populations of all species.Figure 9
Posterior meanx
^
=
E
p
(
x
∣
y
)
[
x
] of the populations of all species obtained in the average simulation run of the pMCMC (left) and the NPMC (right) schemes, in the PO scenario.


## 5.1. Prokaryotic Autoregulatory Model

The prokaryotic autoregulatory model is a SKM that involvesV
=
5 chemical species and K
=
8 reaction equations, r
1
,
…
,
r
K, given by [7] (13)
r
1
:
x
D
N
A
+
x
P
2
→
c
1
x
D
N
A
·
P
2
,
r
2
:
x
D
N
A
·
P
2
→
c
2
x
D
N
A
+
x
P
2
,
r
3
:
x
D
N
A
→
c
3
x
D
N
A
+
x
R
N
A
,
r
4
:
x
R
N
A
→
c
4
x
R
N
A
+
x
P
,
r
5
:
2
x
P
→
c
5
x
P
2
,
r
6
:
x
P
2
→
c
6
2
x
P
,
r
7
:
x
R
N
A
→
c
7
0
,
r
8
:
x
P
→
c
8
0
.We construct theV-dimensional vector containing the population of all species at time instant t as x
(
t
)
=
[
x
R
N
A
(
t
)
,
x
P
(
t
)
,
x
P
2
(
t
)
,
x
D
N
A
·
P
2
(
t
)
,
x
D
N
A
(
t
)
]
⊤. Thus, we obtain a stoichiometry matrix of the form (14)
S
=
0
0
1
0
0
0
-
1
0
0
0
0
1
-
2
2
0
-
1
-
1
1
0
0
1
-
1
0
0
1
-
1
0
0
0
0
0
0
-
1
1
0
0
0
0
0
0 and the hazard vector is given by(15)
h
t
=
c
1
x
D
N
A
x
P
2
,
c
2
x
D
N
A
·
P
2
,
c
3
x
D
N
A
,
c
4
x
R
N
A
,
c
5
x
P
x
P
-
1
2
,
c
6
x
P
2
,
c
7
x
R
N
A
,
c
8
x
P
⊤
,where the time dependance of the population of each species is omitted for notational simplicity.This model involves a conservation law given by the relationx
D
N
A
·
P
2
+
x
D
N
A
=
C, where C is the number of copies of this gene in the genome. We could use this relation to remove x
D
N
A
·
P
2 from the model, replacing any occurrences of the latter in the hazard function with C
-
x
D
N
A, but in this paper we abide by the notation in equation (15). Further details of this model can be found in [1].


## 5.2. Simulation Setup

We have selected most of the simulation parameters following [5]. The true vector of rate parameters which we aim to estimate has been set to (16)
c
=
0.1,0.7,0.35,0.2,0.1,0.9,0.3,0.1
⊤
, which yields log-transformed rate parameters (17)
θ
=
-
2.30,0.36,1.05,1.61,2.30,0.10,1.20,2.30
⊤
. The initial populations and the conservation constant have been set to x
0
=
[
x
1
(
0
)
,
…
,
x
V
(
0
)
]
⊤
=
[
8,8
,
8,5
,
5
]
⊤ and C
=
10, respectively. In all the simulations in this paper we have performed exact sampling from the stochastic model with the Gillespie algorithm to obtain the likelihood approximation via particle filtering. The number of particles for the SMC approximation p
^
J
(
x
∣
θ
,
y
) has been set to J
=
100 for all the simulations (both the pMCMC and the NPMC algorithms converge (as I
→
∞ and M
→
∞, respectively) even for fixed J[6, 17]. We have tested numerically that for J
>
100 performance hardly improves, while running times obviously become larger).Independent uniform priorsU
(
θ
k
;
-
7,2
) have been taken for each θ
k
=
log
⁡
(
c
k
). Opposite to [5], the initial populations x
0 are assumed unknown for the inference algorithm, which employs independent Poisson priors p
(
x
v
(
0
)
)
=
P
(
x
v
(
0
)
;
λ
v
), with the λ
v parameters set to the true initial populations, that is, λ
v
=
x
v
(
0
), v
=
1
,
…
,
V. Following [5], the conservation constant C is assumed to be known in the simulations.We consider two different observation scenarios. In the complete observation (CO) scenario we assume that all speciesx
v, v
=
1
,
…
,
V, are observed at regular time intervals of length Δ
=
1 and corrupted by Gaussian noise with variance σ
2
=
4 (assumed to be known). Thus, the observation matrix is of the form M
=
I
V and the observations are given by (18)
y
n
=
x
n
+
w
n
,
n
=
1
,
…
,
N
. In the CO case the complete vector of observations y
=
[
y
1
⊤
,
…
,
y
N
⊤
]
⊤ has dimension V
N
×
1.In the partial observation scenario (PO) only a linear combination of the proteinsx
P
+
2
x
P
2 is observed, also contaminated by Gaussian noise, i.e., the observation matrix is given by M
=
[
0,1
,
2,0
,
0
] (with dimension 1
×
V) and the observations are generated as (19)
y
n
=
x
2
,
n
+
2
x
3
,
n
+
w
n
,
w
h
e
r
e
w
n
~
N
1
w
n
;
0
,
σ
2
. In the PO case, a vector of scalar observations with dimension N
×
1 is constructed as y
=
[
y
1
,
…
,
y
N
]
⊤.We remark that the assumption of the observation noise terms (eitherw
n or w
n) displaying a Gaussian distribution is not needed to apply neither the pMCMC method nor the NPMC algorithm. Both techniques can be applied as long as the joint likelihood p
(
y
n
∣
x
n
,
θ
) of the state and the parameters can be evaluated, numerically and point-wise (either exactly or approximately; see [22]). In particular, it is possible to plug a Gaussian-process estimate of the likelihood [12] into the NPMC algorithm while running the rest of the algorithm in the same way as described in Section 4 (and still rely on the theoretical guarantees introduced in [22]).


## 5.3. Performance Evaluation

To evaluate the performance of the pMCMC and the NPMC methods we compute, in all the simulation runs, the mean square error (MSE) attained by the sample set that approximates the marginal posterior ofθ, generated by both schemes.For the pMCMC method, we compute the MSE of each parameterθ
k based on the M-size final output (after removing the burn-in period and thinning), as (20)
M
S
E
k
=
1
M
∑
i
=
1
M
θ
k
i
-
θ
k
2
,
k
∈
1
,
…
,
K
.For the NPMC method, we compute the MSE associated with each parameterθ
k, k
=
1
,
…
,
K, based on the unweighted sample set at the l-th iteration {
θ
~
l
i
}
i
=
1
M
,
l
=
1
,
…
,
L, as (21)
M
S
E
l
,
k
=
1
M
∑
i
=
1
M
θ
~
l
,
k
i
-
θ
k
2
=
μ
l
,
k
-
θ
k
2
+
σ
l
,
k
2
, where μ
l
,
k is the k-th component of the mean vector μ
l and the variance term σ
l
,
k
2 is the (
k
,
k
) component of matrix Σ
l.We have chosen the MSE as a figure of merit because the observations are generated synthetically and, therefore, we have access to ground truth values for the rate parameters that can be used to compute MSEs. A direct comparison in terms of the posterior distributions (e.g., using the Kullback-Leibler divergence or the total variation distance) is not possible because the posterior distribution ofθ conditional on the observations y is intractable.However, the MSE cannot be computed in real problems, where the true parametersθ
k are unknown. To monitor the stability and the efficiency of the two sampling schemes based on the generated sample alone, we resort to the so-called normalized effective sample size (NESS), which is often defined differently for MCMC and IS schemes [23].In the MCMC literature, the NESS gives the relative size of an i.i.d. (independent and identically distributed) sample with the same variance as the current sample and thus indicates the loss in efficiency due to the use of a Markov chain [23]. For pMCMC we compute the NESS from the final chain (after removing the burn-in period and thinning) as (22)
M
n
e
f
f
=
1
1
+
2
∑
j
=
1
∞
ρ
^
j
, where ρ
^
j
=
corr
⁡
(
θ
(
0
)
,
θ
(
j
)
) is the average autocorrelation function (ACF) at lag j. For the computation of the NESS, we truncate j when ρ
^
(
j
)
<
0.1.For IS methods, the NESS may be interpreted as the relative size of a sample generated from the target distribution with the same variance as the current sample. Even when high values of the NESS do not guarantee a low approximation error, the NESS is often used as an indicator of the numerical stability of the algorithm [24]. It cannot be evaluated exactly but we may compute an approximation of the NESS at each iteration of the NPMC scheme based on the set of TIWs as (23)
M
l
n
e
f
f
=
1
M
∑
i
=
1
M
w
¯
l
i
2
,
l
=
1
,
…
,
L
.In the next subsection we present, together with other numerical results, some comparisons involving the NESS for the PMCMC and PMC methods. The comparison should be taken with some caution because, although the interpretation of the NESS is the same for both techniques, the estimators have a different form and they are derived in a different way. Nevertheless, we believe that the NESS is a good numerical indicator of degeneracy for the two schemes, because it directly reflects poor mixing (high correlation) of the Markov chain and concentration of the weight in importance samplers, and we plot it with that purpose. See [25] for more details on the NESS and related statistics.


## 5.4. Simulation Results

We consider two simulation scenarios in which a different number of parameters is estimated.


### 5.4.1. Estimation of a Single Rate Parameterθ
1

In this section we present numerical results regarding the approximation of the posterior distributionp
(
θ
1
,
x
∣
θ
∖
1
,
y
) of a single rate parameter θ
1
=
log
⁡
c
1 and the populations x, when the rest of parameters θ
∖
1
=
[
θ
2
,
…
,
θ
K
]
⊤, are assumed known.We compare the pMCMC and the NPMC methods in this simple scenario in order to show that the two algorithms perform almost equivalently in low dimensional problems. This is in a clear contrast with the results presented in Section5.4.2, which show that the NPMC method can be more efficient in higher dimensional settings.We have performedP
=
100 independent simulation runs of the pMCMC and the NPMC schemes in the CO and the PO scenarios, with independent population and observation vectors in each simulation p
=
1
,
…
,
P. Both in the CO and the PO cases, the same true population trajectories x
(
p
) were used, but the observations in the CO scenario, y
C
O
(
p
), and in the PO scenario, y
P
O
(
p
), differ. The number of observation times has been set to N
=
100 and exactly the same data has been used for NPMC and pMCMC.Following [5], as a proposal pdf q
(
θ
⋆
∣
θ
) in the pMCMC scheme we consider a Gaussian random walk update with variance γ
2. We have performed simulations with different γ values and the best results were obtained with γ
2
=
1, which are presented here. A total number of I
=
1
0
4 iterations have been run in each simulation. A final sample of size M
=
1
0
3 has been obtained from each Markov chain by discarding a burn-in period of 1
0
3 samples and thinning the output by a factor of 9.In the NPMC scheme, the number of iterations has been set toL
=
10, the number of samples per iteration is M
=
1
0
3 and theclipping parameter is M
T
=
100. In this way, the computational effort of the two methods is approximately the same, as they both generate I
=
M
L
=
1
0
4 samples in the space of θ.In Figure1 the final MSE obtained by the pMCMC (left) and the NPMC (right) algorithms for each simulation run is depicted versus the final NESS, in the CO and the PO scenarios. Note that the NESS is computed differently for pMCMC and NPMC. It can be observed that both algorithms perform similarly in this case, with an equivalent computational cost. Both algorithms attain on average lower MSE values in the CO scenario, as expected. However, the NESS also takes lower values in the CO case, which indicates a worse mixing of the Markov chains in the pMCMC algorithm and also higher degeneracy in the NPMC algorithm.Figure 1
Performance of the pMCMC (left) and the NPMC (right) methods for the estimation of a unique rate parameter θ
1: MSE (in logarithmic scale) obtained from the final output versus the NESS for each simulation run in the CO and the PO scenario. The big circles (CO) and squares (PO) represent simulation runs with a final mean MSE close to the global average.In Figure2 the evolution of the NESS (left) and the MSE (right) along the iterations of the NPMC algorithm is represented, for the CO and the PO scenarios. It can be observed that both measures attain a steady value by the 5-th iteration, both in the CO and the PO case, which suggests that actually less iterations are sufficient for this problem. Again, we observe that in the CO scenario both the NESS and the MSE reach lower values.Figure 2
Evolution along the iterations of the NPMC algorithm of the average NESS (left) and MSE (right) in the CO and PO scenarios, estimating a single parameter θ
1.Figure3 (left) plots the average ACF of the final pMCMC sample, after removing the burn-in period and thinning the Markov chain by a factor of 9. Particularly high correlations are present in the CO case, leading to a poor NESS. Related to the ACF, the average sample acceptance probability in the pMCMC scheme in the PO scenario is 0.091, while in the CO scenario it is only 0.0034, which means that 910 samples are accepted out of I
=
1
0
4 in the CO case and only 34 in the CO case.Figure 3
Left. Average ACF based on the final sample of size M
=
1
0
3 of the pMCMC scheme in the CO and the PO scenarios, averaged over P
=
100 simulation runs.Right. Marginal posterior pdf estimates p
^
(
θ
1
,
∣
θ
∖
1
,
y
) of an average simulation run, for pMCMC and NPMC in the CO and PO scenarios. The true value θ
1 is also shown.In Figure3 (right) the final pdf estimates p
^
(
θ
1
∣
θ
∖
1
,
y
) of the average simulation runs represented as big circles and crosses in Figure 1 are represented in the CO and the PO scenario, for the pMCMC and the NPMC schemes. For the pMCMC method we have built a Gaussian approximation of the posterior density p
(
θ
1
∣
θ
∖
1
,
y
) based on the final MCMC sample {
θ
1
(
i
)
}
i
=
1
M. For the NPMC method, this approximation corresponds to the proposal pdf for the next iteration L
+
1, i.e., p
^
(
θ
1
∣
θ
∖
1
,
y
)
=
q
L
+
1
(
θ
1
)
=
N
(
θ
1
;
μ
L
,
1
,
σ
L
,
1
2
), where the mean and variance terms μ
L
,
1 and σ
L
,
1
2 are computed as in ∗
. It can be observed in Figure 3 (right) that very similar results are obtained by both algorithms in this scenario. The final MSE values obtained by the pMCMC and the NPMC methods, averaged over P
=
100 simulation runs, are shown in Table 1, together with the MSE corresponding to the prior distribution.Table 1
Final mean and standard deviation (std) values of the MSE forθ
1 in the CO and PO scenarios, for pMCMC and NPMC. The prior values are included for comparison.
mean MSE
std MSE
Prior
6.789
0
PO
pMCMC
0.215
0.171
NPMC
0.195
0.170
CO
pMCMC
0.027
0.026
NPMC
0.022
0.016Figure4 depicts the posterior mean of the populations, x
^
=
E
p
^
(
x
∣
y
)
[
x
], obtained with pMCMC (left) as x
^
=
1
/
M
∑
i
=
1
M
x
(
i
) and with NPMC (right) as x
^
=
∑
i
=
1
M
w
¯
L
(
i
)
x
L
(
i
) in the PO scenario. The results correspond to the particular simulation runs (different for pMCMC and NPMC) identified with big squares in Figure 1 and whose posterior approximations, p
^
(
θ
1
∣
θ
∖
1
,
y
), are shown in Figure 3 (right). It can be observed that, in the PO scenario, the trend of the population of all the species is reasonably identified, even though only a linear combination of the proteins is observed. In the CO scenario the populations of all species are accurately estimated and are not shown for conciseness. Note that the populations of all species are very low, which suggests that the diffusion approximation may perform poorly in this scenario.Figure 4
Posterior mean,x
^
=
E
p
^
(
x
∣
y
)
[
x
], of the populations obtained in a single simulation run of pMCMC (left) and NPMC (right) in the PO scenario (only a linear combination of the proteins is observed, corrupted by noise).The results presented in this section reveal a very similar performance of the two methods in this simple scenario, which suggests that the parameters of pMCMC have been properly tuned. Also in terms of computational complexity pMCMC and NPMC perform very similarly. The execution time per1
0
3 samples (one NPMC iteration and 1
0
3 pMCMC iterations) for the pMCMC scheme is 312 seconds, while for NPMC it is 325 seconds, both in the CO and in the PO cases, on a 3-GHz Intel Core 2 Duo CPU, with 4 GB of RAM. The stochastic forward simulation of the prokaryotic model with the Gillespie algorithm has been implemented in C, and the rest of the code in Matlab R2014b. Out of the execution time per iteration of the NPMC method, only 0.06 second correspond to batch steps 4, 5, and 6, and the rest of operations can be easily parallelized, as previously discussed.However, the pMCMC method does not allow for an easy parallelization, provides a set of highly correlated samples (especially in the CO scenario) and requires the setting of the proposal varianceγ
2, the burn-in and the thinning parameters, which may not be straightforward, and determines the performance of the algorithm. On the contrary, the NPMC scheme provides uncorrelated sets of samples at each iteration and does not require the precise fitting of any parameters. Additionally, the computer simulations suggest that the convergence of the NPMC algorithm may be assessed observing the evolution of the NESS, which usually reaches a steady value simultaneously with the MSE.


### 5.4.2. Estimation of All the Parametersθ
k
,
k
=
1
,
…
,
K

In this section we present simulation results to evaluate the performance of the pMCMC and the NPMC schemes in the approximation of the posterior distribution of the rate parameters and the populations of all species,p
(
θ
,
x
∣
y
), assuming that all the rate parameters are unknown, again in the CO and the PO scenarios.In this case,N
=
200 observation times are assumed for all the simulations. Again, P
=
100 independent simulation runs of each algorithm have been performed. The NPMC scheme has been run for L
=
15 iterations, with M
=
1
0
3 samples per iteration andclipping parameter M
T
=
100. The pMCMC scheme has been run with I
=
15
×
1
0
3 iterations in each simulation run, a burn-in period of 1
0
3 iterations, and thinning the output by a factor of 14. With this setup the computational effort is approximately the same in the two schemes. The variance of the random walk proposal has again been set to γ
2
=
1, since this value seems to yield the best performance.In Figure5 the MSE (in logarithmic scale), averaged over the parameters θ
k, attained by the pMCMC (left) and the NPMC (right) algorithms, is represented versus the NESS, in the CO and PO scenarios. Simulation runs which attained a final MSE close to the global average value are indicated with big circles (CO) and squares (PO) on both plots. It can be observed that the pMCMC method performs similarly in both scenarios, in terms of MSE and NESS, yielding poor results in both cases. On the contrary, the NPMC method provides significantly better MSE results in the CO scenario, where a larger amount of information is available. The NPMC method does not present degradation due to the high degeneracy occurring in the CO scenario.Figure 5
Performance of the pMCMC (left) and the NPMC (right) methods for the estimation of the whole set of rate parameters θ: MSE (in logarithmic scale) versus the final NESS, for each simulation run in the CO and the PO scenario. The big circles (CO) and squares (PO) represent simulation runs with a final mean MSE close to the global average.Figure6 depicts the evolution along the iterations of the NESS (left) and the MSE (right) averaged over P
=
100 independent simulation runs for the NPMC algorithm. Both indices converge to a steady value in a low number of iterations also in this complex scenario. As expected, a significantly higher final MSE is attained in the extremely data-poor PO scenario.Figure 6
Evolution along the NPMC iterations of the average NESS (left) and MSE (right) in the CO and the PO scenario.In Figure7 (left) the average ACF attained by the pMCMC in the CO and the PO cases is represented. Even after thinning the output, the sample correlation is extremely high in both scenarios, which leads to a very low NESS. The acceptance rate is also very low and very long chains are required to obtain reasonable results. In the PO scenario ≈
44 samples are accepted on average in a simulation run of I
=
15
×
1
0
3 samples (acceptance rate 0.0029). In the CO case, only ≈
23 samples are accepted on average (rate 0.0015).Figure 7
Left. Autocorrelations based on the final sample of size 1
0
3 of the pMCMC scheme in the CO and the PO scenarios, averaged over P
=
100 simulation runs.Right. Markov chain provided by the pMCMC method in the PO scenario, corresponding to the average simulation run depicted with a big square in Figure 5 (left).Figure7 (right) depicts the final Markov chain provided by the pMCMC method (after removing the burn-in period and thinning the output) in the average simulation run represented with a big square in Figure 5 (left). It can be observed that the mixing of the chain is very poor, with a total number of accepted samples of 46 (close to the average). Many other simulations, in both the PO and the CO scenarios, provide even lower number of accepted samples and thus very inconsistent results.Figure8 depicts the final Gaussian approximations of the marginal posteriors p
(
θ
k
∣
y
), k
=
1
,
…
,
8, obtained by the pMCMC and the NPMC methods, in the CO and PO scenarios, for the average simulation runs represented as big circles and squares in Figure 5. We can observe that the NPMC method provides a significantly better approximation of the log-rate parameters in the CO scenario, where a larger amount of data is available, which is also clear from Figure 5 (right). However, the pMCMC on average performs similarly in both scenarios, due to the low efficiency of the pMCMC sampling scheme when the dimension of the problem (either K or N) increases.Figure 8
Marginal posterior pdf approximations of each parameterp
^
(
θ
k
∣
y
), k
=
1
,
…
,
K, attained in an average simulation run by the pMCMC and the NPMC, in the CO and in the PO case.In Table2 the MSE of each parameter θ
k averaged over P
=
100 independent simulation runs is shown, as obtained with the pMCMC and the NPMC schemes, for the CO and the PO experiments. In the CO case, NPMC provides homogeneous results for all parameters. On the contrary, in the PO case, some of the parameters (specially θ
5 and θ
6) are significantly poorly estimated, presenting a final MSE close to the initial value (which corresponds to the prior knowledge). The pMCMC scheme presents significantly higher MSE values than NPMC in both observation scenarios and for all parameters θ
k.Table 2
Final MSE for the parametersθ
k
,
k
=
1
,
…
,
K in the CO and PO experiments, averaged over the simulation runs. The last two columns correspond to the mean and standard deviation (std) values of the global MSE (averaged over the parameters). The prior values are included for comparison.
θ
1
θ
2
θ
3
θ
4
θ
5
θ
6
θ
7
θ
8
mean MSE
std MSE
Prior
6.789
11.344
8.853
7.543
6.789
12.484
8.430
6.789
8.628
0
PO
pMCMC
3.412
3.319
5.543
3.200
7.059
8.929
6.799
4.371
5.329
2.926
NPMC
1.246
1.011
2.214
1.490
4.073
7.015
2.311
1.856
2.652
1.020
CO
pMCMC
2.899
2.958
1.676
1.572
1.604
1.547
1.573
1.468
1.912
1.476
NPMC
0.305
0.302
0.162
0.167
0.280
0.280
0.156
0.168
0.228
0.091Figure9 depicts the population posterior mean x
^
=
E
p
(
x
∣
y
)
[
x
] corresponding to the average simulation runs of the pMCMC and the NPMC methods in the PO scenario, represented as big squares in Figure 5. Again, the NPMC method provides more accurate estimates of the unobserved populations than the pMCMC method, especially for x
R
N
A. In the CO scenario both methods provide good approximations of the populations of all species.Figure 9
Posterior meanx
^
=
E
p
(
x
∣
y
)
[
x
] of the populations of all species obtained in the average simulation run of the pMCMC (left) and the NPMC (right) schemes, in the PO scenario.


## 5.4.1. Estimation of a Single Rate Parameterθ
1

In this section we present numerical results regarding the approximation of the posterior distributionp
(
θ
1
,
x
∣
θ
∖
1
,
y
) of a single rate parameter θ
1
=
log
⁡
c
1 and the populations x, when the rest of parameters θ
∖
1
=
[
θ
2
,
…
,
θ
K
]
⊤, are assumed known.We compare the pMCMC and the NPMC methods in this simple scenario in order to show that the two algorithms perform almost equivalently in low dimensional problems. This is in a clear contrast with the results presented in Section5.4.2, which show that the NPMC method can be more efficient in higher dimensional settings.We have performedP
=
100 independent simulation runs of the pMCMC and the NPMC schemes in the CO and the PO scenarios, with independent population and observation vectors in each simulation p
=
1
,
…
,
P. Both in the CO and the PO cases, the same true population trajectories x
(
p
) were used, but the observations in the CO scenario, y
C
O
(
p
), and in the PO scenario, y
P
O
(
p
), differ. The number of observation times has been set to N
=
100 and exactly the same data has been used for NPMC and pMCMC.Following [5], as a proposal pdf q
(
θ
⋆
∣
θ
) in the pMCMC scheme we consider a Gaussian random walk update with variance γ
2. We have performed simulations with different γ values and the best results were obtained with γ
2
=
1, which are presented here. A total number of I
=
1
0
4 iterations have been run in each simulation. A final sample of size M
=
1
0
3 has been obtained from each Markov chain by discarding a burn-in period of 1
0
3 samples and thinning the output by a factor of 9.In the NPMC scheme, the number of iterations has been set toL
=
10, the number of samples per iteration is M
=
1
0
3 and theclipping parameter is M
T
=
100. In this way, the computational effort of the two methods is approximately the same, as they both generate I
=
M
L
=
1
0
4 samples in the space of θ.In Figure1 the final MSE obtained by the pMCMC (left) and the NPMC (right) algorithms for each simulation run is depicted versus the final NESS, in the CO and the PO scenarios. Note that the NESS is computed differently for pMCMC and NPMC. It can be observed that both algorithms perform similarly in this case, with an equivalent computational cost. Both algorithms attain on average lower MSE values in the CO scenario, as expected. However, the NESS also takes lower values in the CO case, which indicates a worse mixing of the Markov chains in the pMCMC algorithm and also higher degeneracy in the NPMC algorithm.Figure 1
Performance of the pMCMC (left) and the NPMC (right) methods for the estimation of a unique rate parameter θ
1: MSE (in logarithmic scale) obtained from the final output versus the NESS for each simulation run in the CO and the PO scenario. The big circles (CO) and squares (PO) represent simulation runs with a final mean MSE close to the global average.In Figure2 the evolution of the NESS (left) and the MSE (right) along the iterations of the NPMC algorithm is represented, for the CO and the PO scenarios. It can be observed that both measures attain a steady value by the 5-th iteration, both in the CO and the PO case, which suggests that actually less iterations are sufficient for this problem. Again, we observe that in the CO scenario both the NESS and the MSE reach lower values.Figure 2
Evolution along the iterations of the NPMC algorithm of the average NESS (left) and MSE (right) in the CO and PO scenarios, estimating a single parameter θ
1.Figure3 (left) plots the average ACF of the final pMCMC sample, after removing the burn-in period and thinning the Markov chain by a factor of 9. Particularly high correlations are present in the CO case, leading to a poor NESS. Related to the ACF, the average sample acceptance probability in the pMCMC scheme in the PO scenario is 0.091, while in the CO scenario it is only 0.0034, which means that 910 samples are accepted out of I
=
1
0
4 in the CO case and only 34 in the CO case.Figure 3
Left. Average ACF based on the final sample of size M
=
1
0
3 of the pMCMC scheme in the CO and the PO scenarios, averaged over P
=
100 simulation runs.Right. Marginal posterior pdf estimates p
^
(
θ
1
,
∣
θ
∖
1
,
y
) of an average simulation run, for pMCMC and NPMC in the CO and PO scenarios. The true value θ
1 is also shown.In Figure3 (right) the final pdf estimates p
^
(
θ
1
∣
θ
∖
1
,
y
) of the average simulation runs represented as big circles and crosses in Figure 1 are represented in the CO and the PO scenario, for the pMCMC and the NPMC schemes. For the pMCMC method we have built a Gaussian approximation of the posterior density p
(
θ
1
∣
θ
∖
1
,
y
) based on the final MCMC sample {
θ
1
(
i
)
}
i
=
1
M. For the NPMC method, this approximation corresponds to the proposal pdf for the next iteration L
+
1, i.e., p
^
(
θ
1
∣
θ
∖
1
,
y
)
=
q
L
+
1
(
θ
1
)
=
N
(
θ
1
;
μ
L
,
1
,
σ
L
,
1
2
), where the mean and variance terms μ
L
,
1 and σ
L
,
1
2 are computed as in ∗
. It can be observed in Figure 3 (right) that very similar results are obtained by both algorithms in this scenario. The final MSE values obtained by the pMCMC and the NPMC methods, averaged over P
=
100 simulation runs, are shown in Table 1, together with the MSE corresponding to the prior distribution.Table 1
Final mean and standard deviation (std) values of the MSE forθ
1 in the CO and PO scenarios, for pMCMC and NPMC. The prior values are included for comparison.
mean MSE
std MSE
Prior
6.789
0
PO
pMCMC
0.215
0.171
NPMC
0.195
0.170
CO
pMCMC
0.027
0.026
NPMC
0.022
0.016Figure4 depicts the posterior mean of the populations, x
^
=
E
p
^
(
x
∣
y
)
[
x
], obtained with pMCMC (left) as x
^
=
1
/
M
∑
i
=
1
M
x
(
i
) and with NPMC (right) as x
^
=
∑
i
=
1
M
w
¯
L
(
i
)
x
L
(
i
) in the PO scenario. The results correspond to the particular simulation runs (different for pMCMC and NPMC) identified with big squares in Figure 1 and whose posterior approximations, p
^
(
θ
1
∣
θ
∖
1
,
y
), are shown in Figure 3 (right). It can be observed that, in the PO scenario, the trend of the population of all the species is reasonably identified, even though only a linear combination of the proteins is observed. In the CO scenario the populations of all species are accurately estimated and are not shown for conciseness. Note that the populations of all species are very low, which suggests that the diffusion approximation may perform poorly in this scenario.Figure 4
Posterior mean,x
^
=
E
p
^
(
x
∣
y
)
[
x
], of the populations obtained in a single simulation run of pMCMC (left) and NPMC (right) in the PO scenario (only a linear combination of the proteins is observed, corrupted by noise).The results presented in this section reveal a very similar performance of the two methods in this simple scenario, which suggests that the parameters of pMCMC have been properly tuned. Also in terms of computational complexity pMCMC and NPMC perform very similarly. The execution time per1
0
3 samples (one NPMC iteration and 1
0
3 pMCMC iterations) for the pMCMC scheme is 312 seconds, while for NPMC it is 325 seconds, both in the CO and in the PO cases, on a 3-GHz Intel Core 2 Duo CPU, with 4 GB of RAM. The stochastic forward simulation of the prokaryotic model with the Gillespie algorithm has been implemented in C, and the rest of the code in Matlab R2014b. Out of the execution time per iteration of the NPMC method, only 0.06 second correspond to batch steps 4, 5, and 6, and the rest of operations can be easily parallelized, as previously discussed.However, the pMCMC method does not allow for an easy parallelization, provides a set of highly correlated samples (especially in the CO scenario) and requires the setting of the proposal varianceγ
2, the burn-in and the thinning parameters, which may not be straightforward, and determines the performance of the algorithm. On the contrary, the NPMC scheme provides uncorrelated sets of samples at each iteration and does not require the precise fitting of any parameters. Additionally, the computer simulations suggest that the convergence of the NPMC algorithm may be assessed observing the evolution of the NESS, which usually reaches a steady value simultaneously with the MSE.


## 5.4.2. Estimation of All the Parametersθ
k
,
k
=
1
,
…
,
K

In this section we present simulation results to evaluate the performance of the pMCMC and the NPMC schemes in the approximation of the posterior distribution of the rate parameters and the populations of all species,p
(
θ
,
x
∣
y
), assuming that all the rate parameters are unknown, again in the CO and the PO scenarios.In this case,N
=
200 observation times are assumed for all the simulations. Again, P
=
100 independent simulation runs of each algorithm have been performed. The NPMC scheme has been run for L
=
15 iterations, with M
=
1
0
3 samples per iteration andclipping parameter M
T
=
100. The pMCMC scheme has been run with I
=
15
×
1
0
3 iterations in each simulation run, a burn-in period of 1
0
3 iterations, and thinning the output by a factor of 14. With this setup the computational effort is approximately the same in the two schemes. The variance of the random walk proposal has again been set to γ
2
=
1, since this value seems to yield the best performance.In Figure5 the MSE (in logarithmic scale), averaged over the parameters θ
k, attained by the pMCMC (left) and the NPMC (right) algorithms, is represented versus the NESS, in the CO and PO scenarios. Simulation runs which attained a final MSE close to the global average value are indicated with big circles (CO) and squares (PO) on both plots. It can be observed that the pMCMC method performs similarly in both scenarios, in terms of MSE and NESS, yielding poor results in both cases. On the contrary, the NPMC method provides significantly better MSE results in the CO scenario, where a larger amount of information is available. The NPMC method does not present degradation due to the high degeneracy occurring in the CO scenario.Figure 5
Performance of the pMCMC (left) and the NPMC (right) methods for the estimation of the whole set of rate parameters θ: MSE (in logarithmic scale) versus the final NESS, for each simulation run in the CO and the PO scenario. The big circles (CO) and squares (PO) represent simulation runs with a final mean MSE close to the global average.Figure6 depicts the evolution along the iterations of the NESS (left) and the MSE (right) averaged over P
=
100 independent simulation runs for the NPMC algorithm. Both indices converge to a steady value in a low number of iterations also in this complex scenario. As expected, a significantly higher final MSE is attained in the extremely data-poor PO scenario.Figure 6
Evolution along the NPMC iterations of the average NESS (left) and MSE (right) in the CO and the PO scenario.In Figure7 (left) the average ACF attained by the pMCMC in the CO and the PO cases is represented. Even after thinning the output, the sample correlation is extremely high in both scenarios, which leads to a very low NESS. The acceptance rate is also very low and very long chains are required to obtain reasonable results. In the PO scenario ≈
44 samples are accepted on average in a simulation run of I
=
15
×
1
0
3 samples (acceptance rate 0.0029). In the CO case, only ≈
23 samples are accepted on average (rate 0.0015).Figure 7
Left. Autocorrelations based on the final sample of size 1
0
3 of the pMCMC scheme in the CO and the PO scenarios, averaged over P
=
100 simulation runs.Right. Markov chain provided by the pMCMC method in the PO scenario, corresponding to the average simulation run depicted with a big square in Figure 5 (left).Figure7 (right) depicts the final Markov chain provided by the pMCMC method (after removing the burn-in period and thinning the output) in the average simulation run represented with a big square in Figure 5 (left). It can be observed that the mixing of the chain is very poor, with a total number of accepted samples of 46 (close to the average). Many other simulations, in both the PO and the CO scenarios, provide even lower number of accepted samples and thus very inconsistent results.Figure8 depicts the final Gaussian approximations of the marginal posteriors p
(
θ
k
∣
y
), k
=
1
,
…
,
8, obtained by the pMCMC and the NPMC methods, in the CO and PO scenarios, for the average simulation runs represented as big circles and squares in Figure 5. We can observe that the NPMC method provides a significantly better approximation of the log-rate parameters in the CO scenario, where a larger amount of data is available, which is also clear from Figure 5 (right). However, the pMCMC on average performs similarly in both scenarios, due to the low efficiency of the pMCMC sampling scheme when the dimension of the problem (either K or N) increases.Figure 8
Marginal posterior pdf approximations of each parameterp
^
(
θ
k
∣
y
), k
=
1
,
…
,
K, attained in an average simulation run by the pMCMC and the NPMC, in the CO and in the PO case.In Table2 the MSE of each parameter θ
k averaged over P
=
100 independent simulation runs is shown, as obtained with the pMCMC and the NPMC schemes, for the CO and the PO experiments. In the CO case, NPMC provides homogeneous results for all parameters. On the contrary, in the PO case, some of the parameters (specially θ
5 and θ
6) are significantly poorly estimated, presenting a final MSE close to the initial value (which corresponds to the prior knowledge). The pMCMC scheme presents significantly higher MSE values than NPMC in both observation scenarios and for all parameters θ
k.Table 2
Final MSE for the parametersθ
k
,
k
=
1
,
…
,
K in the CO and PO experiments, averaged over the simulation runs. The last two columns correspond to the mean and standard deviation (std) values of the global MSE (averaged over the parameters). The prior values are included for comparison.
θ
1
θ
2
θ
3
θ
4
θ
5
θ
6
θ
7
θ
8
mean MSE
std MSE
Prior
6.789
11.344
8.853
7.543
6.789
12.484
8.430
6.789
8.628
0
PO
pMCMC
3.412
3.319
5.543
3.200
7.059
8.929
6.799
4.371
5.329
2.926
NPMC
1.246
1.011
2.214
1.490
4.073
7.015
2.311
1.856
2.652
1.020
CO
pMCMC
2.899
2.958
1.676
1.572
1.604
1.547
1.573
1.468
1.912
1.476
NPMC
0.305
0.302
0.162
0.167
0.280
0.280
0.156
0.168
0.228
0.091Figure9 depicts the population posterior mean x
^
=
E
p
(
x
∣
y
)
[
x
] corresponding to the average simulation runs of the pMCMC and the NPMC methods in the PO scenario, represented as big squares in Figure 5. Again, the NPMC method provides more accurate estimates of the unobserved populations than the pMCMC method, especially for x
R
N
A. In the CO scenario both methods provide good approximations of the populations of all species.Figure 9
Posterior meanx
^
=
E
p
(
x
∣
y
)
[
x
] of the populations of all species obtained in the average simulation run of the pMCMC (left) and the NPMC (right) schemes, in the PO scenario.


## 6. Conclusion

We have investigated the use of Monte Carlo-based Bayesian computation methods for approximating posterior distributions of the parameters and the species populations in stochastic kinetic models. Specifically, we have applied both particle Markov chain Monte Carlo (pMCMC) and nonlinear population Monte Carlo (NPMC) methods, which rely on different sampling and approximation schemes. Both pMCMC and NPMC methods resort to a sequential Monte Carlo approximation of the posterior populations to estimate the unknown interaction parameters. However, the pMCMC constructs a Markov chain of correlated parameter samples while the NPMC algorithm is based on an importance sampling scheme with nonlinearly transformed weights to avoid degeneracy and the numerical problems typically arising in importance samplers when applied to high dimensional problems.We have compared the performance of the two schemes for a challenging autoregulatory feedback network with 5 species and 8 unknown interaction rate parameters. Both methods have been applied without resorting to approximations of the underlying jump process (e.g., the diffusion approximation). We have shown how the NPMC algorithm outperforms the pMCMC method and requires only a moderate computational cost. Besides, the proposed method has a set of important features, common to all PMC schemes, as the sample independence, ease of parallelization, and compared to MCMC schemes, there is no need for convergence (burn-in) periods.Let us note that there is room for improvement of both algorithms in practical applications. To be specific, the pMCMC scheme we have assessed is constructed around a random walk Metropolis algorithm, but other MCMC schemes can be adopted; hence it is possible to further optimize its performance. The same is true for the NPMC method, however. In this case, the step to optimize in practice is the choice of the importance function. In our comparison we have used multivariate Gaussian proposals (because they are straightforward to fit) but other, more efficient proposals can indeed be used within the proposed framework. Also note that a computationally simpler approach to inference in some SKMs is the use of sequential Monte Carlo approximate Bayesian computation (SMC-ABC) techniques [26, 27]. In the SMC-ABC framework the state-space model (for x
n) is only used for forward simulation. There is no attempt to approximate the likelihood of θ, which is replaced by some simple summary statistic. SMC-ABC methods are attractive when computational cost is a major concern, but they can be expected to deliver an inferior performance (as they do not fully exploit the probabilistic model, when it is available).



---

*Source: 7160934-2019-01-20.xml*