Analysis of a nonlinear importance sampling scheme for Bayesian parameter estimation in state-space models

TL;DR

This paper provides a rigorous convergence analysis of a nonlinear importance sampling scheme for Bayesian parameter estimation in state-space models, demonstrating optimal convergence rates even with approximate importance weights.

Contribution

It offers the first theoretical proof of convergence for the nonlinear population Monte Carlo method, including the optimal rate and the property of exact approximation.

Findings

Convergence of the NPMC scheme is almost sure with rate M^{-1/2}.

The scheme achieves optimal Monte Carlo convergence despite constant mean error in importance weights.

Simulation confirms theoretical convergence properties in a target tracking model.

Abstract

The Bayesian estimation of the unknown parameters of state-space (dynamical) systems has received considerable attention over the past decade, with a handful of powerful algorithms being introduced. In this paper we tackle the theoretical analysis of the recently proposed {\it nonlinear} population Monte Carlo (NPMC). This is an iterative importance sampling scheme whose key features, compared to conventional importance samplers, are (i) the approximate computation of the importance weights (IWs) assigned to the Monte Carlo samples and (ii) the nonlinear transformation of these IWs in order to prevent the degeneracy problem that flaws the performance of conventional importance samplers. The contribution of the present paper is a rigorous proof of convergence of the nonlinear IS (NIS) scheme as the number of Monte Carlo samples, $M$, increases. Our analysis reveals that the NIS approximation errors converge to 0 almost surely and with the optimal Monte Carlo rate of $M^{-\frac{1}{2}}$. Moreover, we prove that this is achieved even when the mean estimation error of the IWs remains constant, a property that has been termed {\it exact approximation} in the Markov chain Monte Carlo literature. We illustrate these theoretical results by means of a computer simulation example involving the estimation of the parameters of a state-space model typically used for target tracking.

Equations

$${\pi }_{n,\theta }(A)={\int }_A{\pi }_{n,\theta }(\mathsf{d}\mathbf{x})$$

$${\xi }_{n,\theta }(A)={\int }_A{\xi }_n(\mathsf{d}\mathbf{x})$$

$${\pi }_{0,\theta }^N(\mathsf{d}{\mathbf{x}}_0)=\frac{1}{N}\sum _{i=1}^N{\delta }_{{\mathbf{x}}_0^i}(\mathsf{d}{\mathbf{x}}_0),$$

$$u_n^i=\frac{{\tilde{u}}_n^i}{{\sum }_{j=1}^N{\tilde{u}}_n^j},i=1,...,N.$$

$$(f,\pi )=\int f(\mathbf{x})\pi (\mathsf{d}\mathbf{x})$$

$$\underset{N\to \infty }{\lim }(f,{\pi }_{n,\theta }^N)=(f,{\pi }_{n,\theta })$$

$$||(f,{\pi }_{n,\theta }^N)-(f,{\pi }_{n,\theta })|{|}_p\le \frac{C_n||f|{|}_{\infty }}{\sqrt{N}}$$

$$(f,{\pi }_{n,\theta }^N)=\int f({\mathbf{x}}_n){\pi }_{n,\theta }^N(\mathsf{d}{\mathbf{x}}_n)=\frac{1}{N}\sum _{i=1}^{N}f({\mathbf{x}}_n^i).$$

$${\xi }_{n,\theta }^N(\mathsf{d}{\mathbf{x}}_n)=\frac{1}{N}\sum _{i=1}^N{\delta }_{{\tilde{\mathbf{x}}}_n^i}(\mathsf{d}{\mathbf{x}}_n).$$

$$\ell (\mathbf{y}|\theta )=\prod _{k=1}^n(l_{k,\theta }({\mathbf{y}}_k|\cdot ),{\xi }_{k,\theta }),$$

$$(l_{k,\theta }({\mathbf{y}}_k|\cdot ),{\xi }_{k,\theta })={\int }_{\mathcal{X}}{l}_{k,\theta }({\mathbf{y}}_k|{\mathbf{x}}_k){\xi }_{k,\theta }(\mathsf{d}{\mathbf{x}}_k).$$

$${\ell }^N(\mathbf{y}|\theta )=\prod _{k=1}^n(l_{k,\theta }({\mathbf{y}}_k|\cdot ),{\xi }_{k,\theta }^N)$$

$$p(\theta |\mathbf{y})\propto \ell (\mathbf{y}|\theta )p_0(\theta )$$

$${\mu }_k=\sum _{i=1}^M{\mathsf{w}}_{k-1}^i{\theta }_{k-1}^i\text{mbox}and{\Sigma }_k=\sum _{i=1}^M{\mathsf{w}}_{k-1}^i\left({\theta }_{k-1}^i-{\mu }_k\right){\left({\theta }_{k-1}^i-{\mu }_k\right)}^{\top }.$$

