scispace - formally typeset
Search or ask a question
Proceedings ArticleDOI

Particle filtering for Bayesian parameter estimation in a high dimensional state space model

TL;DR: A novel nested particle filtering scheme for joint Bayesian parameter estimation and tracking of the dynamic variables in a high dimensional state space model-namely a stochastic version of the two-scale Lorenz 96 chaotic system, commonly used as a benchmark model in meteorology and climate science.
Abstract: Researchers in some of the most active fields of science, including, e.g., geophysics or systems biology, have to deal with very-large-scale stochastic dynamic models of real world phenomena for which conventional prediction and estimation methods are not well suited. In this paper, we investigate the application of a novel nested particle filtering scheme for joint Bayesian parameter estimation and tracking of the dynamic variables in a high dimensional state space model-namely a stochastic version of the two-scale Lorenz 96 chaotic system, commonly used as a benchmark model in meteorology and climate science. We provide theoretical guarantees on the algorithm performance, including uniform convergence rates for the approximation of posterior probability density functions of the fixed model parameters.

Summary (1 min read)

I. INTRODUCTION

  • Researchers in some of the most active fields of science have to deal with very large scale stochastic dynamic models of real world phenomena for which conventional prediction and estimation methods are not well suited.
  • The authors investigate a nested particle filtering (PF) scheme for the Bayesian estimation of the dynamic variables and the static parameters of state space models.
  • The latter displays the basic physical features of atmospheric dynamics and, for this reason, the deterministic version of this model is commonly used as a benchmark for data assimilation [5] and parameter estimation techniques [6] in meteorology and climate science.
  • The authors illustrate the performance of the proposed scheme by means of computer simulations on a stochastic two-scale Lorenz 96 model [6] with 16 slow and 160 fast dynamic variables as well as several unknown parameters.
  • Furthermore, the authors obtain explicit convergence rates that link the computational cost of the PF algorithm, the kernel bandwidth and the dimension of the parameter space.

B. Proposed algorithm

  • The authors have carried out 20 independent simulations with N = 200 particles in the outer filter and M = 600 particles in the inner filters used to compute the importance weights of the outer filter.
  • For reference, the authors approximate them in each simulation run using the least squares estimator âLS = arg min.
  • The authors observed how the actual value is closely tracked.
  • This is the case for the remaining dynamic variables, X 2:J,t (not shown).
  • The mean square error of the dynamic variable estimates over the 20 independent simulations (normalised w.r.t. the power of the signals) was ⇡ 0.0313.

VI. CONCLUSIONS

  • The authors have proposed a nested PF scheme to jointly track the dynamic variables and approximate the posterior pdf of the fixed parameters in state space models.
  • The authors have proved a.s. convergence of the pdf approximations and displayed numerical results for a stochastic two-scale Lorenz 96 system.

Did you find this useful? Give us your feedback

Content maybe subject to copyright    Report

Particle filtering for Bayesian parameter estimation
in a high dimensional state space model
Joaqu
´
ın M
´
ıguez, Dan Crisan and In
´
es P. Mari
˜
no
Abstract
Researchers in some of the most active fields of science, including, e.g., geophysics or systems biology, have to
deal with very-large-scale stochastic dynamic models of real world phenomena for which conventional prediction and
estimation methods are not well suited. In this paper, we investigate the application of a novel nested particle filtering
scheme for joint Bayesian parameter estimation and tracking of the dynamic variables in a high dimensional state
space model –namely a stochastic version of the two-scale Lorenz 96 chaotic system, commonly used as a benchmark
model in meteorology and climate science. We provide theoretical guarantees on the algorithm performance, including
uniform convergence rates for the approximation of posterior probability density functions of the fixed model
parameters.
Index Terms
Particle filtering, data assimilation, Bayesian parameter estimation, convergence analysis, kernel density
estimation.
J. M. is with Department of Signal Theory & Communications, Universidad Carlos III de Madrid. E-mail: joaquin.miguez@uc3m.es. Avda.
de la Universidad 30, 28911 Legan
´
es, Madrid (Spain).
D. C. is with Department of Mathematics, Imperial College London. E-mail: d.crisan@imperial.ac.uk. 180 Queens Gate, London SW7 2BZ
(UK).
I. P. M. is with Department of Biology, Geology, Physics & Inorganic Chemistry, Universidad Rey Juan Carlos. E-mail: ines.perez@urjc.es.
Tulip
´
an s/n, 28933 M
´
ostoles, Madrid (Spain).
J. M. acknowledges the financial support of the Spanish Ministry of Economy and Competitiveness (project COMPREHENSION TEC2012-
38883-C02-01) and the regional government of Madrid (program CASI-CAM-CM S2013/ICE-2845).
I. P. M. acknowledges the financial support of the Spanish Ministry of Economy and Competitiveness (projects TEC2012-38883-C02-01
and FIS2013-40653-P).

