# A Point Mass Proposal Method for Bayesian State-Space Model Fitting

Mary Llewellyn, Ruth King, Víctor Elvira and Gordon Ross
School of Mathematics and Maxwell Institute for Mathematical Sciences,
University of Edinburgh, Edinburgh, UK

(January 2023)

###### Abstract

State-space models (SSMs) are often used to model time series data where the
observations depend on an unobserved latent process. However, inference on the
process parameters of an SSM can be challenging, especially when the
likelihood of the data given the parameters is not available in closed-form.
We focus on the problem of model fitting within a Bayesian framework, for
which existing approaches include Markov chain Monte Carlo (MCMC) using
Bayesian data augmentation, sequential Monte Carlo approximation, and particle
MCMC algorithms, which combine sequential Monte Carlo approximations and MCMC
steps. However, these different methods can be inefficient when sample
impoverishment occurs during the sequential Monte Carlo approximation and/or
when the MCMC algorithm mixes poorly. We propose the use of deterministic
hidden Markov models (HMMs) to provide an efficient MCMC with data
augmentation approach, imputing the latent states within the algorithm. Our
approach deterministically approximates the SSM by a discrete HMM, which is
subsequently used as an MCMC proposal distribution for the latent states in
Metropolis-within-Gibbs steps. We demonstrate that the algorithm provides an
efficient alternative method for state-space models with near-chaotic
behaviour.

##  1 Introduction

State-space models (SSMs) describe time series observations as dependent on an
unobserved system process (Durbin and Koopman,, 2012). The unobserved system
process of an SSM consists of a sequence of continuously-valued latent states
that evolve through time but are not observed directly. Instead, observations
are modelled as a function of these unobserved latent states via the
observation process. SSMs have been applied to various problems, including the
modelling of inflation rates (Koopman and Bos,, 2004), neuron responses (Smith
and Brown,, 2003; Lin et al.,, 2019), animal movement (Patterson et al.,,
2008; Auger-Méthé et al.,, 2016), and population dynamics (Buckland et al.,,
2004; King,, 2012; Newman et al.,, 2014).

Fitting an SSM usually refers to inference on the model parameters given the
observations. However, outside of the discrete-state hidden Markov model (HMM)
case (Rabiner,, 1989), or special continuous cases, for example, linear
Gaussian SSMs (Kalman,, 1960; Durbin and Koopman,, 2012), a closed-form
expression for the likelihood does not exist, leading to model fitting
challenges. The various approaches to approximating the likelihood include the
use of linear Gaussian approximations (Julier and Uhlmann,, 1997; Wan and Van
Der Merwe,, 2000), numerical methods using Laplace or HMM approximations (Bucy
and Senne,, 1971; Koyama et al.,, 2010; Langrock and King,, 2013; Thygesen et
al.,, 2017; Herliansyah et al.,, 2022), the use of sequential Monte Carlo to
obtain an unbiased approximation to the likelihood (Beaumont,, 2003; Andrieu
and Roberts,, 2009; Kantas et al.,, 2015; Deligiannidis et al.,, 2018), and
Bayesian data augmentation approaches (Frühwirth-Schnatter,, 2004; Fearnhead,,
2011; Borowska and King,, 2023). A review of these methods, and others, is
given in Newman et al., (2022).

Here, we focus on Bayesian data augmentation approaches, which can be applied
to general SSMs by specifying the latent states as additional auxiliary
variables. In doing so, the data augmentation approach forms the joint
posterior distribution of the model parameters and additional latent states
(Tanner and Wong,, 1987; Hobert,, 2011). The resulting joint (complete data)
likelihood of the parameters and states is typically available in closed-form,
permitting the application of standard MCMC algorithms, such as the Gibbs
sampler, Metropolis-Hastings algorithms, or Metropolis-within-Gibbs
algorithms. An estimate of the marginal posterior distribution of the model
parameters can then be obtained by simply retaining the sampled parameter
values and disregarding the latent state samples to effectively marginalise
over the states. It is, however, typically challenging to efficiently sample
the latent state vector within an Metropolis-Hastings algorithm, since there
are typically many highly correlated latent states, leading to poor mixing
(Frühwirth-Schnatter,, 2004). Updating the states one at a time typically
performs very poorly due to the high correlation between states (Fearnhead,,
2011), whilst designing good block-update proposal distributions can be
challenging (Shephard and Pitt,, 1997).

Historically, data augmentation has been combined with numerical
approximation, sequential Monte Carlo, and/or MCMC to fit the parameters of an
SSM. When an SSM is well-approximated by Gaussian distributions, efficient
proposal distributions for the latent states in a data augmentation approach
include those derived from Laplace approximation (Kristensen et al.,, 2016)
and extensions to the Kalman filter (Giordani and Kohn,, 2010). In the general
case, more computationally intensive methods are usually required, such as the
particle Gibbs (PG) sampler (Lindsten et al.,, 2012; Chopin and Singh,, 2015;
Fearnhead and Meligkotsidou,, 2016) or particle Gibbs with ancestor sampling
(PGAS) (Lindsten et al.,, 2014). PG and PGAS are classes of particle MCMC
method (Andrieu et al.,, 2010) which use sequential Monte Carlo (SMC) to
update all states simultaneously conditional on the parameters. In their use
of SMC, PG and PGAS methods rely on few assumptions about the underlying state
and observation distributions. However, these methods can be computationally
expensive, particularly if many samples are required in the SMC steps due to
sample impoverishment (Rainforth et al.,, 2016; Wang et al.,, 2017). One
approach to overcoming sample impoverishment is the introduction of
deterministic resampling, which uses a deterministic grid over the state-space
to retain samples that may otherwise be discarded in the SMC steps (Li et
al.,, 2012). Although we do not develop a particle MCMC approach here, the use
of a deterministic grid to efficiently sample the latent states is related to
the approach presented in this paper.

We update the model parameters conditional on the imputed latent states, and
we update the latent states conditional on the model parameters in Metropolis-
within-Gibbs steps. We focus on the latter and propose the point mass proposal
Metropolis-Hastings (PMPMH) algorithm, formulating proposals for the latent
states in two steps. The first step uses a discretization of the state space
into a set of pre-defined intervals (a grid) to reduce the SSM to an HMM. We
then sample a sequence of intervals from the approximate distribution of the
states given the model parameters, using forward-filtering backward-sampling
(FFBS) algorithm (Rabiner,, 1989). In the second stage of the algorithm, we
then correct the discretization error imposed by the approximate HMM by
sampling from within the intervals using a continuous proposal distribution,
using these samples to update blocks of states in Metropolis-Hastings steps.
The model parameters are updated conditionally on the latent states using
standard Gibbs or Metropolis-within-Gibbs updates. The first step of the
proposed algorithm is related to point mass smoothing (Bucy and Senne,, 1971;
Kitagawa,, 1987; de Valpine and Hastings,, 2002; Langrock et al.,, 2012;
Langrock and King,, 2013) and embedded HMMs (Neal,, 2003; Shestopaloff and
Neal,, 2013; Finke et al.,, 2016; Shestopaloff and Neal,, 2018), which use the
tractable discrete case of an SSM, an HMM, to sample latent states. However,
here we use the embedding of an HMM via deterministic grid cells to ensure
latent states are proposed with reasonable posterior mass.

The rest of the paper is structured as follows. We start by defining SSMs and
HMMs, and model fitting procedures for SSMs via Bayesian data augmentation in
Section 2. We describe the proposed data augmentation scheme in Section 3, the
PMPMH algorithm, in which we use a grid-based HMM approximation within a
proposal distribution for the latent states. In Section 4, we provide three
particular cases of the algorithm, each given by a different method for
defining the grid cells, and describe when each can be usefully applied.
Finally, we illustrate our approach with two case studies in Section 5 and
discussed the proposed algorithm and extensions in Section 6.

##  2 State-space model inference

Here we define state-space models, the notation used throughout the paper, and
the special case of discrete-valued system processes: HMMs. We then discuss
the challenges of fitting state-space models to data and computational
approaches in Section 2.1, motivating our proposed data augmentation
algorithm.

We consider time series data observed at discrete (regular) time points up to
time $T$, $y_{1:T}=(y_{1},\dots,y_{T})$. An SSM models these data via (_i_) a
system process of unobserved (latent) states, $x_{1:T}=(x_{1},\dots,x_{T})$;
and (_ii_) an observation process linking the latent process with the observed
$y_{1:T}$. We allow the states to be continuous, taking values in some set
$x_{t}\in\chi$ for each $t$, and describe their evolution using initial state
and transition distributions. Similarly, the observation process is described
using an observed state distribution. As standard, we assume a first-order
Markov process for the system process and that the observed data are
conditionally independent given the underlying latent states, that is, the
observation process is only a function of the latent state at time $t$. Given
the set of model parameters, collectively denoted by $\boldsymbol{\theta}$,
the SSM can be described mathematically by

$\displaystyle\mbox{\emph{Initial state distribution}:}\qquad p(x_{1}$
$\displaystyle\mid\boldsymbol{\theta}),$ $\displaystyle\mbox{\emph{State
transition distribution}:}\qquad p(x_{t}$ $\displaystyle\mid
x_{t-1},\boldsymbol{\theta}),\hskip 5.69054ptt=2,\dots,T,$
$\displaystyle\mbox{\emph{Observed state distribution}:}\qquad p(y_{t}$
$\displaystyle\mid x_{t},\boldsymbol{\theta}),\hskip 5.69054ptt=1,\dots,T.$

A particular case of an SSM is an HMM, in which the state space is discrete
and finite. In other words, the states at time $t$, $x_{t}$, take values in
the set $\\{1,\dots,N\\}$. Therefore, when the SSM considered is an HMM, we
simply re-define the initial state and transition probabilities for
$k,n=1,\dots,N$ as

Initial state probabilities: $\displaystyle\qquad
P(x_{1}=n\mid\boldsymbol{\theta}),\hskip 5.69054pt$ State transition
probabilities: $\displaystyle\qquad P(x_{t}=n\mid
x_{t-1}=k,\boldsymbol{\theta}),\hskip 5.69054ptt=2,\dots,T.$

###  2.1 Model fitting