$$\hat{\sf w}_{k}^{j}={\mathcal{T}}_{M}\left(j,\{\tilde{\sf w}_{k}^{l}\}_{l=1}^{M}\right)=\left\{\begin{array}[]{ll}\tilde{\sf w}_{k}^{i_{M_{c}}},&\mbox{if $\tilde{\sf w}_{k}^{j}\geq\tilde{\sf w}_{k}^{i_{M_{c}}}$},\\ \tilde{\sf w}_{k}^{j},&\mbox{if $\tilde{\sf w}_{k}^{j}<\tilde{\sf w}_{k}^{i_{M_{c}}}$},\\ \end{array}\right..$$

$${\hat{\theta }}_{\ast }={\int }_{\mathsf{S}}\theta {\mu }_{\mathbf{y}}(\mathsf{d}\theta ),$$

$$\text{mbox}MSE(\hat{\theta })={\int }_{\mathsf{S}}(\theta -\hat{\theta }{)}^2{\mu }_{\mathbf{y}}(\mathsf{d}\theta ).$$

$${\mu }_{\mathbf{y},k}^M(\mathsf{d}\theta )=\sum _{i=1}^M{\mathsf{w}}_k^i{\delta }_{{\theta }_k^i}(\mathsf{d}\theta ),$$

$$\text{mbox}MSE({\hat{\theta }}_k^M)=\sum _{i=1}^M{\mathsf{w}}_k^i\Vert {\theta }_k^i-{\hat{\theta }}_k^M{\Vert }^2.$$

$$\ell (\theta )\triangleq \ell (\mathbf{y}|\theta )\text{mbox}and{\ell }^N(\theta )\triangleq {\ell }^N(\mathbf{y}|\theta ).$$

$${\tilde{\mathsf{w}}}^i=g^N({\theta }^i)\triangleq \frac{{\ell }^N({\theta }^i)p_0({\theta }^i)}{q({\theta }^i)},$$

$${\hat{\mathsf{w}}}^i=[{\mathcal{T}}^M\circ g^N]({\theta }^i),$$

$$\Vert l{\Vert }_{\infty }=\underset{n\ge 1,{\mathbf{x}}_n\in \mathcal{X},\theta \in \mathsf{S}}{\sup }{l}_{n,\theta }({\mathbf{y}}_n|{\mathbf{x}}_n)<\infty .$$

$${\Vert \frac{p_0}{q}\Vert }_{\infty }=\underset{\theta \in \mathsf{S}}{\sup }\left|\frac{p_0(\theta )}{q(\theta )}\right|<\infty .$$

$$\max \{\ell (\theta ),{\ell }^N(\theta )\}\le \Vert l{\Vert }_{\infty }^R<\infty \text{mbox}andE\left[{\ell }^N(\theta )\right]=\ell (\theta )$$

$$(f,\mu )\triangleq {\int }_{\mathsf{S}}f(\theta )\mu (\mathsf{d}\theta ),$$

$$(f,\mu )\approx (f,{\mu }^M)=\sum _{i=1}^{M}f({\theta }^i){\mathsf{w}}^i,$$

$$|(f,{\mu }^M)-(f,\mu )|\le \frac{V_{f,ϵ}}{M^{\frac{1}{2}-ϵ}}.$$

$$(f,\mu )=\frac{(fg,q)}{(g,q)}$$

$$(f,{\mu }^M)=\frac{(f[{\mathcal{T}}^M\circ g^N],q^M)}{({\mathcal{T}}^M\circ g^N,q^M)}$$

Full text

\usetkzobj

all

1 Introduction

The estimation of the static unknown parameters of state-space dynamic models is a classical problem in statistical signal processing [1, 2, 3, 4, 5, 6] which has also received considerable attention, very recently, from the computational statistics community [7, 8, 9] (see also [10] for a recent survey) partly because of the ubiquity of the problem in science and engineering and partly because of the availability of more powerful computational resources to address it.

The particle Markov chain Monte Carlo (pMCMC) method originally proposed in [7] has been rapidly adopted by researchers in signal processing [11, 12, 6, 13, 14]. This is a Markov chain Monte Carlo (MCMC) algorithm [15] where the target probability density function (pdf) is the posterior density of the unknown parameters conditional on the available observations. This pdf is analytically intractable and, hence, it is approximated (for each element of the chain) via particle filtering [16, 17, 18, 19, 20]. The most popular MCMC schemes (including Metropolis and Metropolis-Hastings algorithms) admit a pMCMC implementation. A key feature of these methods is that they have the so-called exact approximation property. This means that, even if the acceptance test of the MCMC algorithm is only approximate (since the true target pdf is intractable), the stationary distribution of the Markov chain is still actual posterior density of the parameters. While popular, pMCMC procedures suffer from the same limitations as regular MCMC schemes [15, 21]:

Convergence of the chain is purely asymptotic (no convergence rates are known) and potentially very slow (a problem made worse by the particle approximation).

• 2.

The Monte Carlo samples in the chain are correlated, which reduces the accuracy of estimators compared to methods that produce independent samples.

• 3.

If the target pdf is multimodal, MCMC algorithms may get trapped in local maxima of the function.

An alternative to pMCMC methods is to employ schemes based on importance sampling (IS) [21]. This class of techniques includes population Monte Carlo (PMC) [22], the sequential Monte Carlo square (SMC2) of [23] or the nested particle filter of [9]. PMC is an iterative IS scheme in which the proposed functions used to generate Monte Carlo samples (and, hence, to approximate the posterior probability distribution of the unknown parameters) are improved across the iterations of the algorithm. See [24, 25, 26, 27] for recent applications, and new developments, of this methodology in statistical signal processing. SMC2 is a generalisation of the iterative batch importance sampling (IBIS) algorithm of [28]. It mimics the standard particle filter, but the Monte Carlo samples are drawn from the space of the (static) parameters and they are sequentially updated using a pMCMC kernel. All these methods, including SMC2, are batch, meaning that the whole record of observations is typically processed many times. A purely recursive version of the SMC2 algorithm has been proposed in [9]. The reduction in computational complexity, however, is obtained at the expense of a reduction in the convergence rate of the algorithm. It is worth mentioning that all these techniques (including pMCMC) can be fit within the theoretical framework of sequential Monte Carlo samplers introduced in [29].

The key feature of IS-based methods is that the Monte Carlo samples (used to approximate the target distribution) are generated from almost-arbitrary proposal functions and then assigned importance weights (IWs). While this is a very flexible approach, it suffers from the well-known problem of degeneracy of IWs [30, 18, 21, 8]: when the target pdf is concentrated in a very small region of the space of the unknowns, the largest IW tends to be orders of magnitude greater than all other IWs. As a result the IS-based scheme practically yields a degenerate one-sample approximation.

In this paper we address the analysis of the nonlinear population Monte Carlo (NPMC) algorithm proposed in [8]. In the latter scheme, the IWs undergo a nonlinear transformation to control their variance and, in this way, mitigate the degeneracy problem. In [8] it was proved that the approximation of the target distribution produced at each iteration of the NPMC method converges asymptotically, with the number of Monte Carlo samples $M$, and almost surely (a.s.). Therefore, the weight transformation preserves asymptotic convergence, while it has been shown through numerical examples that performance for finite $M$ is consistently improved compared to conventional PMC procedures. The analysis in [8], however

relies on the exact computation of the IWs, which is not feasible for general state-space models,

• 2.

and does not provide explicit convergence rates111Error rates are found in [8] for convergence in probability (not for almost sure convergence) when the IWs are computed exactly.

In this paper we analyse the performance of NPMC methods for the Bayesian estimation the unknown parameters of state space models. Based on some unbiasedness properties of particle filters, we prove that IS with nonlinearly-transformed IWs also yields asymptotic convergence when the weights are approximate, i.e., computed via a particle filter with a fixed computational budget that introduces non-vanishing errors. In other words, we prove that the nonlinear importance sampler enjoys the same exact approximation property as pMCMC and SMC2 algorithms. Moreover, the analysis of this paper also extends considerably the results of [8] by obtaining an explicit (and almost sure) estimation error rate of order $M^{-\frac{1}{2}+\epsilon}$, where $\epsilon>0$ is an arbitrarily small constant. This result holds for approximate weights and under mild assumptions typical of classic IS analyses. It is worth mentioning that the analytical approach developed in this paper can be applied, in a rather natural way, to the study of recently proposed PMC-like algorithms [25, 31] when the target distribution is the posterior density of the parameters of a state space model.

The rest of the paper is organised as follows. The necessary background material, including notation, state-space models and particle filters, is presented in Section 2. The nonlinear IS scheme and its iterative implementation (the NPMC algorithm) are detailed in Section 3 for the case in which the target probability distribution is the posterior distribution of the unknown parameters of a state-space model. In Section 4 we introduce the new analytical results on the convergence of nonlinear importance samplers, which is the main contribution of the paper. We illustrate the exact approximation property, and numerically compare the NPMC algorithm with a pMCMC scheme through computer simulations for a target tracking model in Section 5. Finally, some brief concluding remarks are made in Section 6.

2 Background and problem statement

2.1 State-space model

A Markov state-space model consists of two sequences of random variables (r.v.’s), $\{{\bf x}_{n}\}_{n\geq 0}$ and $\{{\bf y}_{n}\}_{n\geq 1}$. The first sequence, $\{{\bf x}_{n}\}$, is termed the system state. We assume it takes values on some space ${\mathcal{X}}\subseteq\mathbb{R}^{d_{x}}$, hence ${\bf x}_{n}$ is a random $d_{x}\times 1$ vector. The state dynamics are described by a prior probability measure ${\mathcal{K}}_{0}({\sf d}{\bf x}_{0})$ and a sequence of Markov kernels ${\mathcal{K}}_{n,\theta}({\sf d}{\bf x}_{n}|{\bf x}_{n-1})$ that depend on a parameter vector $\theta\in{\sf S}\subset\mathbb{R}^{d_{\theta}}$. In this paper, $\theta$ is assumed unknown and modelled as a random vector, with prior pdf $p_{0}(\theta)$ with respect to (w.r.t.) the Lebesgue measure. The support set of the parameter vector, ${\sf S}$, is assumed to be compact.

The state ${\bf x}_{n}$ cannot be observed directly. Instead, some noisy observations ${\bf y}_{n}\in\mathcal{Y}\subseteq\mathbb{R}^{d_{y}}$, $n=1,2,\ldots$, are collected. We note that ${\bf y}_{n}$ is a $d_{y}\times 1$ vector, with $d_{y}\neq d_{x}$ in general.

We assume that the observations are conditionally independent given the system states and the parameter vector $\theta$, with a conditional pdf w.r.t. the Lebesgue measure, denoted $l_{n,\theta}({\bf y}_{n}|{\bf x}_{n})>0$, which depends on the parameter vector $\theta$ as well.

2.2 The optimal filter and its Monte Carlo approximation

Let ${\bf y}_{1:n}=\{{\bf y}_{1},\ldots,{\bf y}_{n}\}$ denote the sequence of observations collected up the time $n$. The posterior probability measure of the state ${\bf x}_{n}$ conditional on the observations ${\bf y}_{1:n}$ and the parameter vector $\theta$ is denoted $\pi_{n,\theta}$, i.e., for any Borel set $A\subset{\mathcal{X}}$,

[TABLE]

is the posterior probability of the event “${\bf x}_{n}\in A$”, given $\theta$ and ${\bf y}_{1:n}$.

Similarly, $\xi_{n,\theta}$ denotes the posterior probability measure of ${\bf x}_{n}$ conditional on $\theta$ and ${\bf y}_{1:n-1}$ (i.e., not including ${\bf y}_{n}$). This is often referred to as the one-step-ahead predictive measure ([32], Chapter 10). For a Borel set $A\subset{\mathcal{X}}$,

[TABLE]

is the posterior probability of the event “${\bf x}_{n}\in A$”, given $\theta$ and ${\bf y}_{1:n-1}$.

We refer to $\pi_{n,\theta}$ as the optimal filter conditional on the parameter vector $\theta$. It is not possible, in general, to obtain either $\pi_{n,\theta}$ or $\xi_{n,\theta}$ in closed-form (with the notable exception of linear-Gaussian state space models, for which $\pi_{n,\theta}$ and $\xi_{n,\theta}$ are computed recursively and exactly using the Kalman flter [33]) and, therefore, numerical approximation algorithms are needed. One of the most popular schemes is the standard particle filter, also known as bootstrap filter (BF) [16, 34, 18].

The BF with $N$ particles (i.e., Monte Carlo samples on the state space ${\mathcal{X}}$) conditional on a given parameter vector $\theta$ can be briefly outlined as follows.

Initialisation. Draw $N$ samples ${\bf x}_{0}^{1},\ldots,{\bf x}_{0}^{N}$ from the prior distribution ${\mathcal{K}}({\sf d}{\bf x}_{0})$. The particle approximation of $\pi_{0,\theta}({\sf d}{\bf x}_{0})\equiv{\mathcal{K}}_{0}({\sf d}{\bf x}_{0})$ is

[TABLE]

where $\delta_{{\bf x}_{0}^{i}}$ denotes the Dirac delta measure centred at ${\bf x}_{0}^{i}\in{\mathcal{X}}$. 2. 2.

Recursive step. Given the approximation $\pi_{n-1,\theta}^{N}({\sf d}{\bf x}_{n-1})=\frac{1}{N}\sum_{i=1}^{N}\delta_{{\bf x}_{n-1}^{i}}({\sf d}{\bf x}_{n-1})$, take the following steps:

1. (a)

Randomly propagate each particle using the Markov kernel in the model, i.e., draw $\tilde{\bf x}_{n}^{i}$ from ${\mathcal{K}}_{n,\theta}({\sf d}{\bf x}_{n}|{\bf x}_{n-1}^{i})$, $i=1,...,N$. 2. (b)

Compute IWs, $\tilde{u}_{n}^{i}=l_{n,\theta}({\bf y}_{n}|\tilde{\bf x}_{n}^{i})$, for $i=1,...,N$, and 3. (c)

normalise them as

[TABLE] 4. (d)

Resample: draw $N$ times independently from the discrete distribution $\tilde{\pi}_{n,\theta}^{N}({\sf d}{\bf x}_{n})=\sum_{i=1}^{N}u_{n}^{i}\delta_{\tilde{\bf x}_{n}^{i}}({\sf d}{\bf x}_{n})$ and denote the resulting samples as $\{{\bf x}_{n}^{i}\}_{i=1}^{N}$. Construct the unweighted approximation $\pi_{n,\theta}^{N}({\sf d}{\bf x}_{n})=\frac{1}{N}\sum_{i=1}^{N}\delta_{{\bf x}_{n}^{i}}({\sf d}{\bf x}_{n})$.

The resampling step (d) above can be implemented in a number of different ways (see, e.g., [35, 32] or [20] for a brief survey of methods). Here, for simplicity, we have adopted a scheme which is often referred to as multinomial resampling [18, 35] but most asymptotic convergence results hold true for several other schemes as well [36, 32]. The measure-valued r.v. $\pi_{n,\theta}^{N}$ is an approximation of the optimal filter $\pi_{n,\theta}$ (conditional on $\theta$). Let us use the shorthand

[TABLE]

for the integral of a real function $f:{\mathbb{R}}^{d}\rightarrow{\mathbb{R}}$ w.r.t. a measure $\pi$. Under very mild assumptions it can be shown that

[TABLE]

almost surely (a.s.) for any bounded function $f:{\mathcal{X}}\rightarrow{\mathbb{R}}$ [36, 32]. Moreover, if we denote $||f||_{\infty}=\sup|f({\bf x})|$, $E[Z]$ indicates the expected value of a r.v. $Z$ and $||Z||_{p}=(E[|Z|^{p}])^{\frac{1}{p}}$ is its $L_{p}$ norm ($p\geq 1$), then it can be proved [37] that

[TABLE]

where $C_{n}$ is a constant independent of $N$ and

[TABLE]

The algorithm also produces a Monte Carlo approximation of the predictive measure $\xi_{n,\theta}$, namely

[TABLE]

If we write ${\bf y}={\bf y}_{1:n}$ for the complete sequence of observations up to time $n$, it turns out that the conditional pdf of ${\bf y}$ given the parameter vector $\theta$, denoted $\ell({\bf y}|\theta)$, can be written in terms of integrals w.r.t. to the predictive measures $\xi_{k,\theta}$, $k=1,\ldots,n$. To be specific,

[TABLE]

where

[TABLE]

The conditional pdf $\ell({\bf y}|\theta)$ is the likelihood of the parameter vector $\theta$ given the available data ${\bf y}$ and the BF yields the straightforward estimator

[TABLE]

which can be shown to be unbiased (i.e., $E[\ell^{N}({\bf y}|\theta)]=\ell({\bf y}|\theta)$) under very mild assumptions ([36], Theorem 7.4.2).

2.3 Problem statement

Let ${\bf y}_{=}\{{\bf y}_{1},\ldots,{\bf y}_{R}\}$ be the available data set, with $R<\infty$. Our goal is to approximate the probability measure associated to the posterior pdf of the parameter vector, $\theta$, given the data, ${\bf y}$. We denote this pdf as $p(\theta|{\bf y})$ and it is straightforward to show, using Bayes’ theorem, that

[TABLE]

where, we recall, $p_{0}(\theta)$ is the prior pdf of $\theta$.

In the next section, we describe an iterative importance sampling algorithm, originally introduced in [8], for the approximation of $p(\theta|{\bf y}){\sf d}\theta$.

3 Algorithm

The NPMC algorithm of [8] is an iterative importance sampling (IS) scheme that seeks to approximate a target probability distribution, in our case given by the posterior pdf $p(\theta|{\bf y})$, using weighted Monte Carlo samples. It generates a sequence of proposal pdf’s $q_{k}(\theta)$, $k=1,\ldots,K$, from which samples can be drawn and importance weights (IWs) can be computed. This sequence of proposals is expected to yield increasingly better approximations of the target as the algorithm converges. The key feature of the NPMC method, which departs from the classical PMC technique of [22], is to compute a set of transformed importance weights (TIWs) by applying a nonlinear function to the standard IWs. The aim of this transformation is to mitigate the well-known problem of the degeneracy of the IWs (common to many IS methods, see [18, 8]) by controlling the weight variability.

For the case of general state space models, an additional difficulty encountered when trying to estimate the unknown model parameters (denoted $\theta$ in our setup) is that the likelihood $\ell({\bf y}|\theta)$ is intractable. In the last few years, though, it has become a common approach to approximate this likelihood via particle filtering (PF) (see, e.g., [8, 7, 38, 23]). To be specific, we let $\ell^{N}({\bf y}|\theta)$ stand for the approximation of $\ell({\bf y}|\theta)$ computed using a standard bootstrap filter (BF) [16, 39] with $N$ particles (see equation (12) in Section 2.2). One key feature of this approach, that we exploit for our analysis in Section 4, is that $\ell^{N}({\bf y}|\theta)$ can be proved to be an unbiased estimator of $\ell({\bf y}|\theta)$ [36, 40].

The NPMC algorithm applied to a state space model, with $K$ iterations, $M$ Monte Carlo samples per iteration, plain Gaussian proposals $\{q_{k}\}_{k\geq 1}$, and approximate likelihoods is outlined below.

Initialisation. Draw $M$ i.i.d. samples $\theta_{0}^{1},\theta_{0}^{2},\ldots,\theta_{0}^{M}$ from the prior pdf $p_{0}(\theta)$. Then,

compute non-normalised IWs $\tilde{\sf w}_{0}^{i}\propto\ell^{N}({\bf y}|\theta_{0}^{i})$, $i=1,...,M$, 2. 2.

compute TIWs as $\hat{\sf w}_{0}^{i}={\mathcal{T}}_{M}\left(i,\{\tilde{\sf w}_{0}^{j}\}_{j=1}^{M}\right)$, where ${\mathcal{T}}_{M}:\{1,\ldots,M\}\times\{\tilde{\sf w}_{0}^{j}\}_{j=1}^{M}\rightarrow[0,+\infty)$ is a nonlinear transformation, and 3. 3.

normalise the TIWs, ${\sf w}_{0}^{i}=\frac{\hat{\sf w}_{0}^{i}}{\sum_{j=1}^{M}\hat{\sf w}_{0}^{j}}$, $i=1,...,M$.

Iteration. For $k=1,\ldots,K$, take the following steps:

Let $q_{k}(\theta)={\cal N}(\theta|\mu_{k},\Sigma_{k})$ be a multivariate Gaussian pdf with mean vector and covariance matrix obtained, respectively, as

[TABLE]

Note that the random variates $\theta_{k-1}^{i}$, $i=1,...,M$, are $d_{\theta}\times 1$ vectors. The superscript ⊤ denotes transposition. 2. 2.

Draw i.i.d. samples $\theta_{k}^{i}$, $i=1,...,M$, from $q_{k}(\theta)$. 3. 3.

Compute IWs, $\tilde{\sf w}_{k}^{i}=\frac{\ell^{N}({\bf y}|\theta_{k}^{i})p_{0}(\theta_{k}^{i})}{q_{k}(\theta_{k}^{i})}$, $i=1,...,M$. 4. 4.

Compute TIWs, $\hat{\sf w}_{k}^{i}={\mathcal{T}}_{M}\left(i,\{\tilde{\sf w}_{k}^{j}\}_{j=1}^{M}\right)$, $i=1,...,M$, using the same nonlinear map as for $k=0$. 5. 5.

Normalise the TIWs, ${\sf w}_{k}^{i}=\frac{\hat{\sf w}_{k}^{i}}{\sum_{j=1}^{M}\hat{\sf w}_{k}^{j}}$, $i=1,...,M$.

Following [8], the nonlinear map ${\mathcal{T}}_{M}$ of choice is a “clipping” transformation. In particular, let $i_{1},i_{2},...,i_{M}$ be a permutation of the indices $1,2,...,M$ such that the IWs become ordered, namely $\tilde{\sf w}_{k}^{i_{1}}\geq\tilde{\sf w}_{k}^{i_{2}}\geq\cdots\geq\tilde{\sf w}_{k}^{i_{M}}$. The clipping transformation ${\mathcal{T}}_{M}$, with parameter $1\leq M_{c}\leq\sqrt{M}$, flattens the $M_{c}$ largest IWs and makes them equal to the $M_{c}$-th non-normalised IW, $\tilde{\sf w}_{k}^{i_{M_{c}}}$. Specifically, for each $j=1,...,M$, we obtain

[TABLE]

Other choices of ${\mathcal{T}}_{M}$ are possible (e.g., tempering schemes) but clipping has been found particularly effective in practice [8]. The choice of Gaussian proposals (in step 1 of the Iteration) is made merely for simplicity. Other (more efficient) possibilities exist, but we stick to this formulation as it is sufficient for the purpose of this paper.

Given $A\subseteq{\sf S}$, being ${\sf S}$ the support set of the parameter vector $\theta$ described in Section 2, let $\mu_{\bf y}(A)=\int_{A}p(\theta|{\bf y}){\sf d}\theta$ denote the posterior probability measure (conditional on the observed data ${\bf y}$) associated to the parameter vector $\theta$. This measure yields the full probabilistic description of $\theta$ given the available observations. If $\mu_{\bf y}$ is available, then we can compute various types of estimators and assess the associated errors. For example, the posterior-mean estimator is

[TABLE]

and it minimises the mean square error (MSE). For an arbitrary estimator $\hat{\theta}$, the MSE can also be written as an integral w.r.t. $\mu_{\bf y}({\sf d}\theta)$, namely,

[TABLE]

The proposed NPMC algorithm yields a sequence of importance sampling (i.e., weighted Monte Carlo) approximations of $\mu_{\bf y}({\sf d}\theta)$. To be specific, at each iteration $k$ we obtain the random probability measure

[TABLE]

where $\delta_{\theta_{k}^{i}}$ denotes the Dirac delta measure centred at $\theta_{k}^{i}$. Using $\mu_{{\bf y},k}^{M}({\sf d}\theta)$ we can approximate any parameter estimator. For instance, $\hat{\theta}^{M}_{k}=\sum_{i=1}^{M}{\sf w}_{k}^{i}\theta_{k}^{i}$ is the approximation of the posterior mean $\hat{\theta}_{*}$. The corresponding minimum MSE can also be approximately computed as

[TABLE]

In the next section we analyse the convergence of the approximate measure $\mu_{{\bf y},k}^{M}$ as $M\rightarrow\infty$ in a single iteration (i.e., for a given $k$) when the number of particles $N$ used to approximate the likelihood via the BF (i.e., the estimate $\ell^{N}({\bf y}|\theta)$ of $\ell({\bf y}|\theta)$) is kept constant and finite.

4 Analysis

Consider a single iteration $k$ in the NPMC algorithm, with a fixed importance density $q_{k}\equiv q$. We refer to the random measure $\mu_{{\bf y},k}^{M}({\sf d}\theta)=\sum_{i=1}^{M}{\sf w}_{k}^{i}\delta_{\theta_{k}^{i}}({\sf d}\theta)$ computed via the TIWs ${\sf w}_{k}^{i}$, $i=1,..,M$, as a nonlinear importance sampling (NIS) approximation of $\mu_{\bf y}({\sf d}\theta)$. Our aim in this section is to assess whether $\mu_{{\bf y},k}^{M}({\sf d}\theta)$ converges towards the true measure $\mu_{\bf y}({\sf d}\theta)$ or not as $M\rightarrow\infty$. To do this, there are two issues that need to be handled and make the analysis more difficult compared to a conventional IS method (that relies on the standard IWs, rather than the TIWs). These issues are:

• (i)

the distortion in the Monte Carlo approximation due to the clipping of the weights, which introduces additional bias (compared to the use of standard IWs); and

• (ii)

the impossibility to compute the IWs, and hence the TIWs, exactly, since the likelihood $\ell({\bf y}|\theta)$ is intractable and we work with the particle approximation $\ell^{N}({\bf y}|\theta)$ instead.

In [8] it was proved that, when the IWs can be computed exactly, the NIS approximation converges almost surely (a.s.) towards the target probability measure as $M\rightarrow\infty$, which accounts for (i) above222The analysis of [8] does not provide an error rate, though. Such rate is explicitly derived in this paper. The problem of the approximate computation of the weights was partially addressed in [41], for a relatively simple case where the errors in the IWs where assumed deterministic and bounded. However, the estimation problem studied in [41] (parameter estimation for $\alpha$-stable distributions using iid data) did not involve any dynamics and the convergence analysis only showed an upper bound for the approximation errors that included a deterministic constant, namely a non-vanishing term proportional to the approximation error of the IWs.

Here, we show stronger analytical results that ensure the almost sure convergence of the NIS approximation when $M\rightarrow\infty$ and the likelihood function can only be estimated as $\ell^{N}({\bf y}|\theta)$, i.e., using a BF with a finite and fixed number of particles $N$. Under assumptions which are standard in the classical IS theory, we prove that integrals of the form $\int f(\theta)\mu_{{\bf y},k}^{M}({\sf d}\theta)$ converge towards $\int f(\theta)\mu_{{\bf y},k}({\sf d}\theta)$ a.s. as $M\rightarrow\infty$ and provide explicit error rates.

4.1 Notation

Since we focus our attention in the NIS scheme alone, i.e., a single iteration of the proposed algorithm, in the remaining of this section we drop the iteration index $k$. Hence, we assume a fixed importance density $q(\theta)$, from where $M$ independent Monte Carlo samples, $\theta^{1},\theta^{2},\ldots,\theta^{M}$, are drawn. Since the observations ${\bf y}$ are assumed arbitrary but fixed, we drop them from the likelihood notation and write

[TABLE]

Similarly, we simplify the notation for the posterior pdf and write $p(\theta)=p(\theta|{\bf y})$ and $\mu({\sf d}\theta)=\mu_{\bf y}({\sf d}\theta)$. Then, the non-normalised IWs are approximated as

[TABLE]

where we have introduced the weight function $g^{N}\triangleq\ell^{N}p_{0}/q$ as a shorthand. This weight function is a random approximation of the deterministic function $g=\ell p_{0}/q$. The support of $g$ is the same as the support of $q$, $\ell$ and $p_{0}$, denoted ${\sf S}\subseteq\mathbb{R}^{d_{\theta}}$. We assume that $g(\theta)>0$ for every $\theta\in{\sf S}$ as well (a standard assumption in classical IS). It is also apparent that $p\propto gq$, where $p$ is the posterior pdf, and the proportionality constant is independent of $\theta$.

The non-normalised TIWs computed via the clipping function (15) are denoted

[TABLE]

where $\circ$ represents function composition and we omit the index argument of (15) for conciseness (its value is clear from the notation in any case). The normalised TIWs are ${\sf w}^{i}=\frac{\hat{\sf w}^{i}}{\sum_{j=1}^{M}\hat{\sf w}^{j}}$, and they are used to compute the approximate measure $\mu^{M}({\sf d}\theta)=\sum_{i=1}^{M}\delta_{\theta^{i}}({\sf d}\theta){\sf w}^{i}$.

4.2 Assumptions and a preliminary result

Let the state sequence $\{{\bf x}_{n}\}_{n\geq 0}$ take values on ${\mathcal{X}}\subseteq\mathbb{R}^{d_{x}}$. We make the following classical assumptions on the conditional pdf of the observations ${\bf y}_{n}$, $n=1,2,\ldots,R$, the prior density of the parameters, $p_{0}(\theta)$, and the importance function $q(\theta)$.

Assumption 1

The observation sequence ${\bf y}_{1:R}$ is arbitrary but fixed. The functions $l_{n}({\bf y}_{n}|\cdot):{\mathcal{X}}\rightarrow(0,\infty)$, $n=1,2,...,R$, are uniformly bounded, i.e., there exists a finite and positive constant $\|l\|_{\infty}$ such that

[TABLE]

Assumption 2

The ratio of pdf’s $\frac{p_{0}(\theta)}{q(\theta)}$ is bounded on ${\sf S}$, i.e., there exists a positive and finite constant $\left\|\frac{p_{0}}{q}\right\|_{\infty}$ such that

[TABLE]

Remark 1

If the parameter support set ${\sf S}$ is compact, then A.1 and A. 2 hold naturally for most models of practical interest.

The following lemma plays a key role in the asymptotic convergence analysis of the approximation $\mu^{M}({\sf d}\theta)$. It states that $\ell^{N}(\theta)$ is an unbiased estimator of the likelihood $\ell(\theta)$ and enables us to show that the NIS scheme converges when $M\rightarrow\infty$, even if the number of particles $N$ in the approximation $\ell^{N}(\theta)$ remains finite and constant.

Lemma 1

If Assumption 1 holds then

[TABLE]

independently of $N$.

Proof. From the definition of $\ell(\theta)$ in Eq. (10) and its estimator $\ell^{N}(\theta)$ in Eq. (12), it is clear that both $\ell(\theta)\leq\|l\|_{\infty}^{R}$ and $\ell^{N}(\theta)\leq\|l\|_{\infty}^{R}$ when $R$ is the number of available observations. The fact that $\ell^{N}(\theta)$ is unbiased is a consequence of [36, Theorem 7.4.2] (see also [40, Lemma 2] for an alternative proof that does not rely on the Feynmann-Kac framework). \qed

4.3 Asymptotic convergence, error rates and exact approximation

In the sequel we look into the approximation of integrals of the form

[TABLE]

where $f$ is a bounded real function on the parameter space ${\sf S}$. We use $\|f\|_{\infty}\triangleq\sup_{\theta\in{\sf S}}|f(\theta)|<\infty$ to denote the supremum norm of a bounded function, while the set of bounded functions on ${\sf S}$ is denoted $B({\sf S})$. The approximations of interest are

[TABLE]

for any $f\in B({\sf S})$.

The following theorem yields an explicit upper bound for the (random) approximation error $|(f,\mu^{M})-(f,\mu)|$. The bound is proportional to $M^{-\frac{1}{2}+\epsilon}$ (for an arbitrarily small $\epsilon>0$) and, therefore, it vanishes as $M\rightarrow\infty$, independently of the number of particles $N$ used in the approximate likelihoods $\ell^{N}(\theta^{i})$.

Theorem 1

Assume that A.1 and A.2 hold, $M_{c}\leq\sqrt{M}$ and $\int_{\sf S}\ell(\theta)p_{0}(\theta){\sf d}\theta=(\ell,p_{0})>0$. Then, for every $\epsilon\in\left(0,\frac{1}{2}\right)$ (arbitrarily small) and every $f\in B({\sf S})$ there exists a positive and a.s. finite r.v. $V_{f,\epsilon}$, independent of $M$ and $M_{c}$, such that

[TABLE]

In particular, $\lim_{M\rightarrow\infty}|(f,\mu^{M})-(f,\mu)|=0$ a.s.

Proof. Recall the intractable weight function $g=\ell p_{0}/q$ and its random estimator $g^{N}=\ell^{N}p_{0}/q$. The integral of any $f\in B({\sf S})$ w.r.t. the posterior measure $\mu({\sf d}\theta)\propto\ell(\theta)p_{0}(\theta){\sf d}\theta$ can be written as

[TABLE]

by simply noting that $g(\theta)q(\theta)=\ell(\theta)p_{0}(\theta)$. Similarly, for the random measure $\mu^{M}({\sf d}\theta)$ we can write

[TABLE]

where $q^{M}({\sf d}\theta)=\frac{1}{M}\sum_{i=1}^{M}\delta_{\theta^{i}}({\sf d}\theta)$ is the Monte Carlo approximation of the proposal distribution (with pdf $q(\theta)$) and $\circ$ denotes composition of functions, hence $[{\mathcal{T}}^{M}\circ g^{N}](\theta^{i})={\mathcal{T}}^{M}(g^{N}(\theta^{i}))$ is the transformed weight associated to $\theta^{i}$.

Given equations (29) and (30) it is straightforward to show that

[TABLE]

Since $(f,\mu^{M})\leq\|f\|_{\infty}<\infty$ and $(g,q)=(\ell,p_{0})$, where $(\ell,p_{0})>0$ by assumption, Eq. (31) readily yields

[TABLE]

and, therefore, the problem of calculating bounds for $|(f,\mu^{M})-(f,\mu)|$ reduces to the problem of computing bounds for errors of the form

[TABLE]

for $b\in B({\sf S})$.

Choose any $b\in B({\sf S})$. A simple triangle inequality yields

[TABLE]

It is straightforward to obtain an upper bound for the first term on the right hand side of the inequality (34). Indeed, by construction of ${\mathcal{T}}^{M}$ (see Eq. (15)) we readily obtain

[TABLE]

where the inequality follows from the bound $g^{N}\leq\|l\|_{\infty}^{R}\left\|\frac{p_{0}}{q}\right\|_{\infty}$, which is a straightforward consequence of assumptions A.1 and A.2 and the definition of the estimate $\ell^{N}$ produced by the BF (see Eq. (12)).

Finding a suitable bound for the second term on the right hand side of the inequality (34) takes some more effort. Choose, again, any $b\in B({\sf S})$. A simple triangle inequality yields

[TABLE]

Since $q^{M}=\frac{1}{M}\sum_{i=1}^{M}\delta_{\theta^{i}}$, for the second term on the right hand side of (36) we can write

[TABLE]

where the r.v.’s

[TABLE]

are independent, with zero mean (recall the $\theta^{(i)}$’s are i.i.d. draws from $q$) and bounded, because $b$ is bounded and A.1 and A.2 imply that $g<\|l\|_{\infty}^{R}\times\left\|\frac{p_{0}}{q}\right\|_{\infty}<\infty$. Therefore, it is an exercise in combinatorics to show that

[TABLE]

where $\tilde{c}$ is a constant independent of $M$ and $q$. Combining (38) with (37) readily yields

[TABLE]

The inequality (39) implies that there exists an a.s. finite r.v. $\tilde{U}_{b,\epsilon}>0$ such that

[TABLE]

where $0<\epsilon<\frac{1}{2}$ is an arbitrarily small constant independent of $M$ (see [42, Lemma 4.1]).

If we expand the first term on the right hand side of (36) we arrive at

[TABLE]

where the r.v.’s $Z_{N}^{i}=\frac{b(\theta^{i})p_{0}(\theta^{i})}{q(\theta^{i})}\left(\ell^{N}(\theta^{i})-\ell(\theta^{i})\right)$, $i=1,2,...,M$, are independent (because the samples $\theta^{1},\ldots,\theta^{M}$ are independent) and zero mean, as a result of Lemma 1333Note that $E\left[Z_{N}^{i}|\theta^{i}\right]=\frac{b(\theta^{i})p_{0}(\theta^{i})}{q(\theta^{i})}E\left[\ell^{N}(\theta^{i})-\ell(\theta^{i})\right]=0$, because $\ell^{N}(\theta^{i})$ is an unbiased estimator of $\ell(\theta^{i})$, hence $E\left[Z_{N}^{i}\right]=E\left[E\left[Z_{N}^{i}|\theta^{i}\right]\right]=0$.. Since they are also bounded, namely $|Z_{N}^{i}|\leq\|b\|_{\infty}\|l\|_{\infty}^{R}\left\|\frac{p_{0}}{q}\right\|_{\infty}$ as a consequence of A.1 and A.2, it is again an exercise to show that (41) implies

[TABLE]

in the same manner as we obtained the inequality (38). Resorting again to [42, Lemma 4.1], from (42) we deduce that there exists an a.s. finite r.v. $\bar{U}_{b,\epsilon}>0$, independent of $M$, such that

[TABLE]

where $0<\epsilon<\frac{1}{2}$ is an arbitrarily small constant independent of $M$.

Taking together (36), (40) and (43) we arrive at

[TABLE]

where $U_{b,\epsilon}=\tilde{U}_{b,\epsilon}+\bar{U}_{b,\epsilon}\geq 0$ is an a.s. finite r.v. independent of $M$, and $\epsilon\in\left(0,\frac{1}{2}\right)$ can be chosen to be arbitrarily small.

Substituting the inequalities (35) and (44) back into the relation (34) we arrive at the bound

[TABLE]

where the second inequality follows from the assumption $M_{c}\leq\sqrt{M}$ and choosing $\tilde{V}_{b,\epsilon}=2\|l\|_{\infty}^{R}\left\|\frac{p_{0}}{q}\right\|_{\infty}\|b\|_{\infty}+U_{b,\epsilon}$. Since the r.v. $U_{b,\epsilon}$ is a a.s. finite, $\tilde{V}_{b,\epsilon}<\infty$ a.s. as well.

To conclude the proof, we substitute the inequality (45) twice into the relation (32). To be precise, we choose $b=f$ first and use (45) to obtain a bound for the first term on the right hand side of (32). Then, we choose $b=1$ and apply (45) again to find a bound for the second term on the right hand side of (32). As a result, we arrive at

[TABLE]

Since $(\ell,p_{0})>0$ by assumption of Theorem 1, taking

[TABLE]

leads to the desired result and concludes the proof. \qed

Theorem 1 is a general result regarding nonlinear importance sampling. It holds true for any problem involving the approximation of the posterior probability distribution of the unknown parameters of a state space model as long as Assumptions 1 and 2 hold. These assumptions, in turn, are very mild and amount to the classical assumptions in the analysis of standard IS algorithms.

Remark 2

We draw attention to the fact that the error $|(f,\mu^{M})-(f,\mu)|$ vanishes a.s. when $M\rightarrow\infty$ even if the number of particles $N$ in the BF remains fixed and, hence, $\ell^{N}$ does not converge to $\ell$. This property has been coined “exact approximation” in the MCMC literature (see [7]).

5 Computer simulations

5.1 State-space models

In order to illustrate the performance of the NPMC algorithm and the exact approximation property granted by Theorem 1 we have carried out computer simulations for the estimation of the unknown parameters in a problem consisting of the tracking of a target moving over a region monitored by a network of sensors.

5.1.1 Target dynamics

The target moves over a closed rectangular region ${\mathcal{R}}=[-20,+20]\times[-10,+10]$. When it hits the border of ${\mathcal{R}}$, the target bounces back in according to the law of reflection [43]. The state of the system at time $n$ is {\bf x}_{n}=\left[\begin{array}[]{c}{\bf r}_{n}\\ {\bf v}_{n}\\ \end{array}\right]\in\mathbb{R}^{4}, where ${\bf r}_{n}\in{\mathcal{R}}$ is the target position and ${\bf v}_{n}$ its velocity. At time $n=0$, we assume a uniform prior on ${\mathcal{R}}$ for the position and a zero-mean Gaussian distribution for the velocity. To be specific, the prior probability measure is defined as

[TABLE]

where ${\bf I}_{2}$ is the $2\times 2$ identity matrix, ${\mathcal{U}}({\mathcal{R}})$ is the uniform distribution on ${\mathcal{R}}$ and ${\mathcal{N}}({\bf m},{\bf C})$ denotes the Gaussian distribution with mean $\bf m$ and covariance matrix $\bf C$.

At time $n>0$, the state vector ${\bf x}_{n}$ evolves according to a linear-Gaussian equation if the target position remains within the bounded region ${\mathcal{R}}$ but it “reflects” back in when the target reaches a border of ${\mathcal{R}}$. Specifically, let

[TABLE]

where ${\bf u}_{n}\sim{\mathcal{N}}({\bf 0},{\bf C})$ is a Gaussian noise term with [math]-mean and covariance matrix

[TABLE]

$\kappa$ is a time-discretisation step (we assume $\kappa=1$ in our simulations), $\sigma_{u}^{2}$ is a velocity variance parameter, and $\sigma_{z}^{2}$ is a position variance parameter. The latter are assumed known and identical, $\sigma_{u}^{2}=\sigma_{z}^{2}=10^{-2}$. If $\tilde{\bf x}_{n}$ generated in this way is inside ${\mathcal{R}}$, $\tilde{\bf x}_{n}\in{\mathcal{R}}$, then ${\bf x}_{n}=\tilde{\bf x}_{n}$, otherwise ${\bf x}_{n}=f({\bf x}_{n-1})$, where $f$ is the reflection function detailed in A. Note that we do not provide an expression for the kernel ${\mathcal{K}}_{n}({\sf d}{\bf x}_{n}|{\bf x}_{n-1})$ but have just described how to draw samples from it instead. This is enough for the implementation of the bootstrap filter and the NPMC algorithm.

For illustration, Fig. 1 depicts the region ${\mathcal{R}}$ and a sample trajectory (i.e., a sequence of positions ${\bf r}_{0},{\bf r}_{1},\ldots$) which hits the borders of ${\mathcal{R}}$ and is reflected back in at four different times. In the figure, the starting target position is represented by a red diamond, the direction of motion is indicated by arrows and the blue squares represent the position of the sensors used to monitor the target motion.

5.1.2 Observations

There are $J$ sensors deployed in ${\mathcal{R}}$ and, at time $n$, each sensor collects a measurement of the power of the radio signal transmitted by the target. To be specific, the observation recorded by sensor $j$ at time $n$ has the form

[TABLE]

where $P_{t}$ is the power of the transmitted radio signal, ${\bf s}_{j}$ is the location of the $j$th sensor, $||{\bf r}_{n}-{\bf s}_{j}||$ is the distance at time $n$ between the target and the sensor, $\nu>0$ is the path loss exponent, $\rho$ is the sensitivity of the sensor, i.e., the minimum power it can measure (note that $y_{j,n}\rightarrow 10\log(\rho)+\epsilon_{j,n}$ when $||{\bf r}_{n}-{\bf s}_{j}||\rightarrow\infty$) and $\epsilon_{j,n}\sim{\mathcal{N}}(0,\sigma_{\epsilon}^{2})$ is a Gaussian term accounting for observational errors. We assume $\sigma_{\epsilon}^{2}=1$ is a known parameter.

At each time instant $n$, a vector of $J$ observations ${\bf y}_{n}=[y_{1,n},y_{2,n},\ldots,y_{J,n}]^{T}\in\mathbb{R}^{J}$ is collected. The target is observed over $m$ time instants, and hence the available dataset is ${\bf y}={\bf y}_{1:m}$. We set $m=50$ for our computer simulations.

5.1.3 Problem statement

Given the state space model described in Sections 5.1.1 and 5.1.2 above, we aim at estimating the unknown parameters $P_{t}$, $\nu$ and $\rho$. All other parameters (namely the discretisation period $\kappa$ and the relevant variances) are assumed known. For all computer simulations we have set ground truth values $P_{t}=0.8$, $\nu=3$ and $\rho=10^{-5}$ for the parameters to be estimated.

Since $P_{t}>0$ and $\rho>0$, we apply the NPMC algorithm (together with competing algorithms to be described below) to approximate the posterior probability measure $\mu_{{\bf y}}({\sf d}\theta)$ of the vector of unknowns $\theta=[\log P_{t},\nu,\log\rho]^{T}\in\mathbb{R}^{3}$. We assume prior distributions of the form $\log P_{t}\sim{\mathcal{N}}(-0.11,0.22)$, $\nu\sim{\mathcal{N}}(0,4)$ and $\log\rho\sim{\mathcal{N}}(-11.02,0.4)$. Note that, in natural units, the prior mean and variance of $P_{t}$ are $1$ and $0.25$, respectively, while for $\rho$ the prior mean and variance are $2\times 10^{-5}$ and $2\times 10^{-10}$.

The likelihood $\ell({\bf y}|\theta)$ for the model does not have a closed form and, therefore, it is estimated using a BF, for the state space model described in Sections 5.1.1 and 5.1.2, to yield the approximation $\ell^{N}({\bf y}|\theta)$ detailed in Section 2.2.

5.2 Competing methods

We have applied to this problem the NPMC method described in Section 3, a standard PMC procedure and a particle Metropolis-Hastings (pMH) algorithm. The PMC scheme we have used is identical to the NPMC algorithm of Section 3 except that TIWs are not computed, hence all approximations rely on the conventional IWs.

The pMH is a representative of the class of particle MCMC methods [7] that have become popular in the past two years. It generates a Markov chain on the space of the unknown parameter vector $\theta$ according to the following procedure:

Draw $\theta_{0}\sim p_{0}(\theta)$ from the prior distribution of the parameters 2. 2.

At the $r$-th iteration, and given the previous element $\theta_{r-1}$:

1. (a)

Draw a tentative new element $\tilde{\theta}_{r}\sim{\mathcal{N}}(\theta_{r-1},\frac{2}{10}\bf C)$, where both $\bf C=\text{diag}\left(\left[0.22,4,0.4\right]\right)$ and the scale factor $\frac{2}{10}$ have been empirically chosen to optimise the performance of the algorithm. 2. (b)

Compute the (approximate) likelihood $\ell^{N}({\bf y}|\tilde{\theta}_{r})$ and prior density $p_{0}(\tilde{\theta}_{r})$. The acceptance probability for $\tilde{\theta}_{r}$ is

[TABLE] 3. (c)

Draw $u_{r}\sim{\mathcal{U}}(0,1)$. If $u_{r}<\alpha_{r}$ then $\theta_{r}=\tilde{\theta}_{r}$, else $\theta_{r}=\theta_{r-1}$.

When we generate a chain of length $L$ using the procedure above we set a burn-in period of $L\over 2$, hence estimates are computed from the samples $\theta_{\lfloor{L\over 2}\rfloor+1},\ldots,\theta_{L}$ in the chain.

To compare the pMH and PMC-like algorithms on a fair basis, we let $L=M\times K$, where $K$ is the number of iterations of the NPMC and PMC algorithms and $M$ is the number of samples generated per iteration.

All three methods (PMC, NPMC, pMH) rely on a BF with $N$ particles for the computation of $\ell^{N}({\bf y}|\theta)$. The value of $N$ is fixed for all algorithms as $N=400$ unless explicitly stated otherwise.

5.3 Results

Figure 2 shows the evolution of the MSE of the estimators of $\theta$ produced by the PMC, NPMC and pMH algorithms as the number of samples is increased.

The error for the NPMC algorithm is at least one order of magnitude below the errors of the conventional PMC and the pMH algorithms for every tested value of $M$. For $M=200$ samples, for example, the MSE attained by the NPMC is $\approx 1.19\times 10^{-2}$, while for the standard PMC and pMH algorithms the errors are $\approx 2.49\times 10^{-1}$ and $\approx 5.01$, respectively.

Next, we aim at finding out the length of the chain, $L$, required for the pMH algorithm to attain the same performance, in terms of MSE, as the NPMC algorithm. Figure 3 shows the MSE of the pMH method for different chain lengths (equivalently, number of generated samples).

For comparison, the performance of the NPMC algorithm for $M=500$ samples and $K=10$ iterations ($500\times 10=5,000$ Monte Carlo samples overall) is also indicated in the plot. It can be seen that, in the pMH algorithm, chains that are around $500,000$ samples long are required to attain the same MSE as the NPMC algorithm (a 100-fold increase of the computational cost). While the parameters of the pMH scheme may be further tuned to improve this performance, the gap between the algorithms is large enough to conclude that the NMPC method is more efficient in this example.

Finally, we examine the exact approximation property of the NPMC scheme stated by Theorem 1. Figure 4 shows the MSE of the NPMC algorithm versus the number of Monte Carlo samples, $M$, for different values of $N$ (the number of particles used by the BF to approximate the IWs). While Theorem 1 guarantees that the approximation errors vanish as $M\rightarrow\infty$, even if $N$ is fixed, it is reasonable to expect that for a fixed $M<\infty$, greater values of $N$ lead to better performance. This is shown, indeed, by Fig. 4. Note, however, that the difference in performance is very small. For $M=1,000$, the gap between the MSE of the NPMC scheme with $N=400$ and the NPMC scheme with $N=50$ is $\approx 6\times 10^{-3}$.

6 Conclusion

We have rigorously proved, under mild assumptions, that nonlinear importance samplers with clipped IWs converge a.s. with optimal Monte Carlo error rates even when the weights can only be estimated (and have a positive, non-vanishing variance) as long as these estimates are unbiased. Therefore, nonlinear importance samplers can perform exact approximation in the same manner as, e.g., particle MCMC schemes. Besides the theoretical contribution, we have numerically shown that the proposed algorithm can be more efficient than a particle Metropolis-Hastings algorithm of the same complexity for inference on a target tracking model.

Acknowledgments

This research has been partially supported by the Spanish Ministry of Economy and Competitiveness (projects TEC2015- 69868-C2-1-R ADVENTURE and FIS2013-40653-P), the Spanish Ministry of Education, Culture and Sport (mobility award PRX15/00378) and the Office of Naval Research (ONR) Global (Grant Award no. N62909-15-1-2011).

Appendix A Definition of function $f(\cdot)$

Let us denote the upper right, upper left, lower left and lower right vertices of the monitored region by, respectively, ${\bf c}_{0}$, ${\bf c}_{1}$, ${\bf c}_{2}$ and ${\bf c}_{3}$. The sides of the rectangle, obtained by joining adjacent vertices, are denoted ${\bf l}_{0}=\overline{{\bf c}_{1}{\bf c}_{0}}$ (top), ${\bf l}_{1}=\overline{{\bf c}_{1}{\bf c}_{2}}$ (left), ${\bf l}_{2}=\overline{{\bf c}_{2}{\bf c}_{3}}$ (bottom) and ${\bf l}_{3}=\overline{{\bf c}_{3}{\bf c}_{0}}$ (right). With this notation, Algorithm 1 can be used at time $n$ to generate a sample ${\bf x}_{n}=[{\bf r}_{n}^{\top},{\bf v}_{n}^{\top}]^{\top}$ from ${\bf x}_{n-1}=[{\bf r}_{n-1}^{\top},{\bf v}_{n-1}^{\top}]^{\top}$. It accounts for the scenario in which the target hits one of the walls and deals with it by means of the law of reflection [43].

We are implicitly assuming that ${\bf r}_{n}\in{\mathcal{R}}$ in step 5 above. If this is not the case, i.e., $r_{n}\notin{\mathcal{R}}$, then steps 3–5 can be run again to implement a second reflection.