2
Particle filtering for Bayesian parameter estimation
in a high dimensional state space model
I. INTRODUCTION
Researchers in some of the most active fields of science have to deal with very large scale stochastic dynamic
models of real world phenomena for which conventional prediction and estimation methods are not well suited.
In fact, state-of-the-art methods for Bayesian model inference in computational statistics, e.g., the popular particle
Markov chain Monte Carlo (pMCMC) [1] and sequential Monte Carlo square (SMC
2
) [2] algorithms, are batch
techniques and, therefore, they are not well suited to the processing of long observation sequences. Although
some recursive algorithms exist [3], they only yield maximum likelihood (point) estimates for the parameters, and
hence they are subject to various convergence (and complexity) issues when the likelihood is multimodal, contains
singularities or cannot be computed exactly.
In this paper, we investigate a nested particle filtering (PF) scheme for the Bayesian estimation of the dynamic
variables and the static parameters of state space models. The method is similar to the SMC
2
algorithm, but enjoys
a purely recursive structure that makes it better suited for online estimation and dynamical systems (see [4] for
additional details). We apply the new scheme to a stochastic version of the (chaotic) Lorenz 96 system. The latter
displays the basic physical features of atmospheric dynamics and, for this reason, the deterministic version of this
model is commonly used as a benchmark for data assimilation [5] and parameter estimation techniques [6] in
meteorology and climate science. We illustrate the performance of the proposed scheme by means of computer
simulations on a stochastic two-scale Lorenz 96 model [6] with 16 slow and 160 fast dynamic variables as well as
several unknown parameters.
Besides the numerical results, we establish theoretical guarantees on the performance of the proposed scheme.
In particular, we prove that kernel approximations of posterior probability density functions (pdfs) of the model
parameters converge uniformly over the parameter space. Furthermore, we obtain explicit convergence rates that
link the computational cost of the PF algorithm, the kernel bandwidth and the dimension of the parameter space.
The rest of the paper is organised as follows. In Section II we state the Bayesian inference problem to be solved.
The proposed nested particle filtering scheme is outlined in Section III. Kernel density estimators, based on the
output of the nested PF algorithm, are investigated in Section IV. Section V contains computer simulation results
and, finally, Section VI is devoted to the conclusions.
II. PRO BL EM S TATE ME NT
A. Notation
We use upper-case letters to denote random variables (r.v.s), e.g., X, and lower-case letters to indicate their
realisations, e.g., X = x. Given a r.v. X with pdf (x), the integral of a real function f (x) with respect to (w.r.t.)
the probability measure (x)dx is denoted (f, ) ,
R
f(x)(x)dx. We use r.v. independently of the dimension of
X, i.e., X can be a random vector. The set of bounded real functions on the set X is denoted B( X). If f 2 B(X),
then kf k
1
, sup
x2X
|f(x)| < 1.
B. State space models
A state space model can be described in terms of two random sequences, {X
t
}
t0
and {Y
t
}
t1
, and a r.v., ,
taking values in the sets XR
d
x
, YR
d
y
and S R
d
, where d
x
, d
y
and d
are the corresponding dimensions
of X, Y and S, respectively. We refer to the sequence {X
t
}
t0
as the state (or signal) process and we assume that
it is a homogeneous Markov chain, with a priori pdf
0
(x) and Markov transition kernel
(x
t
|x
t1
)dx
t
, indexed
by a realisation =.
The sequence {Y
t
}
t1
is the observation process. Each r.v. Y
t
is assumed to be conditionally independent of all
other observations given X
t
and , namely
P
t
{Y
t
2 A|x
0:t
,,{y
k
}
k6=t
} = P
t
{Y
t
2 A|x
t
,}