We assume that we are primarily interested in inference on the model
parameters, $\boldsymbol{\theta}$, although this may also extend to the latent
states, $x_{1:T}$, depending on the application. In general, inference on the
posterior distribution for $\boldsymbol{\theta}$ requires a closed-form
expression for the observed data (or marginal) likelihood,
$p(y_{1:T}\mid\boldsymbol{\theta})$. The joint likelihood of the model
parameters and latent states is given by

$\displaystyle p(x_{1:T},y_{1:T}\mid\boldsymbol{\theta})$
$\displaystyle=p(x_{1}\mid\boldsymbol{\theta})\prod_{t=2}^{T}p(x_{t}\mid
x_{t-1},\boldsymbol{\theta})\prod_{t=1}^{T}p(y_{t}\mid
x_{t},\boldsymbol{\theta}).$ (1)

A closed-form expression for the marginal likelihood typically only exists in
the linear Gaussian case, where $p(y_{1:T}\mid\boldsymbol{\theta})$ can be
calculated using the Kalman filter (Kalman,, 1960; Durbin and Koopman,, 2012),
or when the state space is discrete and finite (an HMM), where the
marginalisation of Equation (1) amounts to summation over the state space
(Rabiner,, 1989). In general, however, the marginal likelihood is intractable
for continuous SSMs, and standard MCMC implementations cannot be used. One
solution is to use pseudo-marginal likelihood methods (Beaumont,, 2003;
Andrieu and Roberts,, 2009), replacing the marginal likelihood with an
unbiased estimate to formulate a valid Metropolis-Hastings approach, but we do
not focus on these here.

Instead, we consider a Bayesian data augmentation approach. In this case, we
form the joint posterior distribution over the model parameters and latent
states (Tanner and Wong,, 1987; Hobert,, 2011),
$p(x_{1:T},\boldsymbol{\theta}\mid y_{1:T})\propto
p(\boldsymbol{\theta})p(x_{1:T},y_{1:T}\mid\boldsymbol{\theta})$, thus using
the joint likelihood in Equation (1) directly. Since, in general, a closed-
form expression exists for the joint likelihood term up to proportionality,
the proportional likelihood can be evaluated, and standard MCMC approaches can
be used targeting $p(x_{1:T},\boldsymbol{\theta}\mid y_{1:T})$. Samples from
the (marginal) posterior distribution of the parameters,
$p(\boldsymbol{\theta}\mid y_{1:T})$, are then obtained by disregarding the
state samples and retaining the parameter samples.

In order to sample from $p(x_{1:T},\boldsymbol{\theta}\mid y_{1:T})$, a common
approach is to sample the latent states and parameters in turn via their
conditional distributions, simplifying the parameter updates (Carter and
Kohn,, 1994; Frühwirth-Schnatter,, 2004; Durbin and Koopman,, 2012). A
Metropolis-within-Gibbs algorithm can be applied to target the joint
distribution: in each iteration, we update the latent states conditional on
the model parameters; and the model parameters conditional on the latent
states (for which we now have a closed conditional likelihood expression).
However, designing an efficient proposal distribution for the underlying state
vector can be challenging. Ideally, we would be able to design a proposal
distribution to update all states simultaneously and minimise the correlation
between consecutive samples. In general, this state distribution is
intractable, and designing an efficient approximation to the state
distribution can be challenging since it is often high-dimensional in the
number of temporal states. One approach updates all states simultaneously
using a particle Gibbs algorithm (Chopin and Singh,, 2015; Murphy and
Godsill,, 2016; Lindsten et al.,, 2014); values for the entire state vector
are proposed at each iteration from a conditional sequential Monte Carlo
algorithm. In this approach, proposed values for the entire state vector are
always accepted and samples from the correct posterior distribution can be
obtained. However, particle Gibbs methods can suffer from poor mixing sample
impoverishment: degenerating particles being disregarded upon resampling,
impacting mixing and convergence (Wang et al.,, 2017; Rainforth et al.,,
2016). Instead of updating all states simultaneously, for example via particle
Gibbs updates, an alternative approach is to use a single-update algorithm:
updating each latent state individually to reduce the required dimension of
the proposal distribution, but this approach typically leads to very poor
mixing due to correlation in the state process (Fearnhead,, 2011; King,,
2011).

A compromise between single and global updates of the latent states is to use
block updates, simultaneously updating $\ell$ consecutive states (Shephard and
Pitt,, 1997). The posterior conditional distribution of the states in a given
block is, generally, of non-standard form, thus it can be challenging to
define an efficient proposal distribution that accounts for the correlation
between the states. We propose the use of the discrete analogy to SSMs, HMMs,
to define an informative and efficient proposal distributions for blocks of
latent states.

##  3 Point mass proposal Metropolis-Hastings

In this section, we propose the point mass proposal Metropolis-Hastings
(PMPMH) algorithm. The setting is a Metropolis-within-Gibbs algorithm for data
augmentation. The model parameters, $\boldsymbol{\theta}$, are updated
conditionally on $x_{1:T}$ using a standard Metropolis-Hastings or Gibbs
update, and we propose the PMPMH algorithm to sample Metropolis-Hastings
proposals for the latent states given the model parameters.

Step 1 of the algorithm initially converts the SSM into an approximate HMM by
discretizing the state-space at each time point, forming a deterministic grid.
Grid cells are then sampled from the associated (discretized) posterior
distribution, conditional on the model parameters, using the standard forward-
filtering backward-sampling (FFBS) algorithm for HMMs (Rabiner,, 1989) at a
computational cost of $\mathcal{O}(NT)$ for $N$ grid cells at each time point.
Given the sampled grid cells, Step 2 proposes values for $x_{1:T}$ from within
the cells using some specified (bounded) continuous proposal distribution.

For pedagogical purposes, we make a few assumptions in our initial description
of the algorithm, for example, we derive the algorithm for one-dimensional
state-spaces. Extensions to higher-dimensional state-spaces are possible, and
we demonstrate the algorithm in the two-dimensional case with an example in
Section 5. However, we note that there may be practically increased
computational demands for higher-dimensional state-spaces. For notational
simplicity, we present the algorithm for SSMs where the state-space and the
number of grid cells are the same for each time point, that is, the state-
space at all time points is denoted by $\chi$. However, the generalisation to
state-spaces that vary over time is immediate by adapting the notation.
Finally, we initially describe the algorithm in terms of updating the full set
of latent states, with a version of the algorithm for block updates of the
states described later.

###  3.1 Step 1: sampling a grid cell trajectory

In this initial step, we formulate a discrete HMM approximation to the SSM and
sample a trajectory, i.e., a value for $x_{1:T}$, from the discrete
approximation. We then correct for the approximate sampling distribution using
a Metropolis-Hastings step targeting the correct posterior distribution.
First, we propose a partition of the whole state space into a grid: in the
simplest case considered here, the state space at each time $t$, $\chi$, is
partitioned into intervals, forming grid cells. The $N$ grid cells at time $t$
are denoted $I_{t}(n),$ $n=1,\dots,N$ and span $\chi$ with no overlap. The
respective outer grid cells have infinite length if the state space is
unbounded. An example of such a partition is shown in Figure 1.

Figure 1: Example partition of a state space $\chi^{T}$ for each time $t$.
Each grid cell is labelled by the interval it covers.

The grid cells discretize the state space. Thus, the grid cell indices are the
discrete states of an HMM. For the grid cells given by $\\{I_{t}(1),\dots
I_{t}(N)\\}$, the indices are simply $\\{1,\dots N\\}$, with dynamics given by
the SSM density of each grid cell. Mathematically, let $B_{t}$ denote the
random variable of the grid cell indices under the specified grid cell
boundaries so that $B_{t}\in\\{1,\dots,N\\}$ for $t=1,\dots,T$. Then the HMM
probabilities are given for states $k,n\in\\{1,\dots,N\\}$ as

Initial state probabilities: $\displaystyle\qquad
P(B_{1}=n\mid\boldsymbol{\theta}),$ State transition probabilities:
$\displaystyle\qquad P(B_{t}=n\mid B_{t-1}=k,\boldsymbol{\theta}),\hskip
5.69054ptt=2,\dots,T,$ Observed state distribution: $\displaystyle\qquad
p(y_{t}\mid B_{t}=n,\boldsymbol{\theta}),\hskip 5.69054ptt=1,\dots,T.$

These HMM probabilities are therefore, in general, defined by analytically
intractable integrals of the state space equations over the given intervals:

$\displaystyle P(B_{1}=n\mid\boldsymbol{\theta})$
$\displaystyle=\int_{I_{1}(n)}p(x_{1}\mid\boldsymbol{\theta})dx_{1},$
$\displaystyle P(B_{t}=n\mid B_{t-1}=k,\boldsymbol{\theta})$
$\displaystyle=\int_{I_{t}(n)}\int_{I_{t-1}(k)}p(x_{t}\mid
x_{t-1},\boldsymbol{\theta})dx_{t-1}dx_{t},\hskip 5.69054pt2\leq t\leq T,$
$\displaystyle p(y_{t}\mid B_{t}=n,\boldsymbol{\theta})$
$\displaystyle=\int_{I_{t}(n)}p(y_{t}\mid
x_{t},\boldsymbol{\theta})dx_{t},\hskip 5.69054pt1\leq t\leq T.$ (2)

We approximate these integrals using simple (fast) deterministic methods, for
example, the midpoint rule or other Riemann sum methods. Mathematically, we
define $L_{t}(n)$ as the length of the interval $I_{t}(n)$, with the length of
any infinite grid cells set at some arbitrary (finite) length (discussed
further in Appendix B). Within each grid cell, we choose a set of node points,
$\\{\xi_{t}^{s}(n)\\}_{s=1}^{S}$. The set of nodes is defined _a priori_ and
is used to approximate the density functions across each grid cell by a
polynomial. In the simple case where $S=1$ (that is, one node in each grid
cell such as the mid-point), we have that
$\\{\xi_{t}^{s}(n)\\}_{s=1}^{S}=\xi_{t}(n)$, and the HMM probabilities in
Equations (2) for grid cells indices $k,n=1,\dots,N$ are approximated by

$\displaystyle\hat{P}(B_{1}=n\mid\boldsymbol{\theta})$ $\displaystyle\propto
L_{1}(n)p(\xi_{1}(n)\mid\boldsymbol{\theta}),$
$\displaystyle\hat{P}(B_{t}=n\mid B_{t-1}=k,\boldsymbol{\theta})$
$\displaystyle\propto
L_{t}(n)L_{t-1}(k)p(\xi_{t}(n)\mid\xi_{t-1}(k),\boldsymbol{\theta}),\hskip
5.69054ptt=2,\dots,T,$ $\displaystyle\hat{P}(y_{t}\mid
B_{t}=n,\boldsymbol{\theta})$ $\displaystyle\propto
L_{t}(n)p(y_{t}\mid\xi_{t}(n),\boldsymbol{\theta}),\hskip
5.69054ptt=1,\dots,T,$