3
for any Borel set A, where P
t
denotes the probability measure for the triple ({X
n
}
nt
, {Y
n
}
nt
, ). We assume
that the conditional pdf of Y
t
given X
t
, denoted g
(y
t
|x
t
) 0, can be evaluated up to a (possibly unknown)
proportionality constant. Given Y
t
= y
t
, the function g
y
t
(x
t
) , g
(y
t
|x
t
) on the state space is the likelihood of X
t
.
The a priori pdf of the parameter r.v. is denoted µ
0
() (note that X
0
and are a priori independent). The
set
µ
0
,
0
,
, {g
y
t
}
t1
describes a stochastic Markov state-space model in discrete time. Note that the model is
indexed by the d
-dimensional random parameter 2 S.
C. Stochastic Lorenz 96 model
The two-scale Lorenz 96 model is a deterministic system of nonlinear differential equations that displays chaotic
dynamics; see, e.g., [6]. The system dimension, i.e., the number of dynamic variables, can be scaled arbitrarily. A
stochastic version of the model can be easily obtained by converting each differential equation into a stochastic
differential equation driven by an independent and additive Wiener process. Here, for conciseness, we describe the
difference equations that result from the application of the Euler-Maruyama integration method to a model with J
slow variables, Z
j
, j =0,...,J 1, and L fast variables per slow variable,
¯
Z
l
, l =0,...,JL 1. We then obtain
Z
j,n
= Z
j,n1
Z
j1,n1
(Z
j2,n1
Z
j+1,n1
)
+
2
4
f Z
j,n1
hc
b
Lj1
X
l=(j1)L
¯
Z
l,n1
3
5
+
p
U
j,n
, (1)
¯
Z
l,n
=
¯
Z
l,n1
cb
¯
Z
l+1,n1
(
¯
Z
l+2,n1
¯
Z
l1,n1
)
+
cf
b
c
¯
Z
l,n1
+
hc
b
Z
b
l1
L
c,n1
¯
U
l,n
,
where > 0 denotes the discretisation period, n represents discrete time, f is a forcing parameter that controls
the turbulence of the chaotic flow, c determines the time scale of the fast variables, h controls the strength of the
coupling between the fast and slow variables, b determines the amplitude of the fast variables, {U
j,n
,
¯
U
l,n
}
l,j,n0
are sequences of independent and identically distributed (i.i.d.) standard Gaussian r.v.s, and , ¯>0 are scale
parameters.
We assume that observations can only be collected from this system once every n
o
discrete time steps. Therefore,
the observation process has the form
Y
t
=[Z
1,tn
o
,Z
2,tn
o
,...,Z
J,tn
o
]
>
+ V
t
, (2)
where t =1, 2,... and {V
t
}
t1
is a sequence of i.i.d. r.v.s with common pdf N(v
t
;0,
2
y
I
J
), which denotes a
J-dimensional Gaussian density with 0 mean and covariance matrix
2
y
I
J
, where I
J
is the J J identity matrix.
In computer experiments, system (1) is often employed to generate both ground-truth values for the slow variables
{Z
j,n
}
j,n0
and synthetic observations, {Y
t
}
t1
. As a forecast model for the slow variables it is common to use
the difference equation [6]
Z
j,n
= Z
j,n1
Z
j1,n1
(Z
j2,n1
Z
j+1,n1
)
+[f Z
j,n1
`(Z
j,n1
, a)] +
p
U
j,n
, (3)
where a =[a
1
, a
2
]
>
2 R
2
is a (constant) parameter vector and function `(Z
j,n1
, a) 2 R is an ansatz for the
coupling term
hc
b
P
Lj1
l=(j1)L
¯
Z
l,n1
in (1), to be specified later.
Equations (3) and (2) describe a state space model that can be expressed in terms of the general notation in
Section II-B. The state process at time n is
˜
X
n
=[Z
0,n
,...,Z
J1 ,n
]
>
and the transition pdf from time n 1 to
time n is
˜
x
n
|˜x
n1
)=Nx
n
; (˜x
n1
,),
2
x
I
J
), (4)
where =[f, a
>
]
>
2 R
3
,
2
x
=
2
and (˜x
n1
,) 2 R
J
is the deterministic transformation that accounts
for all the terms on the right hand side of (3) except the noise contribution
p
U
j,n
. Since we only collect
observations every n
o
continuous-time units, we need to put the dynamics of the states on the same time scale