each appropriately normalised so that they sum to one and are thus valid
probability mass functions. Once the probability mass function approximations
have been obtained, grid cell indices are proposed from

$\displaystyle\hat{P}(B_{1:T}\mid y_{1:T},\boldsymbol{\theta})$
$\displaystyle\propto\hat{P}(B_{1}\mid\boldsymbol{\theta})\hat{P}(y_{1}\mid
B_{1},\boldsymbol{\theta})\prod_{t=2}^{T}\hat{P}(B_{t}\mid
B_{t-1},\boldsymbol{\theta})\hat{P}(y_{t}\mid B_{t},\boldsymbol{\theta}).$

We sample indices from this distribution, denoted $b^{\prime}_{1:T}$, using
the forward-filtering backward-sampling (FFBS) algorithm in Rabiner, (1989).

###  3.2 Step 2: sampling a point trajectory

Step 2 proposes continuous values for the states which will be subsequently
corrected in Metropolis-Hastings steps. Given the discrete grid cell indices
sampled previously, we propose values for $x_{1:T}$ conditional on the
corresponding intervals
$\\{I_{1}(b^{\prime}_{1}),\dots,I_{T}(b^{\prime}_{T})\\}$. This amounts to
sampling $x_{1:T}$ from the domain $I_{1}(b^{\prime}_{1})\times
I_{2}(b^{\prime}_{2})\times\cdots\times I_{T}(b^{\prime}_{T})$. For
simplicity, we consider proposal distributions for each $x_{t}$ independently
and independently of $\boldsymbol{\theta}$, giving proposal distributions of
the form $q(x_{t}\mid x_{t}\in I_{t}(b^{\prime}_{t}))=q(x_{t}\mid
B_{t}=b^{\prime}_{t})$, for $t=1,\dots,T$. For example, a natural choice is to
propose each $x_{t}$ from a uniform proposal distribution over
$I_{t}(b^{\prime}_{t})$ for finite $I_{t}(b^{\prime}_{t})$, or a truncated
Gaussian distribution if $I_{t}(b^{\prime}_{t})$ has infinite length. However,
any distribution with domain $I_{t}(b^{\prime}_{t})$ is applicable.

###  3.3 The PMPMH proposal distribution

The process of sampling a grid cell trajectory and sampling a value for
$x_{1:T}$ defines an independent proposal distribution on $\chi^{T}$. The
density function of the proposed trajectories is therefore the probability of
selecting the grid cell indices from a specified set of grid cell boundaries,
combined with the density of the (continuous) values within those grid cells.
Let $B_{1:T}$ denote the random variables of the grid cell indices under the
proposed grid, which may change across iterations with intervals
$I_{t}(n),\hskip 5.69054ptn=1,\dots,N,\hskip 5.69054ptt=1,\dots,T$. Then the
density of the subsequently proposed state values, $x^{\prime}_{1:T}$, is
given by

$q(x^{\prime}_{1:T}\mid
y_{1:T},\boldsymbol{\theta})=\hat{P}(B_{1:T}=b^{\prime}_{1:T}\mid
y_{1:T},\boldsymbol{\theta})\prod_{t=1}^{T}q(x^{\prime}_{t}\mid
B_{t}=b^{\prime}_{t}).$ (3)

The resulting proposed trajectory, $x^{\prime}_{1:T}$, is retained according
to the Metropolis-Hastings acceptance probability, given by

$p(x_{1:T},x^{\prime}_{1:T}\mid\boldsymbol{\theta})=\min\left(1,\frac{p(x^{\prime}_{1:T}\mid
y_{1:T},\boldsymbol{\theta})q(x_{1:T}\mid
y_{1:T},\boldsymbol{\theta})}{p(x_{1:T}\mid
y_{1:T},\boldsymbol{\theta})q(x^{\prime}_{1:T}\mid
y_{1:T},\boldsymbol{\theta})}\right),$ (4)

where $p(x_{1:T}\mid y_{1:T},\boldsymbol{\theta})$ denotes the posterior
conditional distribution of $x_{1:T}$ given $y_{1:T}$ and
$\boldsymbol{\theta}$, $q(x^{\prime}_{1:T}\mid y_{1:T},\boldsymbol{\theta})$
is as in Equation (3), and $q(x_{1:T}\mid y_{1:T},\boldsymbol{\theta})$ uses
the same definition for the grid cells, that is, their indices, $b_{1:T}$, are
such that $x_{t}\in I(b_{t})$ for all $t$.

The PMPMH method for obtaining a sample of $x_{1:T}$ conditional on
$\boldsymbol{\theta}$ is summarised in Algorithm 1. For completeness, we also
include the updates of $\boldsymbol{\theta}$ conditional on $x_{1:T}$, which
can be updated using standard Metropolis-Hastings or Gibbs updates conditional
on the current chain value of the states.

Algorithm 1 PMPMH algorithm

Input: initial values $x_{1:T}^{(0)}$, $\boldsymbol{\theta}^{(0)}$ and a
number of iterations $M$. A Gibbs or Metropolis-Hastings sampling scheme to
update $\boldsymbol{\theta}$ from $p(\boldsymbol{\theta}\mid
x_{1:T},y_{1:T})$.

for iterations $m=1,\dots M$ do

1:Update $\boldsymbol{\theta}^{(m)}$ from $p(\boldsymbol{\theta}\mid
x_{1:T}^{(m-1)},y_{1:T})$

2:

3:Define a grid with indices $B_{1:T}$

4:for $t=1,..,T$ do $\triangleright$ (Step 1)

5: for $k=1,...,N$ do

6: calculate $\hat{P}(B_{t}=n\mid y_{1:t},\boldsymbol{\theta}^{(m)})$

7:sample $b^{\prime}_{T}$ from $\\{n,\hat{P}(B_{T}=n\mid
y_{1:T},\boldsymbol{\theta}^{(m)})\\}_{n=1}^{N}$

8:for $t=T-1,...,1$ do

9: sample $b^{\prime}_{t}$ from $\\{n,\hat{P}(B_{t}=n\mid
B_{t+1}=b^{\prime}_{t+1},y_{1:t},\boldsymbol{\theta}^{(m)})\\}_{n=1}^{N}$

10:

11:for $t=1,...,T$ do $\triangleright$ (Step 2)

12: sample $x^{\prime}_{t}$ from $q(x_{t}\mid B_{t}=b^{\prime}_{t})$

13:

14:sample $u\sim U(0,1)$

15:if $u<p(x_{1:T}^{(m-1)},x^{\prime}_{1:T}\mid\boldsymbol{\theta})$ then
$\triangleright$ Equation (4)

16: $x_{1:T}^{(m)}\leftarrow x^{\prime}_{1:T}$

17:else $x_{1:T}^{(m)}\leftarrow x^{(m-1)}_{1:T}$

###  3.4 Block updates of the latent states

Updating all of the states simultaneously may be inefficient if there are many
latent states. We can, instead, update the states in smaller blocks: for a
block of $\ell$ states starting at time $t$, we sample an $x_{t:t+\ell-1}$
from a proposal distribution of the form
$q(x_{t:t+\ell-1}\mid x_{t-1},x_{t+\ell},y_{t:t+\ell-1},\boldsymbol{\theta})$.

A block PMPMH proposal distribution is given by a simple adaptation of the
case where all states are updated simultaneously, conditioning on the current
states at either side of the block; $x_{t-1}$ and $x_{t+\ell}$. Therefore, in
Step 1, we condition on the current grid cell indices, $b_{t-1}$ and
$b_{t+\ell}$, such that $x_{t-1}\in I(b_{t-1})$ and $x_{t+\ell}\in
I(b_{t+\ell})$ under the current grid definition. We then use FFBS to sample
proposed indices from the grid cells in the block. In Step 2, we simply define
the proposal distribution, as before, for each latent state in the given
block, $q(x_{s}\mid B_{s}=b^{\prime}_{s})$, for $s=t,\dots,t+\ell-1$. The
PMPMH algorithm for block updates is summarised in Appendix C.

##  4 Defining the grid cells

We have, so far, given a broad framework for using an approximate HMM as an
Metropolis-Hastings proposal distribution. The method fundamentally relies on
the use of a deterministic grid to discretize the state-space. There are,
however, many different ways to define the grid cells (their size and
location), possibly resulting in substantially different proposal
distributions and computational costs. We see in Section 5 that the efficiency
of the algorithm is highly dependent on the choice of grid cell definition and
resulting proposal distribution. Since the optimal choice depends on the
application, we provide the reader with three approaches for defining the grid
cells and describe when each can be usefully (and efficiently) applied. We
discuss the other practical decisions in Appendix B.

###  4.1 Approach 1: equal grid cells

The first approach sets all finite cells equally sized and the same across all
time points, that is, $I_{t}(n)=I_{s}(n)$ for all $t,s=1,\dots,T$ and
$n=2,\dots,N-1$ (assuming an infinite lower and upper bound on the latent
states). We assume that these equally sized grid cells are centred around the
mean of the data, $\frac{1}{T}\sum_{t=1}^{T}y_{t}$, spanning a range of the
state-space denoted by $\mathcal{S}$, though this approach could be adapted by
centring the grid cells around another function of the data.

If the model parameters do not vary with time, only one transition matrix
needs to be calculated per MCMC iteration, or when the model parameters are
updated. For a fixed number of grid cells, this approach has the lowest
computational cost of the approaches considered. However, since the grid cells
are the same for each time point, but the regions of high posterior density
may change at each time point, the grid cells should cover a large range of
the state-space in order to ensure coverage of these high-density regions and
good mixing. This may require many grid cells and a greater computational cost
if the high posterior regions of the SSM are highly non-uniform over time, and
hence this approach is most efficient for SSMs with uniform such regions.

###  4.2 Approach 2: data-driven quantiles