4
(x
t
|x
t1
)=
Z
...
Z
˜
(x
t
|˜x
tn
o
1
)
n
o
2
Y
i=1
˜
x
tn
o
i
|˜x
tn
o
i1
x
(t1)n
o
+1
|x
t1
)
n
o
1
Y
j=1
d˜x
tn
o
j
. (5)
as the observation process {Y
t
}
t1
in Eq. (2). If we define X
t
=
˜
X
tn
o
2X= R
J
then the transition density from
X
t1
to X
t
follows readily from (4), as shown in Eq. (5) at the top of this page. While
(x
t
|x
t1
) cannot be
evaluated in closed form, it is straightforward to draw a sample from X
t
|x
t1
by simply running Eq. (3) n
o
times,
with starting point x
t1
. The likelihood function is
g
y
t
(x
t
) / exp
(
1
2
2
y
J
X
r=1
(y
r,t
x
i
r
,t
)
2
)
,
which follows from (2) and is independent of for this particular model.
D. Problem statement
Assume that the parameters in (3) are random, with prior pdf µ
0
. Then, the goal is to approximate sequence
of conditional pdfs of the parameters =[F, A
1
, A
2
]
>
given the available observations at each time step t.
We write µ
t
() to denote the pdf of conditional on Y
1:t
= y
1:t
. This pdf can be recursively decomposed as
µ
t
() /
t,
(y
t
)µ
t1
(), where
t,
(y
t
) is the pdf of the r.v. Y
t
conditional on Y
1:t1
= y
1:t1
and =. The
latter density, in turn, can be written as an integral, namely
t,
(y
t
)=(g
y
t
,
t,
), where
t
(x
t,
) is the predictive
pdf of the state vector X
t
conditional on the observations Y
1:t1
= y
1:t1
and the parameter =. It is a well
known result [1], [2], [4] that the sequence of probability measures
t,
(x
t
)dx
t
can be recursively approximated
using a standard PF algorithm, and hence the integral (g
y
t
,
t,
) can be numerically approximated as well.
In next section, we outline a novel PF algorithm that enables the recursive approximation of the pdfs µ
t
,
t =1, 2,..., and produces Bayesian estimates of the state variables X
t
, t =1, 2,..., as a by-product.
III. NESTED PARTICLE FILTERING SCHEME
A. Standard particle filter
Assume = are given parameters. The standard particle filter is a recursive Monte Carlo algorithm for
the approximation of the sequence of predictive probability measures
t,
(x
t
)dx
t
on the state space X, as well
as the associated filtering probability measures
t,
(x
t
)dx
t
, where
t,
is the density of X
t
conditional on given
observations Y
1:t
= y
1:t
and the parameter .
At time t =0, we generate M samples (termed particles), x
(i)
0
0
, i =1,...,M, from the prior density
0
. At
every time t>0, we apply the algorithm below, where
x
0
(dx) denotes the unit delta measure located at x
0
.
Algorithm 1: We take as inputs the parameter , the observation y
t
and the particle approximation of
t1,
(x
t
)dx
t
at time t 1,
M
t1,
(dx
t
)=
1
M
P
M
i=1
x
(i)
t1
(dx
t
).
Computations:
Generate new particles ¯x
(i)
t
(x
t
|x
(i)
t1
) and compute normalised importance weights w
(i)
t
/ g
y
t
x
(i)
t
),
i =1,...,M.
Resample: for i =1,...,M, assign x
(i)
t
x
(j)
t
with probability w
(j)
t
, j 2{1,...,M}.
Outputs: New particle approximations
M
t,
(dx
t
)=
1
M
P
M
i=1
¯x
(i)
t
(dx
t
),
M
t,
(dx
t
)=
1
M
P
M
i=1
x
(i)
t
(dx
t
), and
M
t,
(y
t
)=(g
y
t
,
M
t,
)=
1
M
P
M
i=1
g
y
t
x
(i)
t
).
Given any bounded function f 2 B(X) and any p 1, it can be proved [7] under mild assumptions that
k(f,
M
t,
) (f,
t,
)k
p
c
t,
p
M
,
where c
t,
is constant w.r.t. M. Similar convergence results hold for (f,
M
t,
) and
M
t,
(y
t
).

5
B. Proposed algorithm
We outline a PF scheme for the recursive approximation of the sequence of probability measures µ
t
()d,
t =1, 2,.... See [4] for full details. It is a nested Monte Carlo scheme, where Algorithm 1 is used to compute
importance weights for particles in the parameter space S. We assume that the latter is compact. In particular, we
select S =[f
, f
+
] [a
, a
+
] [a
, a
+
] R
3
for some known and finite bounds f
< f
+
and a
< a
+
.
Algorithm 2: At time t =0, draw N i.i.d. particles
(i)
0
from µ
0
() and NM i.i.d. particles x
(i,j)
0
from
0
(x
0
),
i =1,...,N and j =1,...,M. Let µ
N
0
=
1
N
P
N
i=1
(i)
0
and
M
0,
(i)
0
=
1
M
P
M
j=1
x
(i,j)
0
. Choose a conditional pdf
N
(|
0
) on the parameter space S that satisfies, for at least some p 1,
sup
0
2S
Z
k
0
k
p
(|
0
)d
c
p
N
p
2
. (6)
where c
is some constant independent of N.
Computations: Given µ
N
t1
=
1
N
P
N
i=1
(i)
t1
, the N sets
i,M
t1
= {x
(i,j)
t1
}
M
j=1
and the new observation y
t
, take the
following steps at time t.
Draw N new particles
¯
(i)
t
N
(|
(i)
t1
), i =1,...,N.
Run Algorithm 1 N times, with inputs =
¯
(i)
t
and
i,M
t1
= {x
(i,j)
t1
}
M
j=1
, to obtain
normalised importance weights:
(i)
t
/
M
t,
¯
(i)
t
(y
t
),
updated particles in X: ¯
i,M
t
= {x
(i,j)
t
}
M
j=1
,
for i =1,...,N.
Resample: for i =1,...,N, assign
(i)
t
=
¯
(l)
t
and
i,M
t
l,M
t
with probability
(l)
t
, l 2{1,...,N}.
Outputs: The approximations µ
N
t
(d)=
1
N
P
N
i=1
(i)
t
(d) and sets
i,M
t
= {x
(i,j)
t
}
M
j=1
, for i =1,...,N.
Given µ
N
t
and the collection of sets
i,M
t
, i =1,...,N, it is straightforward to obtain posterior estimates of the
parameters and the dynamic variables, e.g.,
E[|y
1:t
]
1
N
X
(i)
t
,E[X
t
|y
1:t
]
1
NM
N,M
X
i=1,j=1
x
(i,j)
t
.
The inequality (6) simply states that the pdf
N
(|
0
) should have a sufficiently small variance (in all directions).
A simple truncated Gaussian kernel with suitable variance, e.g.,
N
(|
0
) / I
S
()N
;
0
,
c
N
I
3
,
where I
S
()=1if 2 S and 0 otherwise, suffices to make (6) hold for all p 1 (and guarantee convergence [4]).
IV. KERNEL DENSITY ESTIMATORS
We aim at approximating the posterior pdf of the parameters F and A 2 R
2
. If =[F, A
>
]
>
, then we denote
µ
t
()=µ
t
(f, a) and calculate the posterior marginal densities as
µ
t,F
(f)=
Z
µ
t
(f, a)da
t,A
(a)=
Z
µ
t
(f, a)df.
Similary, the particle approximation µ
N
t
(d)=µ
N
t
(df da) yields the approximate marginal probability measures
µ
N
t,F
(df)=
Z
[a
,a
+
]
2
µ
N
t
(df dA)=
1
N
N
X
i=1
f
(i)
t
(df),
µ
N
t,A
(da)=
Z
[f
,f
+
]
µ
N
t
(df dA)=
1
N
N
X
i=1
a
(i)
t
(da),
where we have used the obvious notation
(i)
t
=[f
(i)
t
, a
(i)
>
t
]
>
.