If the high posterior density regions vary over time, a time-inhomogeneous
approximation to the HMM transition matrix may be more efficient than the
previous approach. Hence, we set grid cells in this approach using the
observed data at each time point. In particular, in the implementations of
this paper, we set the grid cell boundaries at the quantiles of one-
dimensional Gaussian distributions centered around each observed data point.
For each time point, $t$, we set the boundaries at the quantiles of $X$ where
$X\sim N(y_{t},{\sigma_{y}}^{2})$. The variance of the Gaussian distribution,
${\sigma_{y}}^{2}$, is set at a scalar value, or as a function of the current
observation process variance within the MCMC iterations and can be used to
control the range that the majority of the grid cells cover via pilot tuning.
The resulting quantiles are simply rounded if integer values are required.
Note that simple extensions to this approach may include using a different
distribution to set the quantiles, provided that its domain is in the state-
space.

When the model parameters are updated, a transition matrix needs to be
recalculated for each time point. It is, however, possible to reduce the
computational cost of this approach by using the same transition matrix across
several time points, centring, for example, on a single data point or the mean
of the corresponding sub-series. By using the observed data at each time
point, we aim to approximately concentrate the grid cells in areas that are
likely to have high posterior mass at each time point. This method may be most
efficient if the observed data used in the centring is a good proxy to the
underlying state process and its dependencies.

###  4.3 Approach 3: latent state quantiles

The previous approaches define the grid cells independently of the state
process, which may result in slow mixing if the latent states of the SSM are
highly correlated over time. Thus, here we set the grid cells using the
current latent states in the MCMC iterations at each time $t$. Let
${x_{t}}^{(m-1)}$ denote the current state in the MCMC iterations at time $t$,
then similarly to the previous approach, we define the grid cell boundaries at
the quantiles of a Gaussian distribution $X$ for each $t$, where $X\sim
N({x_{t}}^{(m-1)},{\sigma_{x}}^{2})$. As with the previous approach, the
variance of the Gaussian distribution used to define the grid cells,
${\sigma_{x}}^{2}$, controls the span of the finite cells. Again, this
variance could be set, for example, as a fixed value (chosen via pilot-
tuning), as a function of the current estimate of the system process variance,
or as a function of the current state within the MCMC iterations. We note that
setting the majority of the grid cells to cover a small range relative to the
high posterior density ranges worked well in practice (see Section 5).

In this approach, the transition matrices need to be recalculated at each time
point for every change in the model parameters but will typically require
fewer grid cells than the other approaches to achieve good mixing.

In general, there is a trade-off between the extent to which the grid cells
are _well-placed_ and the associated computational expense. For example,
equally-sized grid cells are computationally fast but may give lower
acceptance probabilities if they give coarse approximations in high-density
regions. Where equally-sized grid cells give a poor acceptance probabilities,
for example, due to non-uniformity in high posterior regions over time,
acceptance probabilities may be improved using one of the other approaches at
a greater computational expense.

##  5 Numerical illustrations

We demonstrate the proposed point mass proposal Metropolis-Hastings (PMPMH)
algorithm via two case studies. The first is an SSM with a simple one-
dimensional Gaussian mixture state process, demonstrating how the algorithm
can be implemented, and the properties of the algorithm when different
practical decisions are made. We then show how a similar PMPMH implementation
can be used to efficiently sample the latent states of a more challenging
2-dimensional population growth model that can display near-chaotic behaviour.

In each case study, we compare the performance of the algorithm to two
particle Gibbs algorithms: the particle Gibbs (PG) sampler (Chopin and Singh,,
2015; Murphy and Godsill,, 2016) and the particle Gibbs with ancestor sampling
(PGAS) algorithm (Lindsten et al.,, 2014), both with multinomial resampling.
The two methods are particle MCMC algorithms and, like the proposed algorithm,
sample from the joint posterior distribution of the states and parameters by
updating each in turn from their full conditional distributions. The PGAS
algorithm, in particular, is a state-of-the-art method that often improves
upon the mixing properties of the PG algorithm by reducing sample
impoverishment (Nonejad,, 2015; Meent et al.,, 2015; Kantas et al.,, 2015),
though at an increased computational cost. The code used to implement the PG,
PGAS and PMPMH algorithms in each numerical illustration is given in Online
Resource 1 of the Supplementary Information.

###  5.1 Gaussian mixture state-space model

We consider a simple one-dimensional Gaussian mixture state process:

$\displaystyle x_{1}$ $\displaystyle\sim
w_{1}N(1,\sigma^{2}_{\eta_{1}})+(1-w_{1})N(1,\sigma^{2}_{\eta_{2}}),$
$\displaystyle x_{t}\mid x_{t-1}$ $\displaystyle\sim
w_{t}N(x_{t-1},\sigma^{2}_{\eta_{1}})+$
$\displaystyle\quad(1-w_{t})N(x_{t-1},\sigma^{2}_{\eta_{2}}),\hskip
5.69054ptt=2,\dots,T,$ $\displaystyle w_{t}$ $\displaystyle\sim
Bernoulli(p),\hskip 5.69054ptt=1,\dots,T,$

where $p\in[0,1]$ denotes the probability of selecting each mixture component
of the state distribution (the Gaussian distributions with respective
variances $\sigma^{2}_{\eta_{1}}$ or $\sigma^{2}_{\eta_{2}}$). Data,
$y_{1:T}$, are observed according to $y_{t}\mid x_{t}\sim
N(x_{t},\sigma^{2}_{\epsilon})$ and the model parameters are given by
$\boldsymbol{\theta}=(p,\sigma^{2}_{\eta_{1}},\sigma^{2}_{\eta_{2}},\sigma^{2}_{\epsilon})$.
We simulate two data sets from this model, $y_{1:600}^{(1)}$ using
$\boldsymbol{\theta}=(0.9,1,700,1)$ and $T=600$ (Model 1), and
$y_{1:1000}^{(2)}$ using $\boldsymbol{\theta}=(0.99,1,10000,10)$ and $T=1000$
(Model 2), shown in Figure 2.

(a)

(b)

Figure 2: Simulated data from the Gaussian mixture SSM: $y^{(1)}_{1:600}$
using (a) $\boldsymbol{\theta}=(0.9,1,700,1),$ $T=600$, and (b)
$y^{(2)}_{1:1000}$ using $\boldsymbol{\theta}=(0.99,1,10000,10),$ $T=1000$.

Both data sets are simulated with a high-variance second mixture component,
selected with a low probability, resulting in infrequent but large _jumps_ in
the state process. The parameter values for Model 2 result in more occasional
and larger jumps in the state process than the first. We compare the
properties of the algorithm in both cases.

####  5.1.1 PMPMH implementation

We sample from the joint posterior distribution of the states and parameters,
updating $x_{1:T}$ and $\boldsymbol{\theta}$ from their conditional
distributions. The sampling schemes and priors for updating
$\boldsymbol{\theta}$ are given in Appendix A.1. Here, we apply the two stages
of the PMPMH algorithm to update $x_{1:T}$ conditional on
$\boldsymbol{\theta}$. The PMPMH framework presented can be adapted in a
number of ways, for example, by changing the deterministic integration method
used to approximate the HMM probabilities, the number of grid cells, and the
distribution used to propose values from within grid cells. However, we show
in this section that we achieve stable and efficient performance by fixing a
number of these decisions and opting for the simple choices given in Appendix
B.

Instead, we focus on the practical decisions that are relevant to the
efficiency of the algorithm: (a) the choice of each Approach $1-3$ of Section
4 (including the range of the state-space covered by the finite grid cells,
$\mathcal{S}$, which we set using scalar variances), (b) the number of grid
cells, $N$, and (c) the temporal block sizes in which the states are updated,
$\ell$. We test the efficiency of the algorithm under various combinations of
these tuning parameters:

  * (a)

the finite grid cells in Approach 1 over a range of
$\mathcal{S}=150,250,350,450,550$ units for both data sets (compared to the
range of each set of observations, $y_{1:600}^{(1)}$ and $y_{1:1000}^{(2)}$,
equal to $334$ and $478$ units respectively). In Approaches 2 and 3, we use
$\mathcal{S}=1,3,5,7,9,11,13,15$ units for both models,

  * (b)

$N=5,10,20$,

  * (c)

$\ell=1,4,10$, overlapping blocks by one state to improve the mixing of the
states at the “boundaries” of each block as suggested in (Fearnhead,, 2011).

####  5.1.2 Results

For each combination of the tuning parameters listed above, the results are
based on $10$ separate runs of $10,000$ iterations each, taking $1.3-7.5$
hours with $1.6$ GHz CPU for Model 1 (depending on the number of grid cells
and the approach chosen) and $2.2-11$ hours for Model 2. We present the
results for each model in Tables 1 and 2 respectively.

$\mathbf{y^{(1)}_{1:600}}$
---
|  |  $\%$ of range | Approach 2 | Approach 3
N | $\mathbf{\mathcal{S}}$ | of the data | ESS | ESS/s | ESS | ESS/s
$5$ | $1$ | 0.3 | 100 | 0.02 | - | -
| $3$ | 0.9 | 900 | 0.18 | 1000 | 0.20
| $5$ | 1.5 | 700 | 0.14 | 100 | 0.02
| $7$ | 2.1 | 300 | 0.06 | - | -
| $9$ | 2.6 | 200 | 0.04 | - | -
| $\geq 11$ | $\geq 3.3$ | - | - | - | -
$10$ | $1$ | 0.3 | 100 | 0.01 | - | -
| $3$ | 0.9 | 1000 | 0.10 | 1700 | 0.17
| $5$ | 1.5 | 2400 | 0.24 | 200 | 0.02
| $7$ | 2.1 | 2200 | 0.22 | - | -
| $9$ | 2.6 | 2000 | 0.20 | - | -
| $\geq 11$ | $\geq 3.3$ | - | - | - | -
$20$ | $1$ | 0.3 | 100 | 0.00 | - | -
| $3$ | 0.9 | 700 | 0.03 | 2200 | 0.08
| $5$ | 1.5 | 2700 | 0.10 | 1000 | 0.04
| $7$ | 2.1 | 3000 | 0.11 | - | -
| $9$ | 2.6 | 3000 | 0.11 | - | -
| $\geq 11$ | $\geq 3.3$ | - | - | - | -
Table 1: Average effective sample size (ESS) and effective sample size per
second (ESS/s) for Approaches 2 and 3 of the PMPMH algorithm with temporal
blocks of size $\ell=4$. The span of the finite cells is denoted by
$\mathcal{S}$, the number of grid cells, $N$. Dashed lines indicate that
convergence had not occurred within $10,000$ iterations. The ESS is rounded to
the nearest $100$ to account for variation between chains. Since at least one
model parameter was updated at each iteration for each run of the algorithm,
the computational times are the same for Approaches 2 and 3: $5,000$ seconds
for $N=5$, $10,000$ seconds for $N=10$, and $27,000$ seconds for $N=20$, where
the computational time is rounded to the nearest $1000$ seconds to account for
variations in computing speeds.  $\mathbf{y^{(2)}_{1:1000}}$
---
|  |  $\%$ of range | Approach 2 | Approach 3
N | $\mathbf{\mathcal{S}}$ | of the data | ESS | ESS/s | ESS | ESS/s
$5$ | $\leq 3$ | $\leq 0.6$ | - | - | - | -
| $5$ | $1$ | - | - | 200 | 0.03
| $7$ | $1.5$ | - | - | 500 | 0.06
| $\geq 9$ | $\geq 1.9$ | - | - | - | -
$10$ | $\leq 3$ | $\leq 0.6$ | - | - | - | -
| $5$ | $1$ | - | - | 200 | 0.02
| $7$ | $1.5$ | 100 | 0.01 | 500 | 0.04
| $9$ | $1.9$ | 200 | 0.02 | 200 | 0.02
| $11$ | $2.3$ | 100 | 0.01 | - | -
| $\geq 13$ | $\geq 2.7$ | - | - | - | -
$20$ | $\leq 3$ | $\leq 0.6$ | - | - | - | -
| $5$ | $1$ | - | - | 300 | 0.01
| $7$ | $1.5$ | 100 | 0.00 | 600 | 0.02
| 9 | 1.9 | 200 | 0.01 | 400 | 0.01
| 11 | 2.3 | 200 | 0.01 | 400 | 0.01
| $\geq 13$ | $\geq 2.7$ | - | - | - | -
Table 2: Average effective sample size (ESS) and effective sample size per
second (ESS/s) for Approaches 2 and 3 of the PMPMH algorithm with temporal
blocks of size $\ell=4$. The span of the finite cells is denoted by
$\mathcal{S}$, the number of grid cells, $N$. Dashed lines indicate that
convergence had not occurred within $10,000$ iterations. The ESS is rounded to
the nearest $100$ to account for variation between chains. Since at least one
model parameter was updated at each iteration for each run of the algorithm,
the computational times are the same for Approaches 2 and 3: $8,000$ seconds
for $N=5$, $13,000$ seconds for $N=10$, and $38,000$ seconds for $N=20$, where
the computational time is rounded to the nearest $1000$ seconds to account for
variations in computing speeds.

Increasing the number of grid cells, $N$, results in a more accurate
approximation to $p(x_{1:T}\mid y_{1:T},\boldsymbol{\theta})$ in the region
covered by the finite grid cells, improving mixing but at a higher
computational cost. The equally-spaced grid cell Approach 1 performed poorly
on both simulated models compared to the quantile-based methods of Approaches
2 and 3, requiring a large span for both models ($350$ and $550$ units
respectively) due to a large range of values in the high posterior regions,
and thus large numbers of grid cells and a large computational cost, to
provide a reasonable HMM approximation.

Under Approaches 2 and 3, using blocks of size $\ell=10$ or a span for the
finite cells greater than $11$ units required more than $N=20$ grid cells for
convergence and resulted in much more costly and less efficient
implementations than the other sets of tuning parameters. Further, using
single-site updates ($\ell=1$) gave poor mixing when compared to blocks of
size $\ell=4$ due to the correlation between consecutive states. The results
when using $\ell=1,10$, a finite span greater than $11$ units and $N>20$ are
thus excluded from Tables 1 and 2 for brevity.

The results for Approaches 2 and 3 in Table 1 use blocks of size $\ell=4$, and
where convergent, converge within $1000-5000$ iterations using the Brooks-
Gelman-Rubin (BGR) diagnostics in Gelman and Rubin, (1992); Brooks and Gelman,
(1998). In the implementations relating to Model 1, both approaches for
defining the grid cell boundaries yield similar levels of efficiency when
considering the effective sample size (ESS) per second, likely since the small
observation error means that both approaches focus grid cells in roughly the
same region of the state space around the current state. This is potentially
also the reason that both methods required a small span for the finite cells
and a small number of grid cells for convergence relative to the range of the
data ($1-9$ vs. $344$ units, around $0.3-2.6\%$ of the range of the data): as
is true of proposal distributions that make local moves, grid cells can be
focused over a smaller range of the data compared to global-move approaches
and still achieve good acceptance probabilities and mixing. The smaller range
also means that as few as $5$ grid cells can be used to achieve good
convergence properties (compared to $50$ in Approach 1) via a good HMM
approximation over the region, reducing the computational cost of the
approaches.

In comparison to Model 2, where the observation process variance is larger,
Approach 2 is different from a local-move proposal distribution and generally
exhibits poor mixing and convergence in comparison to Approach 3. Conversely,
the local moves of Approach 3 are effective at achieving convergence and
exhibit relatively stable performance even for $N=5$. The efficient ranges for
the finite grid cells are similar to those in the implementations of the first
model as a percentage of the range of the data (here, $1-2.3\%$ of the range
of the data).

We also fitted both of the models using the PG and PGAS algorithms with $5$ to
$1000$ particles, and various combinations of resampling thresholds based on
the standard percentage ESS criterion (Cappé et al.,, 2005). The PG sampler
did not converge for either model using as many as $1000$ particles, resulting
in computational times of around 14 hours for Model 1 and 28 hours for Model 2
on $1.6$ GHz CPU. Conversely, the PGAS sampler converged using as few as $5$
particles for both models. On the whole, the PGAS sampler gave greater levels
of efficiency than the PMPMH algorithm for both of the models, achieving an
average ESS per second of around $4.45$ for Model 1 and $1.49$ for Model 2.
However, this efficiency was not uniform across all states. The average ESS
per second of the states simulated according to the second mixture was $0.25$
to $0.41$ for Model 1 and $0.004$ to $0.008$ for Model 2. In contrast, where
convergent, the PMPMH algorithm was more robust to the mixture associated with
the state, with the ESS per second of states in the second mixture minimally
$99\%$ of those quoted in Tables 1 and 2 (the average ESS per second in the
second mixture-distributed states ranging from $0.01$ to $0.24$ for Model 1
and $0.01$ to $0.06$ for Model 2). We now investigate the comparative
performance of the algorithm on a model that can display near-chaotic
behaviour.

###  5.2 Nicholson’s blowfly model

We consider Nicholson’s Blowfly model for chaotic population growth in Wood,
(2010). The population counts over time, denoted by $N_{1:T}$, arise from two
correlated survival and birth processes, denoted $S_{1:T}$ and $R_{\tau+1:T}$,
$\tau>0$, respectively. Following Wood, (2010), we let
$\exp(-\delta\epsilon_{t})$ denote the daily survival probability with
associated environmental error term $\epsilon_{t}$, such that
$\epsilon_{t}\sim\Gamma(\beta_{\epsilon},\beta_{\epsilon})$,
$\beta_{\epsilon}>0$. The survival component of the system process is given
for $t=1,\dots,T$ by

$S_{t}\sim Binom(N_{t-1},\exp(-\delta\epsilon_{t})),$

Letting $e_{t}$ denote an environmental noise term in the reproductive process
such that $e_{t}\sim\Gamma(\beta_{e},\beta_{e})$, $\beta_{e}>0$, the
reproductive component of the system process is given for $t=\tau+1,\dots,T,$
by

$R_{t}\sim
Po\left(PN_{t-\tau-1}\exp(-\frac{N_{t-\tau-1}}{N_{0}})e_{t}\right),$

where $N_{t}=S_{t}+R_{t}$ for $t=\tau+1,\dots,T$ and $N_{t}=S_{t}$ for
$t=1,\dots,\tau$. We let $\tau=5$ be the known time lag between population
count and subsequent birth count, and $N_{0}=50$ is the known initial
population count. Further, $\epsilon_{1:T}$ and $e_{\tau+1:T}$ are known. The
survival and birth processes, $S_{1:T}$ and $R_{\tau+1:T}$, are unknown latent
states, with observed population counts $y_{t}\sim Po(\phi N_{t})$, for all
$t$. Further, the model parameters
$\boldsymbol{\theta}=(\delta,P,\beta_{\epsilon},\beta_{e},\phi)$ are unknown
and positively valued. Figure 3 shows the simulated data used in this case
study, $y_{1:T}$, using $\boldsymbol{\theta}=(0.7,50,1,0.1,1)$ and $T=300$.

Figure 3: Simulated blowfly population count data using $T=300$ and
$\boldsymbol{\theta}=(0.7,50,1,0.1,1)$.

####  5.2.1 PMPMH implementation

We devise an Metropolis-Hastings scheme targeting the joint posterior
distribution of the unknown latent states and model parameters,
$p(S_{1:T},R_{\tau+1:T},\boldsymbol{\theta}\mid y_{1:T})$. The model
parameters are updated from their full conditional distribution,
$p(\boldsymbol{\theta}\mid S_{1:T},R_{\tau+1:T},y_{1:T})$, using the priors
and sampling schemes given in Appendix A.2. In this section, we describe how
the PMPMH algorithm can efficiently update the latent states, $S_{1:T}$ and
$R_{\tau+1:T}$, conditional on $\boldsymbol{\theta}$.

First, the state process is now two-dimensional at each time point. In order
to simplify the design of the grid used within the PMPMH algorithm, we build
separate PMPMH proposal distributions targeting the full conditional
distributions in each state dimension, $p(S_{1:T}\mid
R_{\tau+1:T},y_{1:T},\boldsymbol{\theta})$ and $p(R_{\tau+1:T}\mid
S_{t:T},y_{1:T},\boldsymbol{\theta})$, respectively. To sample from each of
the full conditional distributions in turn, a number of the practical
decisions remain the same as in the previous example of Section 5.1. We first
fix several of the practical decisions to those in Appendix B and use blocks
of size $\ell=4$, overlapping blocks by one state to reduce the correlation
between states at the block boundaries. We test the performance of the
algorithm for $N\in\\{10,20,50\\}$ grid cells.