Citations
More filters
Journal ArticleDOI
TL;DR: A state dimension reduction method is proposed to overcome the particle degeneracy problem of particle filter which is used to fusion GNSS and 5G D2D measurements and shows that the proposed integrated methodology outperforms the nonintegrated one.
Abstract: Global navigation satellite system (GNSS) is not suitable for the dense urban or indoor environments as the satellite signals are very weak. Meanwhile, positioning is an important application of the fifth-generation (5G) communication system. GNSS/5G integrated positioning system becomes a promising research topic with the development of 5G standard. This paper focuses on the integrated methodology of GNSS and device to device (D2D) measurements in 5G communication system. We analyze the characteristics of this type of integrated system and propose a high-efficiency D2D positioning measure protocol, named crossover multiple-way ranging, which consumes less communication resources. Then, to deal with the high-dimensional state space in the integrated system, a state dimension reduction method is proposed to overcome the particle degeneracy problem of particle filter which is used to fusion GNSS and 5G D2D measurements. Three integrated algorithms in different scenarios have been proposed: the first one is the integrated algorithm when the range measurements can be measured directly. The second one is the integrated algorithm with unknown time skew and offset of each mobile terminal. The third one is the integrated algorithm in GNSS-denied environment which is prevalent in urban and indoor applications. The simulation and experimental results show that our proposed integrated methodology outperforms the nonintegrated one.

65 citations


Cites background from "Particle filtering for Bayesian par..."

  • ...Thus, the state space is a high-dimensional variable which means it needs huge amount of particles to draw the joint distribution [21]....

    [...]

Posted Content
01 Jan 2013
TL;DR: The SMC^2 algorithm proposed in this paper is a sequential Monte Carlo algorithm, defined in the theta-dimension, which propagates and resamples many particle filters in the x-dimension.
Abstract: We consider the generic problem of performing sequential Bayesian inference in a state-space model with observation process y, state process x and fixed parameter theta. An idealized approach would be to apply the iterated batch importance sampling (IBIS) algorithm of Chopin (2002). This is a sequential Monte Carlo algorithm in the theta-dimension, that samples values of theta, reweights iteratively these values using the likelihood increments p(y_t|y_1:t-1, theta), and rejuvenates the theta-particles through a resampling step and a MCMC update step. In state-space models these likelihood increments are intractable in most cases, but they may be unbiasedly estimated by a particle filter in the x-dimension, for any fixed theta. This motivates the SMC^2 algorithm proposed in this article: a sequential Monte Carlo algorithm, defined in the theta-dimension, which propagates and resamples many particle filters in the x-dimension. The filters in the x-dimension are an example of the random weight particle filter as in Fearnhead et al. (2010). On the other hand, the particle Markov chain Monte Carlo (PMCMC) framework developed in Andrieu et al. (2010) allows us to design appropriate MCMC rejuvenation steps. Thus, the theta-particles target the correct posterior distribution at each iteration t, despite the intractability of the likelihood increments. We explore the applicability of our algorithm in both sequential and non-sequential applications and consider various degrees of freedom, as for example increasing dynamically the number of x-particles. We contrast our approach to various competing methods, both conceptually and empirically through a detailed simulation study, included here and in a supplement, and based on particularly challenging examples.

50 citations

01 Apr 2012
TL;DR: In this article, the authors study numerical methods available for estimating closure parameters in chaotic models and discuss three techniques: off-line likelihood calculations using filtering methods, the state augmentation method, and the approach that utilizes summary statistics from long model simulations.
Abstract: Abstract. Many dynamical models, such as numerical weather prediction and climate models, contain so called closure parameters. These parameters usually appear in physical parameterizations of sub-grid scale processes, and they act as "tuning handles" of the models. Currently, the values of these parameters are specified mostly manually, but the increasing complexity of the models calls for more algorithmic ways to perform the tuning. Traditionally, parameters of dynamical systems are estimated by directly comparing the model simulations to observed data using, for instance, a least squares approach. However, if the models are chaotic, the classical approach can be ineffective, since small errors in the initial conditions can lead to large, unpredictable deviations from the observations. In this paper, we study numerical methods available for estimating closure parameters in chaotic models. We discuss three techniques: off-line likelihood calculations using filtering methods, the state augmentation method, and the approach that utilizes summary statistics from long model simulations. The properties of the methods are studied using a modified version of the Lorenz 95 system, where the effect of fast variables are described using a simple parameterization.