There are two main differences in this implementation compared to the previous
example. The first accounts for the near-chaotic state processes and the large
range of the data ($1.59\times 10^{5}$ units). We apply the current state-
centred Approach 3 of Section 4, which requires fewer grid cells for good
mixing properties due to the lower span for the finite cells required to
achieve efficient (more local) moves. Within this approach, we permit large
variability around large values for the state process by adjusting the
variance of the Gaussian distributions used to set the boundaries of the grid
cells: we set the variance proportional to the current state at each time
point. We try factors of proportionality of $0.1,0.25,0.5,1$, with finite grid
cells between the $q$ and $1-q$ quantiles with $q=0.01,0.1,0.2$, resulting in
(average) spans for the finite cells ranging from $\mathcal{S}=52-400$. The
second difference is that we ensure the grid cells are over a discrete space
bounded at zero by rounding the quantiles used to determine the grid cells,
setting the lower boundary to zero if needed.

####  5.2.2 Results

We assessed the performance of the PMPMH algorithm using $10$ independent MCMC
simulations of $50,000$ iterations and present the performance metrics in
Table 3. Where the implementations of the PMPMH algorithm converged,
convergence was achieved within $4,000-18,000$ iterations, taking $11-27$
hours using $1.6$ GHz CPU. Fairly consistent and efficient performance (in
terms of ESS per second, Table 3) is achieved using $20-50$ grid cells with
the finite cells covering a region of the state-space $0.05-0.25\%$ of the
range of the data ($\mathcal{S}\in[82,400]$).

Max N | Average S |  |
---|---|---|---
| $\mathbf{\mathcal{S}}$ | R | ESS | ESS/s
10 | $\leq 76$ | $\leq 70$ | - | -
| 82-100 | 75-90 | 2500 | 0.05
| 100-200 | 90-180 | 3200 | 0.07
| $>200$ | $>180$ | - | -
20 | $\leq 76$ | $\leq 70$ | - | -
| 82-100 | 75-90 | 4400 | 0.09
| 100-200 | 90-180 | 5400 | 0.09
| 200-300 | 180-260 | 7000 | 0.11
| 300-400 | 260-300 | 6700 | 0.10
50 | $\leq 76$ | $\leq 70$ | - | -
| 82-100 | 75-90 | 4400 | 0.06
| 100-200 | 90-180 | 5700 | 0.06
| 200-300 | 180-260 | 7300 | 0.08
| 300-400 | 260-300 | 6900 | 0.07
Table 3: Average ESS and ESS/s across state and parameter samples for various
sets of tuning parameters under Approach 3 with blocks of size $\ell=4$. Since
the span of the finite cells varies across time points, we provide the average
span provided by the MCMC output. The ESS has been rounded to the nearest 100
to account for the variability within the chains.

Approach 3 mixes effectively by producing a reasonable HMM approximation in a
relatively small region around the current state and is also computationally
cheap to implement. The efficiency of the algorithm reduces at an upper bound
of $50$ grid cells, indicating that the improved mixing properties of the
algorithm are not justified by the extra cost when compared with using fewer
grid cells. Conversely, $10$ grid cells, although lower in computational cost,
do not produce a sufficiently accurate approximation to the posterior in the
region of the current state. We also note that, similarly to the previous case
study, implementations relying on _extremely_ small moves around the state
space (finite cells $<0.05\%$ of the range of the data) gave poor mixing and
convergence properties.

The PGAS sampler, when applied to the model using both joint updating of the
latent states and updating each latent state process in turn from their full
conditional distributions, did not converge in $50,000$ iterations with $1000$
particles, remaining in a range of the state space unrepresentative of the
posterior due to sample impoverishment and taking $160$ hours with $1.6$ GHz
CPU.

##  6 Discussion

We provide a novel approach for fitting general state-space models to observed
data. The idea utilises tractable HMM approximations to efficiently update the
unobserved latent states in an Metropolis-Hastings algorithm. We demonstrate
the generality of the proposed approach by its application to two problems,
including a near-chaotic problem that is generally challenging. The proposed
PMPMH algorithm is demonstrated to provide reliable posterior estimates within
reasonable computational time frames and compared to the PGAS algorithm often
applied to such problems. The flexibility of the PMPMH approach via the tuning
parameters, including the location, size, and number of grid cells, provides
an adaptable and efficient algorithm. Using simple methods for placing the
grid cells, such as equally-sized grid cells, can be efficient since they are
computationally cheap. However, such simple methods may require a large number
of grid cells to achieve good mixing if the range of the data is large. A
current state-centered approach for approximating the HMM is useful,
especially when the state space covers a large range, and using a data-
centered approach is useful provided the range of the data is small.

The PMPMH approach presents some interesting points for future research. One
way to reduce the overall computational cost of the algorithm is to reduce the
number of times the HMM approximation is calculated across iterations. This is
valid for approximation methods that do not depend on the current state or
model parameters. However, for deterministic grids that adapt to previous
samples (for example, Approaches 2 and 3), the resulting posterior estimate
may become biased since the chain is no longer Markovian (Haario et al.,,
1999). However, Haario et al., (1999) also show that the bias introduced is
negligible in some applications and that unbiased samples can be obtained by
introducing some pre-determined stopping criterion for updating the grid, that
is, specifying a number of iterations beyond which the grid is no longer
updated. The computational gains from updating the HMM approximation less
frequently could, however, be balanced with potentially reduced mixing
properties. For example, local-move grids using a low span for finite cells
are more sensitive to the frequency of the HMM approximation since they rely
on proposed moves being made in the region of the current state. However,
updating the grid less frequently may give an efficient approach when applied
to more global-move samplers.

One may also attempt to find efficient ways to adapt the HMM approximations.
With any such adaptations, there is a trade-off between the improved mixing
properties and the computational cost arising from the complexity of the
method. For example, in Step 1 of the algorithm, more complex numerical
integration methods could be used to approximate the SSM density in each grid
cell. Such methods include increasing the number of points in each grid cell
used for approximation calculations, simple linear functions joining
evaluations at midpoints, or other fast evaluations that include the gradient
in each grid cell. There may also be ways to improve the efficiency of the
within-cell proposal distributions in Step 2, including approaches that
approximate the posterior density in each grid cell where the conditional
posteriors are irregular or vary in an unsystematic way depending on the value
of the parameter. These approaches could be valuable to explore for more
complex SSMs, in which the improved mixing may justify the added computational
cost. Further, these approaches to improving mixing in each step may permit
the use of fewer grid cells, decreasing the overall computational cost, which
may also permit greater scalability of the HMM to higher dimensional latent
spaces. Although we demonstrated the scalability of the proposed algorithm to
higher dimensional spaces by updating each state dimension conditionally on
the remaining state dimensions, higher dimensional spaces are a particular
challenge if individual state dimensions are highly correlated, resulting in
poor mixing. This scalability of grid-based methods to higher dimensional
spaces is a challenge to the proposed algorithm and an active area of
research.

Supplementary material

GitHub repository containing the R code for the implementations of the PMPMH,
PG and PGAS algorithms in Section 5:
https://github.com/mallewellyn/PMPMH.

Acknowledgements

Funding: M.L. was supported by studentship funding from the Engineering and
Physical Sciences Research Council (EPSRC). R.K. was supported by the
Leverhulme Research Fellowship RF-2019-299. V.E. was supported by the
Leverhulme Research Fellowship RF-2021-593.

##  Appendix A Parameter prior distributions

###  A.1 Gaussian Mixture State-Space Model Parameters in Section 5.1

In Section 5.1, we conditionally update the parameters from both
$p(\boldsymbol{\theta}\mid x_{1:T},y_{1:T}^{(1)})$ and
$p(\boldsymbol{\theta}\mid x_{1:T},y_{1:T}^{(2)})$ using single-site conjugacy
in the observation variance parameter. For both Model 1 and Model 2, we assign
vague priors to enable comparison with the simulated parameter for diagnostic
purposes. For Model 1,

$\sigma_{\epsilon}^{2}\sim\Gamma^{-1}(2,2),$

resulting in a Gibbs sampler for $\sigma_{\epsilon}^{2}$. We assign further
high-variance priors to the other parameters in Model 1:

$\displaystyle p$ $\displaystyle\sim U(0,1),$
$\displaystyle\sigma_{\eta_{1}}^{2}$ $\displaystyle\sim\Gamma^{-1}(2,2),$
$\displaystyle\sigma_{\eta_{2}}^{2}$ $\displaystyle\sim\Gamma^{-1}(2,700),$

where $p$, $\sigma_{\eta_{1}}^{2}$ and $\sigma_{\eta_{2}}^{2}$ are proposed at
iteration $m$ from uniformly-distributed random walk proposals over intervals
of length $0.3$, $2$, and $160$ units respectively. For Model 2, we similarly
assign the independent priors:

$\displaystyle\sigma_{\epsilon}^{2}$ $\displaystyle\sim\Gamma^{-1}(2,10),$
$\displaystyle p$ $\displaystyle\sim U(0,1),$
$\displaystyle\sigma_{\eta_{1}}^{2}$ $\displaystyle\sim\Gamma^{-1}(2,2),$
$\displaystyle\sigma_{\eta_{2}}^{2}$ $\displaystyle\sim\Gamma^{-1}(2,700),$

where the prior for $\sigma_{\epsilon}^{2}$ again results in Gibbs steps for
this parameter. We propose $p$, $\sigma_{\eta_{1}}^{2}$ and
$\sigma_{\eta_{2}}^{2}$ at iteration $m$ from a uniform random walk proposal
distribution over $0.02$, $0.5$, and $20000$ units respectively. All intervals
for the random walk proposal distributions are set using pilot tuning.

###  A.2 Nicholson’s Blowfly Model Parameters in
Section 5.2

We give details for the priors and sampling schemes for the model parameters,
$\boldsymbol{\theta}$, in Section 5.2. Once again, we make use of the
conjugate priors for single site updates of $\boldsymbol{\theta}$ where
possible and assign vague priors to help diagnose convergence:

$\displaystyle\delta$ $\displaystyle\sim\Gamma(0.007,0.01),$ $\displaystyle P$
$\displaystyle\sim\Gamma(50,1),$ $\displaystyle\beta_{\epsilon}$
$\displaystyle\sim\Gamma^{-1}(100,100),$ $\displaystyle\beta_{e}$
$\displaystyle\sim\Gamma^{-1}(10,1),$ $\displaystyle\phi$
$\displaystyle\sim\Gamma(0.01,0.01).$