32 citations

16 Apr 2022
TL;DR: A new recursive methodology for Bayesian inference that aims at estimating the static parameters and tracking the dynamic variables ofMulti-scale systems, where variables of interest evolve in different time-scales and live in different state-spaces is introduced.
Abstract: Multi-scale problems, where variables of interest evolve in different time-scales and live in different state-spaces. can be found in many fields of science. Here, we introduce a new recursive methodology for Bayesian inference that aims at estimating the static parameters and tracking the dynamic variables of these kind of systems. Although the proposed approach works in rather general multiscale systems, for clarity we analyze the case of a heterogeneous multi-scale model with 3 time-scales (static parameters, slow dynamic state variables and fast dynamic state variables). The proposed scheme, based on nested filtering methodology of [26], combines three intertwined layers of filtering techniques that approximate recursively the joint posterior probability distribution of the parameters and both sets of dynamic state variables given a sequence of partial and noisy observations. We explore the use of sequential Monte Carlo schemes in the first and second layers while we use an unscented Kalman filter to obtain a Gaussian approximation of the posterior probability distribution of the fast variables in the third layer. Some numerical results are presented for a stochastic two-scale Lorenz 96 model with unknown parameters.
References
More filters
Journal ArticleDOI
TL;DR: It is shown here how it is possible to build efficient high dimensional proposal distributions by using sequential Monte Carlo methods, which allows not only to improve over standard Markov chain Monte Carlo schemes but also to make Bayesian inference feasible for a large class of statistical models where this was not previously so.
Abstract: Summary. Markov chain Monte Carlo and sequential Monte Carlo methods have emerged as the two main tools to sample from high dimensional probability distributions. Although asymptotic convergence of Markov chain Monte Carlo algorithms is ensured under weak assumptions, the performance of these algorithms is unreliable when the proposal distributions that are used to explore the space are poorly chosen and/or if highly correlated variables are updated independently. We show here how it is possible to build efficient high dimensional proposal distributions by using sequential Monte Carlo methods. This allows us not only to improve over standard Markov chain Monte Carlo schemes but also to make Bayesian inference feasible for a large class of statistical models where this was not previously so. We demonstrate these algorithms on a non-linear state space model and a Levy-driven stochastic volatility model.

1,869 citations


"Particle filtering for Bayesian par..." refers background or methods in this paper

  • ...It is a well known result [1, 2, 4] that the sequence of probability measures ξt,θ(xt)dxt can be recursively approximated using a standard PF algorithm, and hence the integral (gt , ξt,θ) can be numerically approximated as well....

    [...]

  • ..., the popular particle Markov chain Monte Carlo (pMCMC) [1] and sequential Monte Carlo square (SMC(2)) [2] algorithms, are batch techniques and, therefore, they are not well suited to the processing of long observation...

    [...]

Journal ArticleDOI
TL;DR: In this paper, a fully nonlinear particle filter is proposed for higher dimensional problems by exploiting the freedom of the proposal density inherent in particle filtering, which can be applied to high dimensional problems.
Abstract: Almost all research fields in geosciences use numerical models and observations and combine these using data-assimilation techniques. With ever-increasing resolution and complexity, the numerical models tend to be highly nonlinear and also observations become more complicated and their relation to the models more nonlinear. Standard data-assimilation techniques like (ensemble) Kalman filters and variational methods like 4D-Var rely on linearizations and are likely to fail in one way or another. Nonlinear data-assimilation techniques are available, but are only efficient for small-dimensional problems, hampered by the so-called ‘curse of dimensionality’. Here we present a fully nonlinear particle filter that can be applied to higher dimensional problems by exploiting the freedom of the proposal density inherent in particle filtering. The method is illustrated for the three-dimensional Lorenz model using three particles and the much more complex 40-dimensional Lorenz model using 20 particles. By also applying the method to the 1000-dimensional Lorenz model, again using only 20 particles, we demonstrate the strong scale-invariance of the method, leading to the optimistic conjecture that the method is applicable to realistic geophysical problems. Copyright c � 2010 Royal