This results in single-site Gibbs updates for $P$ and $\phi$. For $\delta$,
$\beta_{\epsilon}$ and $\beta_{e}$ and we use a random walk Metropolis-
Hastings step with uniform proposal distributions over intervals of length
$0.03$, $0.5$, and $0.05$ respectively. The parameters of the proposal
distributions were chosen via pilot tuning.

##  Appendix B Practical decisions

We make a number of practical decisions in order to implement the PMPMH
algorithm in our case studies. We opt for simple choices for most of these and
discuss these here.

###  B.1 Deterministic integration method

In Section 3.1, we focus on the simple case for approximating the HMM
probabilities in each grid cell: midpoint integration using $S=1$. However,
the Riemann sum integral approximation method scales simply for higher order
polynomials where $S\geq 2$, or these methods can be easily replaced by more
complex numerical integration strategies. However, the complexity of these
methods should be balanced with their associated computational cost to ensure
efficiency in this step. We therefore test the efficiency of the algorithm
using midpoint integration ($S=1$)

The _midpoint_ and length of the cells must be defined to apply midpoint
integration, even when the $n=1$ or $n=N$ grid cells have an infinite range.
In these infinite cells, we simply set this artificial length at the average
length of the finite cells, and we set the artificial midpoint at half that
length away from the corresponding upper or lower boundary of the finite
cells. Further, to ensure that the HMM probabilities are sufficiently high to
avoid the proposal distribution resulting in a ‘near-reducible’ Markov chain,
we set a lower bound on the HMM probabilities. In all the implementations of
this paper, we use a lower bound of $0.01$ and renormalize so that they sum to
one.

###  B.2 Proposal distributions within the grid cells

Once we have sampled a set of grid cells, indexed by $b^{\prime}_{1:T}$ as per
Section 3.2, we sample point values from within the grid cells using simple
proposal distributions for each $t\geq 1$. In all implementations, we sample
from uniform distributions with domain in the finite grid cells.

To sample values for the state in infinite grid cells in the first case study
of Section 5.1, we sample from a truncated Gaussian distribution with mean
equal to the finite boundary and a variance of $5$. The variance of the
infinite-cell distribution is relatively low so that proposals in this grid
cell are mostly concentrated around the boundary of the finite cells,
increasing the density in the tails of the proposal distribution. In the
second case study, we sample from a Poisson distribution with mean parameter
equal to $2$ (again to ensure a relatively heavy-tailed proposal
distribution), shifted to the lower bound of the upper (infinite) grid cell.

##  Appendix C Pseudocode for the PMPMH algorithm with block updates

Here we present the pseudocode describing the two stages of the PMPMH
algorithm with the states updated in blocks. This supplements the description
of the block-updating procedure in Section 3.4.

Note that, in this case, the acceptance probability for a proposed set of
$\ell$ states starting at time $t$, $x^{\prime}_{t:t+\ell-1}$, is given by

$\displaystyle
p(x_{t:t+\ell-1},x^{\prime}_{t:t+\ell-1}\mid\boldsymbol{\theta})=$
$\displaystyle\quad\frac{p(x^{\prime}_{t:t+\ell-1}\mid
x_{t-1},x_{t+\ell},y_{t:t+\ell-1},\boldsymbol{\theta})}{p(x_{t:t+\ell-1}\mid
x_{t-1},x_{t+\ell},y_{t:t+\ell-1},\boldsymbol{\theta})}\times$
$\displaystyle\quad\quad\frac{q(x_{t:t+\ell-1}\mid
x_{t-1},x_{t+\ell},y_{t:t+\ell-1},\boldsymbol{\theta})}{q(x^{\prime}_{t:t+\ell-1}\mid
x_{t-1},x_{t+\ell},y_{t:t+\ell-1},\boldsymbol{\theta})}$ (C.1)

where in the first block, $p(x_{t:t+\ell-1}\mid
x_{t-1},x_{t+\ell},y_{t:t+\ell-1},\boldsymbol{\theta})=p(x_{t:t+\ell-1}\mid
x_{t+\ell},y_{t:t+\ell-1},\boldsymbol{\theta})$, and similarly for the
proposal density. Likewise, in the last block,
$p(x_{t:t+\ell-1}\mid
x_{t-1},x_{t+\ell},y_{t:t+\ell-1},\boldsymbol{\theta})=p(x_{t:t+\ell-1}\mid
x_{t-1},y_{t:t+\ell-1},\boldsymbol{\theta})$.

Algorithm 2 PMPMH algorithm: block updates

Input: initial values $x_{1:T}^{(0)}$, $\boldsymbol{\theta}^{(0)}$, a number
of iterations $M$. A blocking strategy for the states: a block size $\ell$ and
a set of starting points for each block
$\\{1,\dots,D\\}\subset\\{1,\dots,T\\}$. A Gibbs or Metropolis-Hastings
sampling scheme to update $\boldsymbol{\theta}$ from
$p(\boldsymbol{\theta}\mid x_{1:T},y_{1:T})$.

For $1\leq t\leq\ell$ (block $d=1$), define

$\hat{P}(B_{t}=k\mid
B_{d-1}=b_{d-1},y_{d:t},\boldsymbol{\theta})=\hat{P}(B_{t}=k\mid
y_{d:t},\boldsymbol{\theta}),$

$\hat{P}(B_{t}=k\mid
B_{d-1}=b_{d-1},B_{t+1}=b_{t+1},y_{d:d+\ell-1},\boldsymbol{\theta})=\hat{P}(B_{t}=k\mid
B_{t+1}=b_{t+1},y_{d:d+\ell-1},\boldsymbol{\theta}).$

Also, for $t=T$ in block $D$ define

$\hat{P}(B_{t}=k\mid
B_{d-1}=b_{d-1},B_{t+1}=b_{t+1},y_{d:t},\boldsymbol{\theta})=\hat{P}(B_{T}=k\mid
B_{d-1}=b_{d-1},y_{d:t},\boldsymbol{\theta}).$

for iterations $m=1,\dots M$ do

1:Update $\boldsymbol{\theta}^{(m)}$ from $p(\boldsymbol{\theta}\mid
x_{1:T}^{m-1},y_{1:T})$

2:

3:for $d=1,\dots,D$ do

4: Define intervals with indices $B_{\max(1,d-1):\min(d+\ell,T)}$

5: for $t=d,\dots d+\ell-1$ do $\triangleright$ (Step 1)

6: for $k=1,\dots N$ do

7: calculate $\hat{P}(B_{t}=k\mid
B_{d-1}=b_{d-1},y_{d:t},\boldsymbol{\theta}^{(m)})$

8:

9: sample $b^{\prime}_{d+\ell-1}$ from

10: $\\{k,\hat{P}(B_{d+\ell-1}=k\mid
B_{d+\ell}=b_{d+\ell},B_{d-1}=b_{d-1},y_{d:d+\ell-1},\boldsymbol{\theta}^{(m)})\\}_{k=1}^{N}$

11: for $t=d+\ell-2,\dots,d$ do

12: sample $b^{\prime}_{t}$ from

13: $\\{k,\hat{P}(B_{t}=k\mid
B_{t+1}=b^{\prime}_{t+1},B_{d-1}=b_{d-1},y_{d:t},\boldsymbol{\theta}^{(m)})\\}_{k=1}^{N}$

14:

15: for $t=d,\dots,d+\ell-1$ do$\triangleright$ (Step 2)

16: sample $x^{\prime}_{t}$ from $q(x_{t}\mid B_{t}=b^{\prime}_{t})$

17: sample $u\sim U(0,1)$

18: if
$u<p(x_{d:d+\ell-1}^{(m-1)},x^{\prime}_{d:d+\ell-1}\mid\boldsymbol{\theta})$
then $\triangleright$ Equation (C.1)

19: $x_{d:d+\ell-1}^{(m)}\leftarrow x^{\prime}_{d:d+\ell-1}$

20: else $x_{d:d+\ell-1}^{(m)}\leftarrow x^{(i-1)}_{d:d+\ell-1}$

## References

  * Andrieu et al., (2010) Andrieu, C., Doucet, A., and Holenstein, R. (2010).  Particle Markov chain Monte Carlo methods.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 72(3):269–342.
  * Andrieu and Roberts, (2009) Andrieu, C. and Roberts, G. O. (2009).  The Pseudo-Marginal Approach for Efficient Monte Carlo Computations.  Annals of Statistics, 37(2):697–725.
  * Auger-Méthé et al., (2016) Auger-Méthé, M., Field, C., Albertsen, C. M., Derocher, A. E., Lewis, M. A., Jonsen, I. D., and Mills Flemming, J. (2016).  State-space models’ dirty little secrets: even simple linear Gaussian models can have estimation problems.  Scientific Reports, 6(1):26677.
  * Beaumont, (2003) Beaumont, M. A. (2003).  Estimation of Population Growth or Decline in Genetically Monitored Populations.  Genetics, 164(3):1139–1160.
  * Borowska and King, (2023) Borowska, A. and King, R. (2023).  Semi-complete data augmentation for efficient state-space model fitting.  Journal of Computational and Graphical Statistics, in press.
  * Brooks and Gelman, (1998) Brooks, S. P. and Gelman, A. (1998).  General Methods for Monitoring Convergence of Iterative Simulations.  Journal of Computational and Graphical Statistics, 7(4):434–455.
  * Buckland et al., (2004) Buckland, S. T., Newman, K. B., Thomas, L., and Koesters, N. B. (2004).  State-space models for the dynamics of wild animal populations.  Ecological Modelling, 171(1):157–175.
  * Bucy and Senne, (1971) Bucy, R. S. and Senne, K. D. (1971).  Digital synthesis of non-linear filters.  Automatica, 7(3):287–298.
  * Cappé et al., (2005) Cappé, O., Moulines, E., and Rydén, T. (2005).  Inference in Hidden Markov Models.  Springer Verlag, New York.
  * Carter and Kohn, (1994) Carter, C. K. and Kohn, R. (1994).  On Gibbs Sampling for State Space Models.  Biometrika, 81(3):541–553.
  * Chopin and Singh, (2015) Chopin, N. and Singh, S. S. (2015).  On particle Gibbs sampling.  Bernoulli, 21(3):1855–1883.
  * de Valpine and Hastings, (2002) de Valpine, P. and Hastings, A. (2002).  Fitting population models incorporating process noise and observation error.  Ecological Monographs, 72(1):57–76.
  * Deligiannidis et al., (2018) Deligiannidis, G., Doucet, A., and Pitt, M. K. (2018).  The correlated pseudomarginal method.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 80(5):839–870.
  * Durbin and Koopman, (2012) Durbin, J. and Koopman, S. J. (2012).  Time Series Analysis by State Space Methods: Second Edition.  Oxford University Press.
  * Fearnhead, (2011) Fearnhead, P. (2011).  MCMC for State-Space Models.  In Handbook of Markov Chain Monte Carlo, pages 513–529. Chapman & Hall/CRC.
  * Fearnhead and Meligkotsidou, (2016) Fearnhead, P. and Meligkotsidou, L. (2016).  Augmentation schemes for particle MCMC.  Statistics and Computing, 26(6):1293–1306.
  * Finke et al., (2016) Finke, A., Doucet, A., and Johansen, A. (2016).  On embedded hidden Markov models and particle Markov chain Monte Carlo methods.  arXiv.  https://doi.org/10.48550/arXiv.1610.08962.
  * Frühwirth-Schnatter, (2004) Frühwirth-Schnatter, S. (2004).  Efficient Bayesian parameter estimation.  In State Space and Unobserved Component Models: Theory and Applications, pages 123–151. Cambridge University Press.
  * Gelman and Rubin, (1992) Gelman, A. and Rubin, D. B. (1992).  Inference from Iterative Simulation Using Multiple Sequences.  Statistical Science, 7(4):457–472.
  * Giordani and Kohn, (2010) Giordani, P. and Kohn, R. (2010).  Adaptive Independent Metropolis—Hastings by Fast Estimation of Mixtures of Normals.  Journal of Computational and Graphical Statistics, 19(2):243–259.
  * Haario et al., (1999) Haario, H., Saksman, E., and Tamminen, J. (1999).  Adaptive proposal distribution for random walk Metropolis algorithm.  Computational Statistics, 14(3):375–395.
  * Herliansyah et al., (2022) Herliansyah, R., King, R., and King, S. E. (2022).  Laplace approximations for individual heterogeneity capture-recapture models.  Journal of Agricultural, Biological, and Environmental Statistics, 22(1):401–418.
  * Hobert, (2011) Hobert, J. P. (2011).  The Data Augmentation Algorithm: Theory and Methodology.  In Handbook of Markov Chain Monte Carlo, pages 253–293. Chapman & Hall/CRC.
  * Julier and Uhlmann, (1997) Julier, S. J. and Uhlmann, J. K. (1997).  New extension of the Kalman filter to nonlinear systems.  In Signal Processing, Sensor Fusion, and Target Recognition VI, pages 182–193.
  * Kalman, (1960) Kalman, R. E. (1960).  A New Approach to Linear Filtering and Prediction Problems.  Journal of Basic Engineering, 82(1):35–45.
  * Kantas et al., (2015) Kantas, N., Doucet, A., Singh, S. S., Maciejowski, J., and Chopin, N. (2015).  On Particle Methods for Parameter Estimation in State-Space Models.  Statistical Science, 30(3):328–351.
  * King, (2011) King, R. (2011).  Statistical ecology.  In Handbook of Markov Chain Monte Carlo, pages 419–447. Chapman & Hall/CRC.
  * King, (2012) King, R. (2012).  A review of bayesian state-space modelling of capture-recapture-recovery data.  Interface Focus, 2(2):190–204.
  * Kitagawa, (1987) Kitagawa, G. (1987).  Non-Gaussian State-Space Modeling of Nonstationary Time Series.  Journal of the American Statistical Association, 82(400):1032–1041.
  * Koopman and Bos, (2004) Koopman, S. J. and Bos, C. S. (2004).  State Space Models with a Common Stochastic Variance.  Journal of Business & Economic Statistics, 22(3):346–357.
  * Koyama et al., (2010) Koyama, S., Pérez-Bolde, L., Shalizi, C., and Kass, R. (2010).  Approximate Methods for State-Space Models.  Journal of the American Statistical Association, 105(489):170–180.
  * Kristensen et al., (2016) Kristensen, K., Nielsen, A., Berg, C. W., Skaug, H., and Bell, B. M. (2016).  TMB: Automatic Differentiation and Laplace Approximation.  Journal of Statistical Software, 70(5):1–21.
  * Langrock and King, (2013) Langrock, R. and King, R. (2013).  Maximum likelihood estimation of mark-recapture-recovery models in the presence of continuous covariates.  The Annals of Applied Statistics, 7(3):1709–1732.
  * Langrock et al., (2012) Langrock, R., MacDonald, I. L., and Zucchini, W. (2012).  Some nonstandard stochastic volatility models and their estimation using structured hidden Markov models.  Journal of Empirical Finance, 19(1):147–161.
  * Li et al., (2012) Li, T., Sattar, T. P., and Sun, S. (2012).  Deterministic resampling: Unbiased sampling to avoid sample impoverishment in particle filters.  Signal Processing, 92(7):1637–1645.
  * Lin et al., (2019) Lin, A., Zhang, Y., Heng, J., Allsop, S. A., Tye, K. M., Jacob, P. E., and Ba, D. (2019).  Clustering Time Series with Nonlinear Dynamics: A Bayesian Non-Parametric and Particle-Based Approach.  In Proceedings of the Twenty-Second International Conference on Artificial Intelligence and Statistics, pages 2476–2484.
  * Lindsten et al., (2014) Lindsten, F., Jordan, M. I., and Schön, T. B. (2014).  Particle Gibbs with Ancestor Sampling.  Journal of Machine Learning Research, 15(63):2145–2184.
  * Lindsten et al., (2012) Lindsten, F., Schön, T. B., and Jordan, M. (2012).  On the use of backward simulation in the particle gibbs sampler.  In 2012 IEEE International Conference on Acoustics, Speech and Signal Processing, pages 3845–3848.
  * Meent et al., (2015) Meent, J.-W., Yang, H., Mansinghka, V., and Wood, F. (2015).  Particle Gibbs with Ancestor Sampling for Probabilistic Programs.  In Proceedings of the Eighteenth International Conference on Artificial Intelligence and Statistics, pages 986–994.
  * Murphy and Godsill, (2016) Murphy, J. and Godsill, S. J. (2016).  Blocked Particle Gibbs Schemes for High Dimensional Interacting Systems.  IEEE Journal of Selected Topics in Signal Processing, 10(2):328–342.
  * Neal, (2003) Neal, R. M. (2003).  Markov Chain Sampling for Non-linear State Space Models Using Embedded Hidden Markov Models.  arXiv.  https://doi.org/10.48550/arXiv.math/0305039.
  * Newman et al., (2014) Newman, K. B., Buckland, S. T., Morgan, B. J. T., King, R., Borchers, D. L., Cole, D. J., Besbeas, P., Gimenez, O., and Thomas, L. (2014).  Modelling Population Dynamics: Model Formulation, Fitting and Assessment using State-space Methods.  Springer.
  * Newman et al., (2022) Newman, K. B., King, R., Elvira, V., de Valpine, P., McCrea, R. S., and Morgan, B. J. T. (2022).  State-space models for ecological time series data: Practical model-fitting.  Methods in Ecology and Evolution, 14(1):26 – 42.
  * Nonejad, (2015) Nonejad, N. (2015).  Particle Gibbs with ancestor sampling for stochastic volatility models with: heavy tails, in mean effects, leverage, serial dependence and structural breaks.  Studies in Nonlinear Dynamics & Econometrics, 19(5):561–584.
  * Patterson et al., (2008) Patterson, T. A., Thomas, L., Wilcox, C., Ovaskainen, O., and Matthiopoulos, J. (2008).  State-space models of individual animal movement.  Trends in Ecology & Evolution, 23(2):87–94.
  * Rabiner, (1989) Rabiner, L. (1989).  A tutorial on hidden Markov models and selected applications in speech recognition.  Proceedings of the IEEE, 77(2):257–286.
  * Rainforth et al., (2016) Rainforth, T., Naesseth, C. A., Lindsten, F., Paige, B., van de Meent, J.-W., Doucet, A., and Wood, F. (2016).  Interacting Particle Markov Chain Monte Carlo.  In Proceedings of the 33rd International Conference on Machine Learning, pages 2616–2625.
  * Shephard and Pitt, (1997) Shephard, N. and Pitt, M. K. (1997).  Likelihood Analysis of Non-Gaussian Measurement Time Series.  Biometrika, 84(3):653–667.
  * Shestopaloff and Neal, (2013) Shestopaloff, A. Y. and Neal, R. M. (2013).  MCMC for non-linear state space models using ensembles of latent sequences.  arXiv.  https://doi.org/10.48550/arXiv.1305.0320.
  * Shestopaloff and Neal, (2018) Shestopaloff, A. Y. and Neal, R. M. (2018).  Sampling Latent States for High-Dimensional Non-Linear State Space Models with the Embedded HMM Method.  Bayesian Analysis, 13(3):797–822.
  * Smith and Brown, (2003) Smith, A. C. and Brown, E. N. (2003).  Estimating a state-space model from point process observations.  Neural Computation, 15(5):965–991.
  * Tanner and Wong, (1987) Tanner, M. A. and Wong, W. H. (1987).  The Calculation of Posterior Distributions by Data Augmentation.  Journal of the American Statistical Association, 82(398):528–540.
  * Thygesen et al., (2017) Thygesen, U. H., Albertsen, C. M., Berg, C. W., Kristensen, K., and Nielsen, A. (2017).  Validation of ecological state space models using the Laplace approximation.  Environmental and Ecological Statistics, 24(2):317–339.
  * Wan and Van Der Merwe, (2000) Wan, E. A. and Van Der Merwe, R. (2000).  The unscented kalman filter for nonlinear estimation.  In Proceedings of the IEEE 2000 Adaptive Systems for Signal Processing, Communications, and Control Symposium, pages 153–158.
  * Wang et al., (2017) Wang, X., Li, T., Sun, S., and Corchado Rodríguez, J. (2017).  A Survey of Recent Advances in Particle Filters and Remaining Challenges for Multitarget Tracking.  Sensors, 17(12):2707.
  * Wood, (2010) Wood, S. N. (2010).  Statistical inference for noisy nonlinear ecological dynamic systems.  Nature, 466(7310):1102–1104.