400 citations


"Particle filtering for Bayesian par..." refers methods in this paper

  • ...The latter displays the basic physical features of atmospheric dynamics and, for this reason, the deterministic version of this model is commonly used as a benchmark for data assimilation [5] and parameter estimation techniques [6] in meteorology and climate science....

    [...]

Journal ArticleDOI
08 Nov 2004
TL;DR: A detailed overview of particle methods, a set of powerful and versatile simulation-based methods to perform optimal state estimation in nonlinear non-Gaussian state-space models, is provided.
Abstract: Particle methods are a set of powerful and versatile simulation-based methods to perform optimal state estimation in nonlinear non-Gaussian state-space models. The ability to compute the optimal filter is central to solving important problems in areas such as change detection, parameter estimation, and control. Much recent work has been done in these areas. The objective of this paper is to provide a detailed overview of them.

352 citations


"Particle filtering for Bayesian par..." refers background in this paper

  • ...Although some recursive algorithms exist [3], they only yield maximum likelihood (point) estimates for the parameters, and hence they are subject to various convergence (and complexity) issues when the likelihood is multimodal, contains singularities or cannot be computed exactly....

    [...]

Posted Content
TL;DR: The SMC2 algorithm is proposed, a sequential Monte Carlo algorithm, defined in the θ‐dimension, which propagates and resamples many particle filters in the x‐ dimension, which explores the applicability of the algorithm in both sequential and non‐sequential applications and considers various degrees of freedom.
Abstract: We consider the generic problem of performing sequential Bayesian inference in a state-space model with observation process y, state process x and fixed parameter theta. An idealized approach would be to apply the iterated batch importance sampling (IBIS) algorithm of Chopin (2002). This is a sequential Monte Carlo algorithm in the theta-dimension, that samples values of theta, reweights iteratively these values using the likelihood increments p(y_t|y_1:t-1, theta), and rejuvenates the theta-particles through a resampling step and a MCMC update step. In state-space models these likelihood increments are intractable in most cases, but they may be unbiasedly estimated by a particle filter in the x-dimension, for any fixed theta. This motivates the SMC^2 algorithm proposed in this article: a sequential Monte Carlo algorithm, defined in the theta-dimension, which propagates and resamples many particle filters in the x-dimension. The filters in the x-dimension are an example of the random weight particle filter as in Fearnhead et al. (2010). On the other hand, the particle Markov chain Monte Carlo (PMCMC) framework developed in Andrieu et al. (2010) allows us to design appropriate MCMC rejuvenation steps. Thus, the theta-particles target the correct posterior distribution at each iteration t, despite the intractability of the likelihood increments. We explore the applicability of our algorithm in both sequential and non-sequential applications and consider various degrees of freedom, as for example increasing dynamically the number of x-particles. We contrast our approach to various competing methods, both conceptually and empirically through a detailed simulation study, included here and in a supplement, and based on particularly challenging examples.

258 citations


"Particle filtering for Bayesian par..." refers background or methods in this paper

  • ...It is a well known result [1, 2, 4] that the sequence of probability measures ξt,θ(xt)dxt can be recursively approximated using a standard PF algorithm, and hence the integral (gt , ξt,θ) can be numerically approximated as well....

    [...]

  • ..., the popular particle Markov chain Monte Carlo (pMCMC) [1] and sequential Monte Carlo square (SMC(2)) [2] algorithms, are batch techniques and, therefore, they are not well suited to the processing of long observation...

    [...]

  • ...1: (a)-(b) Estimates of μt,F(f) over the interval [2, 30]....

    [...]

Frequently Asked Questions (1)
Q1. What have the authors contributed in "Particle filtering for bayesian parameter estimation in a high dimensional state space model" ?

In this paper, the authors investigate the application of a novel nested particle filtering scheme for joint Bayesian parameter estimation and tracking of the dynamic variables in a high dimensional state space model –namely a stochastic version of the two-scale Lorenz 96 chaotic system, commonly used as a benchmark model in meteorology and climate science. The authors provide theoretical guarantees on the algorithm performance, including uniform convergence rates for the approximation of posterior probability density functions of the fixed model parameters.