# Maximum Likelihood Recursive State Estimation using the Expectation
Maximization Algorithm

Mohammad S. Ramadan<EMAIL_ADDRESS>Robert R. Bitmead
<EMAIL_ADDRESS>Department of Mechanical & Aerospace Engineering,
University of California, San Diego, La Jolla CA 92093-0411, USA, and
Electrical & Computer Engineering Department, San Diego State University, San
Diego CA 92182, USA  Department of Mechanical & Aerospace Engineering,
University of California, San Diego, La Jolla CA 92093-0411, USA.

###### Abstract

A Maximum Likelihood recursive state estimator is derived for non-linear and
non-Gaussian state-space models. The estimator combines a particle filter to
generate the conditional density and the Expectation Maximization algorithm to
compute the maximum likelihood state estimate iteratively. Algorithms for
maximum likelihood state filtering, prediction and smoothing are presented.
The convergence properties of these algorithms, which are inherited from the
Expectation Maximization algorithm, are proven and examined in two examples.
It is shown that, with randomized reinitialization, which is feasible because
of the algorithm simplicity, these methods are able to converge to the Maximum
Likelihood Estimate (MLE) of multimodal, truncated and skewed densities, as
well as those of disjoint support.

###### keywords:

Maximum Likelihood, State Estimation, Expectation Maximization, Particle
Filter .

,

##  1 Introduction

We derive a maximum likelihood recursive state estimator for nonlinear systems
based on: knowledge of the system equations and process noise and measurement
noise densities; the operation of a particle filter to propagate conditional
densities; and the Expectation Maximization (EM) algorithm [2, 14]. If the
state at time $k$ is denoted $\zeta_{k}$, and given the sequence of outputs up
to time $k$, namely $\mathcal{Y}_{0:k}=\\{y_{0},y_{1},\dots,y_{k}\\}$, the
algorithm provides an estimate of the maximizer of the posterior state density
$p(\zeta_{k}|\mathcal{Y}_{0:k})$. This contrasts with Least-Squares
estimators, which produce the conditional mean of this density. The
contribution of the paper is to provide an easily implementable algorithm
which is recursive and which yields this Maximum Likelihood Estimate (MLE). We
propagate the conditional densities, both prior and posterior, using a
particle filter, and the EM algorithm is used to compute the MLE of the
posterior density from the particles of the prior density, the current
measurement and the known noise densities. The observation is that knowledge
of the noise densities permits computation of the gradient of successive
functions whose maximizers are non-decreasing in their likelihood with respect
to the posterior density, which follows from the desirable convergence
properties of the EM algorithm. If the system under consideration has Gaussian
additive noises and linear measurement equation, these successive functions
are quadratic and the algorithm simplifies into a fixed-point iteration. The
paper clarifies these details in relation to the particle filter, the EM
algorithm and MLE.

In many scenarios, the conditional state density can: be multi-modal as when
tracking a group of objects [13, 19]; have skewness, limited support, one-
sidedness as in some Stochastic Volatility Models [20, 9]; or be governed by
inequality state constraints [8], leading to truncated or disjoint support
such as might arise in stochastic model Predictive Control. The justification
for the choice of the conditional mean as the point estimate of a filter may
be weakened or render an estimate which is infeasible under the density. In
such scenarios, the MLE may have more reasonable grounds for use. Since
finding the MLE requires optimization over the posterior density, the
recursive procedure involves two stages: the propagation of the posterior
density, for which we rely on a particle filter; and, an iterative
optimization to find the maximizer. However, the density provided by the
particle filter is a collection of point masses, usually with the same
singleton mass, and so is not suited to Newton-like optimization. A novelty of
this paper and its algorithms is that the previous time’s density plus the
current measurement are combined with the known continuous densities of the
process and measurement noises to yield a gradient based maximizer.
Effectively, this procedure yields a smooth interpolation of the particle
filter’s posterior density and of its gradient.

In [11], a piecewise continuous cumulative distribution function is
constructed from the weighted particles approximating the filtered state
density $p(\zeta_{k}|\mathcal{Y}_{0:k})$, limited to the case when the state-
space is one-dimensional. In [1], this solution is extended to a multi-
dimensional state-space. However, it has a computational complexity of
$\mathcal{O}(N^{2})$, where $N$ is the number of particles, and a nonstandard
particle filtering scheme, as reported in [4]. Earlier, in the construction of
the auxiliary particle filter [15], it was shown that it is possible to
construct a smooth approximation of $p(\zeta_{k}|\mathcal{Y}_{0:k})$ under the
name of “the empirical filtered density” using the weighted particles
approximating $p(\zeta_{k-1}|\mathcal{Y}_{0:k-1})$ and the current
measurement.We extend this construction with the computation of an approximate
gradient, as shall be explained shortly. This constructed smooth density could
be directly maximized using various search methods and the new calculation
permits gradient search for the MLE. We introduce the EM algorithm as an
iterative procedure to achieve this, perhaps with multiple initial starting
values in the multimodal case. This yields a sequence of estimates with non-
decreasing likelihood which converges to a stationary point [22, 2]. This will
be supported by computational examples in Sections 3.4 and 3.5. In effect, the
EM algorithm smooths and recenters the empirical filtered density at the
E-step before maximizing at the M-step. The iterations refine and localize the
smoothing to approach the MLE. Ultimately, this is limited by the quality of
the empirical density generated by the particle filter. Without such
smoothing, it is not possible to find the maximizer of the particle filter’s
point-mass density.

The EM algorithm [2] is versatile with many variations [12, 21, 5] and
applications primarily in parameter estimation. It consists of two major
steps: the expectation (E) step, in which an approximant to the log-likelihood
function is constructed, followed by the maximization (M) step, which seeks to
find the maximizer of this function. Hence, EM iterations deliver a sequence
of functions for which their maximizers are non-decreasing with respect to the
likelihood function under consideration [22]. Equipped with Monte-Carlo
particle methods, the EM algorithm is widely used in parameter estimation for
nonlinear non-Gaussian systems [18, 16, 10, 6]. We examine the use of EM in
the context of state estimation; to our knowledge this is a novel application
and different in principle, since, unlike parameters, the components of the
state are dynamically changing according to the tailored system model.

The particle filter [3], itself, is a common feature of nonlinear estimation
and operates by propagating the conditional densities as characterized by a
collection of samples, called particles, distributed with this density. For a
basic particle filter, the time update pass, i.e. posterior-to-prior
densities, uses the state equation and samples from the process noise density.
This is followed by the measurement update, prior-to-posterior densities,
relying on weighting and possibly resampling these particles based on the
output equation and the output noise density.

In this paper, we use the weighted particles’ approximation of
$p(\zeta_{k-1}|\mathcal{Y}_{0:k-1})$ and the current observation $y_{k}$ with
the EM algorithm to construct a sequence of functions of $\zeta_{k}$, whose
maximizers are non-decreasing in their likelihood with respect to
$p(\zeta_{k}|\mathcal{Y}_{0:k})$. In Section 2, the underlying assumptions on
the state-space model under consideration are elucidated. In Section 3, the EM
state filtering (EMSF) algorithm is introduced, and implemented in two
examples 3.4 and 3.5. First, it is applied to state filtering of a linear
Gaussian state-space model, and since the conditional mean is equivalent to
the conditional mode of such systems [17], it is compared to the Kalman filter
solution, and shown to be convergent to the conditional mode (the global
maximizer) from a random initial guess. This establishes the efficacy of the
algorithm and permits statements concerning its complexity. In the second
example, it is used to track the MLE of a multi-modal state distribution of a
nonlinear state-space model. It is also shown, in Section 3.6 that this new EM
state filter (EMSF) algorithm reduces to a convergent fixed-point iteration
method under mild assumption. The extensions of this algorithm to the contexts
of state prediction and smoothing are provided in Sections 4 and 5.

##  2 Problem Statement

Consider the following nonlinear dynamical system affine in the noises111The
notation $\sim$ here means “is distributed as the density on the right”.:

$\displaystyle\zeta_{k+1}$
$\displaystyle=f(\zeta_{k})+w_{k},\,\zeta_{k}\in\mathbb{R}^{r_{\zeta}},$ (1a)
$\displaystyle y_{k}$
$\displaystyle=g(\zeta_{k})+v_{k},\,y_{k}\in\mathbb{R}^{r_{y}},$ (1b)
$\displaystyle w_{k}$ $\displaystyle\sim W_{k}(.),\,v_{k}\sim V_{k}(.),$
$\displaystyle\zeta_{0}$ $\displaystyle\sim p_{0}(.)$

where,

  * •

the stochastic processes $\\{w_{k}\\}$ and $\\{v_{k}\\}$ are white and
independent, and their probability density functions, $W_{k}(\cdot)$ and
$V_{k}(\cdot)$ are known and continuously differentiable for all $k$,

  * •

$\zeta_{k}$ is the state vector, with initial value, $\zeta_{0}$, having a
known initial distribution which is independent of $\\{w_{k}\\}$ and
$\\{v_{k}\\}$,

  * •

$y_{k}$ is the measurement vector,

  * •

the functions $f$ and $g$ are known continuously differentiable functions.

Given the observation sequence,
$\mathcal{Y}_{0:k}=\\{y_{0},y_{1},\dots,y_{k}\\}$, the problem is to find the
MLE of the state-vector $\zeta_{k}$ using the log-likelihood function $\log
p(\zeta_{k}|\mathcal{Y}_{0:m})$, for $m=k$ (filtering), $m>k$ (smoothing), and
$m<k$ (prediction).

##  3 ML state filtering by the expectation maximization algorithm

Given, at time $k$:

  * •

the problem statement and assumptions as in Section 2,

  * •

the measurement $y_{k}$,

  * •

the filtered density at time $k-1$, i.e. $p(\zeta_{k-1}|\mathcal{Y}_{0:k-1})$,

  * •

initial time-$k$ state estimate $\zeta^{0}_{k}$,

we set iteration counter $i=0$ and seek to find $\zeta_{k}^{i+1}$ such that
$\log p(\zeta_{k}^{i+1}|\mathcal{Y}_{0:k})\geq\log
p(\zeta_{k}^{i}|\mathcal{Y}_{0:k})$ [22]. This is achieved by two major steps:
the expectation step (E-step), followed by the maximization step (M-step). We
continue with these iterations in $i$ until we satisfy a convergence
criterion, which is guaranteed by the non-decreasing property. At this stage,
we either restart with a different initial condition, $\zeta_{k}^{0}$, or move
to time $k+1$.

###  3.1 E-step

The E-step is performed by evaluating the following expectation,

$\displaystyle Q_{\zeta}(\zeta_{k},\zeta_{k}^{i})$
$\displaystyle=\mathbb{E\,}_{\zeta_{k-1}|\zeta_{k}^{i},\mathcal{Y}_{0:k}}\\{\log
p(\zeta_{k}|\zeta_{k-1},\mathcal{Y}_{0:k})|\zeta_{k}^{i},\mathcal{Y}_{0:k}\\}$
$\displaystyle\hskip-42.67912pt=\int\log
p(\zeta_{k}|\zeta_{k-1},\mathcal{Y}_{0:k})p(\zeta_{k-1}|\zeta_{k}^{i},\mathcal{Y}_{0:k})\,d\zeta_{k-1}$
(2)

whose calculation in terms of available densities will be presented shortly in
the following theorem.

Before presenting the theorem, we first derive a lemma detailing a conditional
independence result, which will be used repeatedly in the proofs and
derivations.

######  Lemma 1.

For system (1) at time $n\geq k$, $y_{n}$ conditioned on $\zeta_{k}$ is
conditionally independent from both
$\zeta_{0:k-1}=\\{\zeta_{0},\zeta_{1},\ldots,\zeta_{k-1}\\}$ and
$\mathcal{Y}_{0:k-1}$, that is

$\displaystyle
p(y_{n}|\zeta_{k},\zeta_{0:k-1},\mathcal{Y}_{0:k-1})=p(y_{n}|\zeta_{k})$

Proof. From (1) and for $n\geq k$, $y_{n}$ is a function of $\zeta_{k}$,
$w_{k:n-1}$ and $v_{n}$. Since $w_{k:n-1}$ and $v_{n}$ are independent from
$\zeta_{0:k-1}$ and $\mathcal{Y}_{0:k-1}$, we have the result. $\square$

######  Theorem 1.

Given the dynamical system (1), the filtered density
$p(\zeta_{k-1}|\mathcal{Y}_{0:k-1})$, the current measurement $y_{k}$, and
$\zeta_{k}^{i}$, the integral in (2) is, apart from a positive multiplicative
constant and an additive constant both independent from $\zeta_{k}$, equal to

$\displaystyle Q_{\zeta}(\zeta_{k},\zeta_{k}^{i})=$
$\displaystyle\int\Big{[}\log V_{k}(y_{k}-g(\zeta_{k}))+\log
W_{k-1}(\zeta_{k}-f(\zeta_{k-1}))\Big{]}\times$ $\displaystyle
W_{k-1}(\zeta_{k}^{i}-f(\zeta_{k-1}))p(\zeta_{k-1}|\mathcal{Y}_{0:k-1})\,d\zeta_{k-1},$
(3)

and its gradient is given by

$\displaystyle\nabla_{\zeta_{k}}Q_{\zeta}(\zeta_{k},\zeta_{k}^{i})$
$\displaystyle=\int\Big{[}\nabla_{\zeta_{k}}\log V_{k}(y_{k}-g(\zeta_{k}))+$
$\displaystyle\hskip-28.45274pt\nabla_{\zeta_{k}}\log
W_{k-1}(\zeta_{k}-f(\zeta_{k-1}))\Big{]}\times$
$\displaystyle\hskip-28.45274ptW_{k-1}(\zeta_{k}^{i}-f(\zeta_{k-1}))p(\zeta_{k-1}|\mathcal{Y}_{0:k-1})d\zeta_{k-1}.$
(4)

Proof. We use Bayes rule and Lemma 1, to yield

$\displaystyle p(\zeta_{k}|\zeta_{k-1},\mathcal{Y}_{0:k})$
$\displaystyle=\frac{p(y_{k}|\zeta_{k},\zeta_{k-1},\mathcal{Y}_{0:k-1})p(\zeta_{k}|\zeta_{k-1},\mathcal{Y}_{0:k-1})}{p(y_{k}|\zeta_{k-1},\mathcal{Y}_{0:k-1})}$
$\displaystyle=\frac{p(y_{k}|\zeta_{k})p(\zeta_{k}|\zeta_{k-1})}{p(y_{k}|\zeta_{k-1})}.$

Thus,

$\displaystyle\log p(\zeta_{k}|\zeta_{k-1},\mathcal{Y}_{0:k})$
$\displaystyle=\log p(y_{k}|\zeta_{k})+\log p(\zeta_{k}|\zeta_{k-1})$
$\displaystyle\hskip-28.45274pt-\log p(y_{k}|\zeta_{k-1})$ $\displaystyle=\log
V_{k}(y_{k}-g(\zeta_{k}))$ $\displaystyle\hskip-28.45274pt+\log
W_{k-1}(\zeta_{k}-f(\zeta_{k-1}))-\tilde{C}$ (5)

where $\tilde{C}=\log p(y_{k}|\zeta_{k-1})$ is independent from $\zeta_{k}$
and so integrates to an additive constant in the expectation (2).

Similarly, by Bayes rule and Lemma 1,

$\displaystyle p(\zeta_{k-1}|\zeta_{k}^{i},\mathcal{Y}_{0:k})$
$\displaystyle=p(\zeta_{k-1}|\zeta_{k}^{i},\mathcal{Y}_{0:k-1},y_{k}),$
$\displaystyle=\frac{p(y_{k}|\zeta_{k}^{i},\zeta_{k-1},\mathcal{Y}_{0:k-1})p(\zeta_{k-1}|\zeta_{k}^{i},\mathcal{Y}_{0:k-1})}{p(y_{k}|\zeta_{k}^{i},\mathcal{Y}_{0:k-1})},$
$\displaystyle=\frac{p(y_{k}|\zeta_{k}^{i})p(\zeta_{k-1}|\zeta_{k}^{i},\mathcal{Y}_{0:k-1})}{p(y_{k}|\zeta_{k}^{i})},$
$\displaystyle=p(\zeta_{k-1}|\zeta_{k}^{i},\mathcal{Y}_{0:k-1}),$

and, applying Bayes rule again,

$\displaystyle p(\zeta_{k-1}|\zeta_{k}^{i},\mathcal{Y}_{0:k})$
$\displaystyle=\frac{p(\zeta_{k}^{i}|\zeta_{k-1})p(\zeta_{k-1}|\mathcal{Y}_{0:k-1})}{p(\zeta_{k}^{i}|\mathcal{Y}_{0:k-1})}.$

Hence,

$\displaystyle p(\zeta_{k-1}|\zeta_{k}^{i},\mathcal{Y}_{0:k})$
$\displaystyle=$
$\displaystyle\hskip-14.22636pt\frac{W_{k-1}(\zeta_{k}^{i}-f(\zeta_{k-1}))p(\zeta_{k-1}|\mathcal{Y}_{0:k-1})}{\hat{C}}$
(6)

where $\hat{C}=p(\zeta_{k}^{i}|\mathcal{Y}_{0:k-1})$ is independent from
$\zeta_{k-1}$ and leads to a positive multiplicative constant in the
expectation (2).

Substituting (5) and (6) into (2) and neglecting the additive and positive
multiplicative constants above yields

$\displaystyle Q_{\zeta}(\zeta_{k},\zeta_{k}^{i})$
$\displaystyle=\int\Big{\\{}\big{[}\log p(y_{k}|\zeta_{k})+\log
p(\zeta_{k}|\zeta_{k-1})\big{]}\times$ $\displaystyle\hskip
42.67912ptp(\zeta_{k}^{i}|\zeta_{k-1})p(\zeta_{k-1}|\mathcal{Y}_{0:k-1})\Big{\\}}\,d\zeta_{k-1}$
$\displaystyle\hskip-42.67912pt=\int\Big{\\{}\big{[}\log
V_{k}(y_{k}-g(\zeta_{k}))+\log W_{k-1}(\zeta_{k}-f(\zeta_{k-1}))\big{]}\times$
$\displaystyle
W_{k-1}(\zeta_{k}^{i}-f(\zeta_{k-1}))p(\zeta_{k-1}|\mathcal{Y}_{0:k-1})\Big{\\}}\,d\zeta_{k-1}.$

The smoothness assumptions in Section 2 permit exchanging differentiation and
integration to arrive at the gradient formula.

$\displaystyle\nabla_{\zeta_{k}}Q_{\zeta}(\zeta_{k},\zeta_{k}^{i})$
$\displaystyle=\nabla_{\zeta_{k}}\int\Big{\\{}\big{[}\log
V_{k}(y_{k}-g(\zeta_{k}))+$ $\displaystyle\hskip-28.45274pt\log
W_{k-1}(\zeta_{k}-f(\zeta_{k-1}))\big{]}\times$
$\displaystyle\hskip-28.45274ptW_{k-1}(\zeta_{k}^{i}-f(\zeta_{k-1}))p(\zeta_{k-1}|\mathcal{Y}_{0:k-1})\Big{\\}}\,d\zeta_{k-1}$
$\displaystyle=\int\Big{\\{}\big{[}\nabla_{\zeta_{k}}\log
V_{k}(y_{k}-g(\zeta_{k}))$
$\displaystyle\hskip-28.45274pt+\nabla_{\zeta_{k}}\log
W_{k-1}(\zeta_{k}-f(\zeta_{k-1}))\big{]}\times$
$\displaystyle\hskip-28.45274ptW_{k-1}(\zeta_{k}^{i}-f(\zeta_{k-1}))p(\zeta_{k-1}|\mathcal{Y}_{0:k-1})\Big{\\}}\,d\zeta_{k-1}.$

$\square$

The gradient $\nabla_{\zeta_{k}}Q_{\zeta}(\zeta_{k},\zeta_{k}^{i})$ can be
used within a gradient search algorithm to evaluate the M-step, which follows
immediately.

###  3.2 M-step

The M-step maximizes the expectation in the E-step w.r.t $\zeta_{k}$, choosing
$\zeta_{k}^{i+1}$ to be the maximizer,

$\displaystyle\zeta_{k}^{i+1}=\text{arg}\max_{\zeta_{k}\in\mathbb{R}^{r_{\zeta}}}Q(\zeta_{k},\zeta_{k}^{i})$
(7)

This choice of $\zeta_{k}^{i+1}$ will be shown in Theorem 2 to be non-
decreasing in its likelihood compared to $\zeta_{k}^{i}$. A proof of the
general statement of EM algorithm can be found in [22]. We present the proof
in the context of state filtering.

######  Theorem 2.

The sequence $\\{\zeta_{k}^{i}\\}_{i\geq 0}$, commencing at
$\zeta_{k}^{0}\in\mathbb{R}^{r_{\zeta}}$ and satisfying the recursion (7), has
the property that $\\{\log p(\zeta_{k}^{i}|\mathcal{Y}_{0:k})\\}_{i\geq 0}$ is
monotonically non-decreasing.

Proof. The log-likelihood function $\log p(\zeta_{k}|\mathcal{Y}_{0:k})$ can
be expressed as

$\displaystyle\log p(\zeta_{k}|\mathcal{Y}_{0:k})$ $\displaystyle=\log
p(\zeta_{k},\zeta_{k-1}|\mathcal{Y}_{0:k})-\log
p(\zeta_{k-1}|\mathcal{Y}_{0:k},\zeta_{k})$ $\displaystyle=\log
p(\zeta_{k}|\zeta_{k-1},\mathcal{Y}_{0:k})+\log
p(\zeta_{k-1}|\mathcal{Y}_{0:k})$ $\displaystyle\hskip 14.22636pt-\log
p(\zeta_{k-1}|\mathcal{Y}_{0:k},\zeta_{k}).$

Taking the expectation with $p(\zeta_{k-1}|\zeta_{k}^{i},\mathcal{Y}_{0:k})$
as in (2),

$\displaystyle\log p(\zeta_{k}|\mathcal{Y}_{0:k})$
$\displaystyle=Q_{\zeta}(\zeta_{k},\zeta_{k}^{i})+\mathbb{E}_{\zeta_{k-1}|\zeta_{k}^{i},\mathcal{Y}_{0:k}}\\{\log
p(\zeta_{k-1}|\mathcal{Y}_{0:k})\\}$ $\displaystyle\hskip
14.22636pt-\mathbb{E}_{\zeta_{k-1}|\zeta_{k}^{i},\mathcal{Y}_{0:k}}\\{\log
p(\zeta_{k-1}|\mathcal{Y}_{0:k},\zeta_{k})\\}.$

Hence, the log-likelihood of $\zeta_{k}^{i}$,

$\displaystyle\log p(\zeta_{k}^{i}|\mathcal{Y}_{0:k})$
$\displaystyle=Q_{\zeta}(\zeta_{k}^{i},\zeta_{k}^{i})+\mathbb{E}_{\zeta_{k-1}|\zeta_{k}^{i},\mathcal{Y}_{0:k}}\\{\log
p(\zeta_{k-1}|\mathcal{Y}_{0:k})\\}$ $\displaystyle\hskip
14.22636pt-\mathbb{E}_{\zeta_{k-1}|\zeta_{k}^{i},\mathcal{Y}_{0:k}}\\{\log
p(\zeta_{k-1}|\mathcal{Y}_{0:k},\zeta_{k}^{i})\\}.$

Subtracting these two equations,

$\displaystyle\log p(\zeta_{k}|\mathcal{Y}_{0:k})-\log
p(\zeta_{k}^{i}|\mathcal{Y}_{0:k})=Q_{\zeta}(\zeta_{k},\zeta_{k}^{i})-Q_{\zeta}(\zeta_{k}^{i},\zeta_{k}^{i})$
$\displaystyle+\mathbb{E}_{\zeta_{k-1}|\zeta_{k}^{i},\mathcal{Y}_{0:k}}\\{\log
p(\zeta_{k-1}|\mathcal{Y}_{0:k},\zeta_{k}^{i})-\log
p(\zeta_{k-1}|\mathcal{Y}_{0:k},\zeta_{k})\\}$

where the term,

$\displaystyle\mathbb{E}_{\zeta_{k-1}|\zeta_{k}^{i},\mathcal{Y}_{0:k}}\\{\log
p(\zeta_{k-1}|\mathcal{Y}_{0:k},\zeta_{k}^{i})-\log
p(\zeta_{k-1}|\mathcal{Y}_{0:k},\zeta_{k})\\}$
$\displaystyle=\int_{\mathbb{R}^{r_{\zeta}}}\frac{\log
p(\zeta_{k-1}|\mathcal{Y}_{0:k},\zeta_{k}^{i})}{\log
p(\zeta_{k-1}|\mathcal{Y}_{0:k},\zeta_{k})}p(\zeta_{k-1}|\zeta_{k}^{i},\mathcal{Y}_{0:k})d\zeta_{k-1}$
$\displaystyle\geq 0.$

Since the right-hand side is the Kullback-Leibler divergence metric between
these two densities, this quantity is non-negative. (See Example 7.7 in [7].)
Thus,

$\displaystyle\log p(\zeta_{k}|\mathcal{Y}_{0:k})-\log
p(\zeta_{k}^{i}|\mathcal{Y}_{0:k})\geq
Q_{\zeta}(\zeta_{k},\zeta_{k}^{i})-Q_{\zeta}(\zeta_{k}^{i},\zeta_{k}^{i}).$

It follows from the definition of $\zeta_{k}^{i+1}$ as the maximizer of
$Q_{\zeta}(\zeta_{k},\zeta_{k}^{i})$ that

$\displaystyle
Q_{\zeta}(\zeta_{k}^{i+1},\zeta_{k}^{i})-Q_{\zeta}(\zeta_{k}^{i},\zeta_{k}^{i})\geq
0.$

Therefore,

$\displaystyle\log p(\zeta_{k}^{i+1}|\mathcal{Y}_{0:k})-\log
p(\zeta_{k}^{i}|\mathcal{Y}_{0:k})\geq 0,$

which completes the proof. $\square$

######  Corollary 1.

If the log-likelihood function $\log p(\zeta_{k}|\mathcal{Y}_{0:k})$ is
unimodal with the MLE being the only stationary point, then
$\\{\zeta_{k}^{i}\\}_{i\geq 0}$ converges to this MLE.

This parallels Corollary 1 in [22].

The EM algorithm for ML state filtering at time-step $k$ follows.

Algorithm 1

1:Given the current measurement $y_{k}$, the prior filtered density
$p(\zeta_{k-1}|\mathcal{Y}_{0:k-1})$ and prior MLE $\zeta_{k-1}^{*}$, set
$i\leftarrow 0$ and let
$\zeta_{k}^{0}=\zeta_{k}^{i}=f(\zeta_{k-1}^{*})+\mathbb{E}w_{k}$;

2:Evaluate the integral:  $\displaystyle Q_{\zeta}(\zeta_{k},\zeta_{k}^{i})$
$\displaystyle=\int\Big{\\{}\Big{[}\log V_{k}(y_{k}-g(\zeta_{k}))+$
$\displaystyle\hskip 42.67912pt\log
W_{k-1}(\zeta_{k}-f(\zeta_{k-1}))\Big{]}\times$
$\displaystyle\hskip-28.45274ptW_{k-1}(\zeta_{k}^{i}-f(\zeta_{k-1}))p(\zeta_{k-1}|\mathcal{Y}_{0:k-1})\Big{\\}}\,d\zeta_{k-1}$
(8)

3:Evaluate $\zeta_{k}^{i+1}$:
$\zeta_{k}^{i+1}=\text{arg}\max_{\zeta_{k}\in\mathbb{R}^{r_{\zeta}}}Q_{\zeta}(\zeta_{k},\zeta_{k}^{i})$
(9)

4:If $\\{\zeta_{k}^{i}\\}_{i\geq 0}$ satisfies a convergence criterion: set
$\zeta_{k}^{*}=\zeta_{k}^{i}$ and terminate. Else: $i\leftarrow i+1$ and go to
Step 2.

As guaranteed by Theorem 2, Algorithm 1 generates a sequence of estimates
$\\{\zeta_{k}^{i}\\}_{i\geq 0}$ which is non-decreasing in its likelihoods.

Although the densities $V_{k}$ and $W_{k}$ are assumed to be known, the
integral (8) in Algorithm 1 requires the knowledge of
$p(\zeta_{k-1}|\mathcal{Y}_{0:k-1})$ as well. In general, outside the linear
Gaussian case, such an integral cannot be analytically evaluated. Hence, a
numerical approximation has to be considered. In this paper, we turn to using
a particle filter.

###  3.3 EM with a particle filter

The EM algorithm, so far, relies on the provision of the filtered density
$p(\zeta_{k-1}|\mathcal{Y}_{0:k-1})$ for the calculation at the E-step. A
particle filter [3] provides a version of this information through a set of
particles, i.e. a collection of values of
$\\{\tilde{\zeta}^{j}_{k-1}\\}_{j=1}^{N}$, and possibly their corresponding
importance weights $\\{\omega_{k-1|k-1}^{j}\\}_{j=1}^{N}$. The density is
approximated by

$p(\zeta_{k-1}|\mathcal{Y}_{0:k-1})\approx\sum_{j=1}^{N}\omega_{k-1|k-1}^{j}\delta(\zeta_{k-1}-\tilde{\zeta}^{j}_{k-1}),$
(10)

where $\delta(\cdot)$ is the delta function. If the particles are produced by
resampling then each of the importance weights $\omega_{k-1|k-1}^{j}=1/N$.
Using these particles in (10), the expectation
$Q_{\zeta}(\zeta_{k},\zeta_{k}^{i})$, in (8), can be replaced by its particle
approximation $\hat{Q}^{N}_{\zeta}(\zeta_{k},\zeta_{k}^{i})$. Algorithm 1 can
then be recast using the particles for the E-step as follows.

Algorithm 2 EM State Filter (EMSF) Algorithm

1:Given $y_{k}$, the particles $\\{\tilde{\zeta}^{j}_{k-1}\\}_{j=1}^{N}$,
their normalized importance weights $\\{\omega_{k-1|k-1}^{j}\\}_{j=1}^{N}$ and
$\zeta_{k-1}^{*}$, set $i\leftarrow 0$ and let
$\zeta_{k}^{i}=f(\zeta_{k-1}^{*})+\mathbb{E}w_{k}$;

2:Evaluate the summation:
$\displaystyle\hat{Q}^{N}_{\zeta}(\zeta_{k},\zeta_{k}^{i})$
$\displaystyle=\sum_{j=1}^{N}\Big{[}\log V_{k}(y_{k}-g(\zeta_{k}))$
$\displaystyle+\log
W_{k-1}(\zeta_{k}-f(\tilde{\zeta}_{k-1}^{j}))\Big{]}\times$ $\displaystyle
W_{k-1}(\zeta_{k}^{i}-f(\tilde{\zeta}_{k-1}^{j}))\omega_{k-1|k-1}^{j}$ (11)

3:Evaluate $\zeta_{k}^{i+1}$:
$\zeta_{k}^{i+1}=\text{arg}\max_{\zeta_{k}\in\mathbb{R}^{r_{\zeta}}}\hat{Q}^{N}_{\zeta}(\zeta_{k},\zeta_{k}^{i})$
(12)

4:If $\\{\zeta_{k}^{i}\\}_{i\geq 0}$ satisfies a convergence criterion: set
$\zeta_{k}^{*}=\zeta_{k}^{i}$ and terminate. Else: $i\leftarrow i+1$ and go to
Step 2.

We apply EMSF in the following computed examples.

###  3.4 Example 1: Linear Gaussian state-space model

In the linear Gaussian state-space model case, all random variables under
consideration are Gaussian. Hence, the conditional mean is equivalent to the
conditional mode [17]. In the following example, the Kalman filter, which is
designed to track the propagation of the conditional mean, is used in
comparison with the EMSF. It will be shown that the EMSF, starting from a
random initialization $\zeta_{k}^{0}$, converges to the Kalman Filter
solution, which is guaranteed by Corollary 1.

Consider the following model,

$\displaystyle\zeta_{k+1}$ $\displaystyle=F\zeta_{k}+w_{k},$ (13)
$\displaystyle y_{k}$ $\displaystyle=H\zeta_{k}+v_{k},$ $\displaystyle v_{k}$
$\displaystyle\sim\mathcal{N}(0,R),\,w_{k}\sim\mathcal{N}(0,S)$
$\displaystyle\zeta_{0}$
$\displaystyle\sim\mathcal{N}(\zeta_{0|-1},\Sigma_{0|-1})$

where the sequences $\\{w_{k}\\}$ and $\\{v_{k}\\}$ and the random variable
$\zeta_{0}$ follow the independence and whiteness assumptions as in Section 2.
The term $\mathcal{N}(\mu,\Sigma)$ denotes the multivariate normal density of
mean $\mu$ and covariance $\Sigma$, and $\mathcal{N}(x|\mu,\Sigma)$ denotes
the value of this density when evaluated at $x$.

We use the following particle approximation,

$p(\zeta_{k-1}|\mathcal{Y}_{0:k-1})\approx\sum_{j=1}^{N}\omega_{k-1|k-1}^{j}\delta(\zeta_{k-1}-\tilde{\zeta}^{j}_{k-1}).$
(14)

The summation in Algorithm 2, ignoring constant terms superfluous to the
maximization, can be simplified to,

$\displaystyle\hat{Q}^{N}_{\zeta}(\zeta_{k},\zeta_{k}^{i})$
$\displaystyle=\sum_{j=1}^{N}\Big{[}-(y_{k}-H\zeta_{k})^{T}R^{-1}(y_{k}-H\zeta_{k})$
$\displaystyle\hskip
14.22636pt-(\zeta_{k}-F\tilde{\zeta}_{k-1}^{j})^{T}S^{-1}(\zeta_{k}-F\tilde{\zeta}_{k-1}^{j})\Big{]}$
$\displaystyle\hskip 14.22636pt\times\lambda^{j}\omega_{k-1|k-1}^{j}$ (15)

where
$\lambda^{j}=W_{k-1}(\zeta_{k}^{i}-F\tilde{\zeta}_{k-1}^{j})=\mathcal{N}(\zeta_{k}^{i}-F\tilde{\zeta}_{k-1}^{j}|0,S)$.
Since (15) is quadratic in $\zeta_{k}$, the maximizer $\zeta_{k}^{i+1}$ can be
found analytically,

$\displaystyle\zeta_{k}^{i+1}\text{ is such that
}\nabla_{\zeta_{k}}Q^{N}_{\zeta}(\zeta_{k}^{i+1},\zeta_{k}^{i})=0.$

Hence,

$\displaystyle\zeta_{k}^{i+1}$
$\displaystyle=\Big{[}H^{T}R^{-1}H+S^{-1}\Big{]}^{-1}\frac{1}{\sum_{l=1}^{N}\lambda^{l}\omega_{k-1|k-1}^{l}}\times$
$\displaystyle\Big{\\{}H^{T}R^{-1}y_{k}+S^{-1}F\sum_{j=1}^{N}\tilde{\zeta}_{k-1}^{j}\lambda^{j}\omega_{k-1|k-1}^{j}\Big{\\}},$
(16)

or, using the matrix inversion lemma,

$\displaystyle\zeta_{k}^{i+1}$
$\displaystyle=\Big{[}F\sum_{j=1}^{N}\tilde{\zeta}_{k-1}^{j}\lambda^{j}\omega_{k-1|k-1}^{j}+$
$\displaystyle\hskip-28.45274pt\mathcal{B}(y_{k}-HF\sum_{j=1}^{N}\tilde{\zeta}_{k-1}^{j}\lambda^{j}\omega_{k-1|k-1}^{j})\Big{]}\frac{1}{\sum_{l=1}^{N}\lambda^{l}\omega_{k-1|k-1}^{l}}$
(17)

with

$\displaystyle\mathcal{B}=SH^{T}(HSH^{T}+R)^{-1}.$

Note that, while (17) bears a passing resemblance to the Kalman filter
recursion as (16) does to the information filter, the $\lambda^{j}$ terms and
$\hat{Q}^{N}_{\zeta}(\zeta_{k},\zeta^{i}_{k})$ in (15) are dependent on the
previous iterate $\zeta^{i}_{k}$ and the update follows the maximizer.

Now take

$\displaystyle\zeta_{k}=\begin{bmatrix}\zeta_{k,1}\\\ \zeta_{k,2}\\\
\zeta_{k,3}\end{bmatrix},\;F=\begin{bmatrix}0.66&-1.31&-1.11\\\
0.07&0.73&-0.06\\\ 0.00&0.08&0.80\end{bmatrix},$ $\displaystyle
H=\begin{bmatrix}0&1&1\end{bmatrix},\,S=\begin{bmatrix}0.2&0&0\\\ 0&0.3&0\\\
0&0&0.5\end{bmatrix},\,R=0.1,$ $\displaystyle\zeta_{0|-1}=\begin{bmatrix}0\\\
0\\\ \end{bmatrix},\;\;\Sigma_{0|-1}=\begin{bmatrix}0.3&0&0\\\ 0&0.3&0\\\
0&0&0.3\end{bmatrix}.$

A simulation over $100$ time-steps, with $N=2000$ particles, is conducted with
the results plotted in Figure 1. The EMSF recursion at each time-step $k$, as
described by (17), is initialized by
$\zeta_{k}^{0}\sim\mathcal{N}(0,I_{3\times 3})$. The results show that the
EMSF converges to the Kalman Filter solution, which property follows from
Corollary 1.

Figure 1: The Kalman Filter solution in black, and the EMSF solution in orange
dots.

Figure 2, depicts the convergence of the iterates at three individual times of
$\zeta_{k,1}$, where the algorithm stops whenever the maximum absolute
relative error is below $0.005$, i.e.
$\max\,_{q\in\\{1,2,3\\}}|(\zeta_{k,q}^{i+1}-\zeta_{k,q}^{i})/\zeta_{k,q}^{i+1}|\leq
0.5\%$.

The EMSF algorithm is of $\mathcal{O}(N)$ computational complexity, per
iteration, with $N$ the number of particles. Ten simulations were conducted
using Matlab, an uncompiled interpretive program, on a computer with a 2.5 GHz
Intel Core i5-7200U processor and 8.00 GB of RAM. The running times for the
one hundred time points averaged over the ten runs were: 3 milliseconds for
the Kalman filter, 3.4 seconds for the bootstrap particle filter, and 26.4
seconds for EMSF including the bootstrap particle filter.

Figure 2: The convergence of the sequences $\\{\zeta_{k,1}^{i}:i\geq 0\\}$ for
$k=15,55,90$ in green, orange and purple, respectively.

###  3.5 Example 2: Nonlinear state-space model with multi-modal/skewed state
distribution

In this example, the EMSF is tested with a nonlinear state-space model with
highly skewed and/or multimodal state densities. The EMSF solution is compared
to the conditional mean of the bootstrap particle filter. This example also
serves as an introduction to the next section, in which the EMSF is shown to
reduce to a convergent fixed-point iteration for an important class of
nonlinear systems.

Consider the following nonlinear state-space model,

$\displaystyle\zeta_{k+1}$
$\displaystyle=\alpha_{k}\tanh{(\pi\zeta_{k})}+w_{k}$ (18) $\displaystyle
y_{k}$ $\displaystyle=\frac{1}{2}\zeta_{k}+v_{k}$ $\displaystyle w_{k}$
$\displaystyle\sim\mathcal{N}(0,\frac{1}{5}),\,v_{k}\sim\mathcal{N}(0,1)$
$\displaystyle\zeta_{0}$ $\displaystyle\sim\mathcal{N}(0,1)$

where $\alpha_{k}=(1+0.5\sin(2\pi k/20))$. The summation (8) in Algorithm 2,
after ignoring constant terms, can be simplified to,

$\displaystyle\hat{Q}^{N}_{\zeta}(\zeta_{k},\zeta_{k}^{i})$
$\displaystyle=\sum_{j=1}^{N}\Big{[}-(y_{k}-\frac{1}{2}\zeta_{k})^{2}-$
$\displaystyle\hskip-28.45274pt5(\zeta_{k}-\alpha_{k-1}\tanh(\pi\tilde{\zeta}_{k-1}^{j}))^{2}\Big{]}\times\lambda^{j}\omega_{k-1|k-1}^{j}.$
(19)

where
$\lambda^{j}=\mathcal{N}(\zeta_{k}^{i}-\alpha_{k-1}\tanh(\pi\tilde{\zeta}_{k-1}^{j})|0,1/5)$.
Since (19) is quadratic in $\zeta_{k}$, the unique maximizer is given by

$\displaystyle\zeta_{k}^{i+1}$
$\displaystyle=\frac{1}{10.5}\frac{1}{\sum_{l=1}^{N}\lambda^{l}\omega_{k-1|k-1}^{l}}\times$
$\displaystyle\sum_{l=1}^{N}\big{[}y_{k}+10\alpha_{k}\tanh{(\pi\tilde{\zeta}_{k-1}^{j})}\lambda^{j}\omega_{k-1|k-1}^{j}\big{]}.$
(20)

The simulation is carried out over $40$ time-steps, and $N=500$ particles. The
EMSF, at each time $k$, is defined by the recursion (20) and is initialized by
$5$ independent initial guesses $\zeta_{k}^{0}\sim\mathcal{U}[-2,2]$. The
number of iterations is set to be $\leq 10$, and the point of highest final
likelihood is taken over the five initial conditions. Notice that

$\displaystyle p(\zeta_{k}|\mathcal{Y}_{0:k})$
$\displaystyle=p(\zeta_{k}|y_{k},\mathcal{Y}_{0:k-1})$ $\displaystyle\propto
p(y_{k}|\zeta_{k}).p(\zeta_{k}|\mathcal{Y}_{0:k-1})$
$\displaystyle=p(y_{k}|\zeta_{k})\int
p(\zeta_{k}|\zeta_{k-1})p(\zeta_{k-1}|\mathcal{Y}_{0:k-1})\,d\zeta_{k-1}.$

By using the particle approximation of the filtered density
$p(\zeta_{k-1}|\mathcal{Y}_{0:k-1})$, we have

$\displaystyle p(\zeta_{k}|\mathcal{Y}_{0:k})$ $\displaystyle\approx
p(y_{k}|\zeta_{k})\sum_{j=1}^{N}p(\zeta_{k}|\tilde{\zeta}_{k-1}^{j})\omega_{k-1|k-1}^{j}$

which is the empirical filtered density, and can be used to compare the
likelihoods of the convergence points resulting from multiple initial guesses.
For this example,

$\displaystyle p(\zeta_{k}|\mathcal{Y}_{0:k})$
$\displaystyle\approx\mathcal{N}(y_{k}-\frac{1}{2}\zeta_{k}|0,1)\times$
$\displaystyle\sum_{j=1}^{N}\mathcal{N}(\zeta_{k}-\alpha_{k-1}\tanh(\pi\tilde{\zeta}_{k-1}^{j})|0,\frac{1}{5})\omega_{k-1|k-1}^{j}$

where the particles $\\{\tilde{\zeta}_{k-1}^{j}\\}$ and their weights
$\\{\omega_{k-1|k-1}^{j}\\}$ follow the definition as in the EMSF algorithm.
The results are shown in Figures 3 to 5.

Figure 3: The shading represents the empirical density function (up to a
positive multiplicative constant). In black is the output of the EMSF
algorithm, and in blue is the conditional mean of the particle filter. Figure
4: The bi-modal empirical density function of $\zeta_{k}$ at $k=10$ (apart
from a multiplicative constant). xEMSF is the solution of the EMSF algorithm
and xCM is the conditional mean of the bootstrap particle filter. Figure 5:
The highly skewed empirical density function of $\zeta_{k}$ at $k=25$. xEMSF
is the solution of the EMSF algorithm and xCM is the conditional mean of the
bootstrap particle filter.

In this example, the EMSF reduces to a fixed-point iteration. This leads to
the following result for state-space models with Gaussian noises and linear
measurement equations, in which this simplification of EMSF is possible.

###  3.6 Gaussian state-space models with linear measurement equations

In the case of a nonlinear state-space model with linear measurement equation
and affine Gaussian noises, Algorithm 2 reduces to a fixed-point iteration
method which is convergent to a stationary point of the empirical filtered
density.

######  Theorem 3.

For the nonlinear state-space system (1) defined in Section 2 with
$g(\zeta_{k})=H\zeta_{k}$ for $H\in\mathbb{R}^{r_{y}\times r_{\zeta}}$ and
$w_{k}\sim\mathcal{N}(0,S_{k})$, $v_{k}\sim\mathcal{N}(0,R_{k})$, Steps 2 and
3 of Algorithm 2 reduce to the following fixed-point iteration.

$\displaystyle\zeta_{k}^{i+1}$
$\displaystyle=\frac{1}{\sum_{l=1}^{N}\lambda^{l}\omega_{k-1|k-1}^{l}}\times\Big{[}\sum_{j=1}^{N}f(\tilde{\zeta}_{k-1}^{j})\lambda^{j}\omega_{k-1|k-1}^{j}+$
$\displaystyle\hskip
28.45274pt\mathcal{B}_{k}(y_{k}-H\sum_{j=1}^{N}f(\tilde{\zeta}_{k-1}^{j})\lambda^{j}\omega_{k-1|k-1}^{j})\Big{]},$

where

$\displaystyle\mathcal{B}_{k}=S_{k-1}H^{T}(HS_{k-1}H^{T}+R_{k})^{-1},$

and

$\displaystyle\lambda^{j}=\mathcal{N}(\zeta_{k}^{i}-f(\tilde{\zeta}_{k-1}^{j})|0,S_{k-1}),\quad
j=1,2,\ldots N.$

Proof. The summation (8) in Algorithm 2 reduces to

$\displaystyle\hat{Q}^{N}_{\zeta}(\zeta_{k},\zeta_{k}^{i})$
$\displaystyle=\sum_{j=1}^{N}\Big{[}-(y_{k}-g(\zeta_{k}))^{T}R_{k}^{-1}(y_{k}-g(\zeta_{k}))$
$\displaystyle\hskip-42.67912pt-(\zeta_{k}-f(\tilde{\zeta}_{k-1}^{j}))^{T}S_{k-1}^{-1}(\zeta_{k}-f(\tilde{\zeta}_{k-1}^{j}))\Big{]}\times\lambda^{j}\omega_{k-1|k-1}^{j}$

which is quadratic in $\zeta_{k}$. Hence, the unique maximizer of
$Q^{N}_{\zeta}(\zeta_{k},\zeta_{k}^{i})$ can be found by setting
$\nabla_{\zeta_{k}}Q^{N}_{\zeta}(\zeta_{k},\zeta_{k}^{i})$ to zero. Thus,

$\displaystyle\zeta_{k}^{i+1}$
$\displaystyle=\Big{[}H^{T}R_{k}^{-1}H+S_{k-1}^{-1}\Big{]}^{-1}\frac{1}{\sum_{l=1}^{N}\lambda^{l}\omega_{k-1|k-1}^{l}}\times$
$\displaystyle\Big{\\{}H^{T}R_{k}^{-1}y_{k}+S_{k-1}^{-1}\sum_{j=1}^{N}f(\tilde{\zeta}_{k-1}^{j})\lambda^{j}\omega_{k-1|k-1}^{j}\Big{\\}}.$

Applying the matrix inversion lemma, as in Section 3.4, completes the proof.
$\square$

We have the following ML state filtering algorithm.

Algorithm 3 Fixed-Point State Filter (FPSF) Algorithm for linear measurements
and affine Gaussian noises

1:Subject to the assumptions of Theorem 3 and given: $y_{k}$, the particles
$\\{\tilde{\zeta}^{j}_{k-1}\\}_{j=1}^{N}$, their normalized importance weights
$\\{\omega_{k-1|k-1}^{j}\\}_{j=1}^{N}$ and the initial guess $\zeta_{k}^{0}$,
set $i~{}\leftarrow~{}0$;

2:Do:  $\displaystyle\zeta_{k}^{i+1}$
$\displaystyle=\frac{1}{\sum_{l=1}^{N}\lambda^{l}\omega_{k-1|k-1}^{l}}\times\Big{[}\sum_{j=1}^{N}f(\tilde{\zeta}_{k-1}^{j})\lambda^{j}\omega_{k-1|k-1}^{j}+$
$\displaystyle\hskip
34.1433pt\mathcal{B}_{k}(y_{k}-H\sum_{j=1}^{N}f(\tilde{\zeta}_{k-1}^{j})\lambda^{j}\omega_{k-1|k-1}^{j})\Big{]}$
where  $\displaystyle\mathcal{B}_{k}=S_{k-1}H^{T}(HS_{k-1}H^{T}+R_{k})^{-1}$
and
$\displaystyle\lambda^{j}=\mathcal{N}(\zeta_{k}^{i}-f(\tilde{\zeta}_{k-1}^{j})|0,S_{k-1}),\quad
j=1,2,\ldots N$

3:If $\\{\zeta_{k}^{i}\\}_{i\geq 0}$ satisfies a convergence criterion,
terminate. Else: $i~{}\leftarrow~{}i+1$ and go to step 2.

##  4 ML state prediction by the expectation maximization algorithm

Given,

  * •

the problem statement and assumptions as in Section 2,

  * •

the particle approximation of the time-$m$ filtered density
$p(\zeta_{m}|\mathcal{Y}_{0:m})$, for $m\in\\{1,2,\ldots,k-1\\}$, i.e.

$p(\zeta_{m}|\mathcal{Y}_{0:m})\approx\sum_{j=1}^{N}\omega_{m|m}^{j}\delta(\zeta_{m}-\tilde{\zeta}^{j}_{m|m})$

  * •

initial state estimate $\zeta^{0}_{k}$,

the time-update of the bootstrap particle filter can be used to propagate the
particles and achieve the particle approximation

$\displaystyle p(\zeta_{k-1}|\mathcal{Y}_{0:m})$
$\displaystyle\approx\sum_{j=1}^{N}\omega_{m|m}^{j}\delta(\zeta_{k-1}-\tilde{\zeta}^{j}_{k-1|m}).$
(21)

The importance weights are unchanged because no extra measurements are
available after time-$m$. The EM algorithm can then be used to evaluate the
MLE of $\zeta_{k}$ of the predicted density $p(\zeta_{k}|\mathcal{Y}_{0:m})$,
using the particle approximation (21). The E-step is performed by evaluating
the following expectation

$\displaystyle Q^{p}_{\zeta}(\zeta_{k},\zeta_{k}^{i})$
$\displaystyle=\mathbb{E\,}_{\zeta_{k-1}|\zeta_{k}^{i},\mathcal{Y}_{0:m}}\\{\log
p(\zeta_{k}|\zeta_{k-1},\mathcal{Y}_{0:m})|\zeta_{k}^{i},\mathcal{Y}_{0:m}\\},$
$\displaystyle=\int\Big{\\{}\log p(\zeta_{k}|\zeta_{k-1},\mathcal{Y}_{0:m})$
$\displaystyle\hskip
28.45274ptp(\zeta_{k-1}|\zeta_{k}^{i},\mathcal{Y}_{0:m})\Big{\\}}\,d\zeta_{k-1},$

which, analogous to the proof of Theorem 1, is equivalent to

$\displaystyle Q^{p}_{\zeta}(\zeta_{k},\zeta_{k}^{i})$
$\displaystyle=\int\Big{\\{}\big{[}\log
W_{k-1}(\zeta_{k}-f(\zeta_{k-1}))\big{]}\times$
$\displaystyle\hskip-28.45274ptW_{k-1}(\zeta_{k}^{i}-f(\zeta_{k-1}))p(\zeta_{k-1}|\mathcal{Y}_{0:m})\Big{\\}}\,d\zeta_{k-1}.$
(22)

The particle approximation (21) can be used to construct $Q^{p,N}_{\zeta}$, an
approximate of (22) up to additive and multiplicative constants which are
independent from $\zeta_{k}$,

$\displaystyle\hat{Q}^{p.N}_{\zeta}(\zeta_{k},\zeta_{k}^{i})$
$\displaystyle=\sum_{j=1}^{N}\Big{[}\log
W_{k-1}(\zeta_{k}-f(\tilde{\zeta}_{k-1|m}^{j}))\Big{]}$ $\displaystyle\times
W_{k-1|m}(\zeta_{k}^{i}-f(\tilde{\zeta}_{k-1|m}^{j}))\omega_{m|m}^{j}.$

And the M-step is performed by maximizing this approximate over $\zeta_{k}$.
Hence, the EM algorithm for ML state prediction,

Algorithm 4 EM State Predictor (EMSP) Algorithm

1:Given the particles $\\{\tilde{\zeta}^{j}_{k-1|m}\\}_{j=1}^{N}$, their
normalized importance weights $\\{\omega_{m|m}^{j}\\}_{j=1}^{N}$ and
$\zeta_{k-1}^{*}$, set $i\leftarrow 0$ and let
$\zeta_{k}^{i}=f(\zeta_{k-1}^{*})+\mathbb{E}w_{k}$;

2:Evaluate the summation:
$\displaystyle\hat{Q}^{p.N}_{\zeta}(\zeta_{k},\zeta_{k}^{i})$
$\displaystyle=\sum_{j=1}^{N}\Big{[}\log
W_{k-1}(\zeta_{k}-f(\tilde{\zeta}_{k-1|m}^{j}))\Big{]}$ $\displaystyle\times
W_{k-1}(\zeta_{k}^{i}-f(\tilde{\zeta}_{k-1|m}^{j}))\omega_{m|m}^{j}$ (23)

3:Evaluate $\zeta_{k}^{i+1}$:
$\zeta_{k}^{i+1}=\text{arg}\max_{\zeta_{k}\in\mathbb{R}^{r_{\zeta}}}\hat{Q}^{p,N}_{\zeta}(\zeta_{k},\zeta_{k}^{i})$
(24)

4:If $\\{\zeta_{k}^{i}\\}_{i\geq 0}$ satisfies a convergence criterion: set
$\zeta_{k}^{*}=\zeta_{k}^{i}$ and terminate. Else: $i\leftarrow i+1$ and go to
step 2.

Similar versions of Theorem 2 and Corollary 1 hold for the sequence
$\\{\zeta_{k}^{i}\\}_{i\geq 0}$ generated by Algorithm 4, i.e. the sequence
$\\{\log p(\zeta_{k}^{i}|\mathcal{Y}_{m})\\}_{i\geq 0}$ is non-decreasing, and
if $\log p(\zeta_{k}|\mathcal{Y}_{m})$ is unimodal and its MLE is the only
stationary point, then $\\{\zeta_{k}^{i}\\}_{i\geq 0}$ of Algorithm 4
converges to the MLE.

##  5 ML state smoothing by the expectation maximization algorithm

In contrast to prediction and filtering, smoothed state estimation involves
considerably greater complexity and, in the MLE environment, greater demands
of the particle filter. This leads to an EM algorithm which has complexity
$\mathcal{O}(N^{2})$. The algorithm EMSS, however, follows directly, as we now
show.

Given,

  * •

the problem statement and assumptions as in Section 2,

  * •

the measurements sequence $\mathcal{Y}_{0:n}$, for $n>k+1$,

  * •

The particle approximations of the following densities:

$\displaystyle\hskip-14.22636ptp(\zeta_{k-1}|\mathcal{Y}_{0:k-1})$
$\displaystyle\approx\sum_{j=1}^{N}\omega_{k-1|k-1}^{j}\delta(\zeta_{k-1}-\tilde{\zeta}^{j}_{k-1|k-1}),$
(25) $\displaystyle p(\zeta_{k}|\mathcal{Y}_{0:k})$
$\displaystyle\approx\sum_{j=1}^{N}\omega_{k|k}^{j}\delta(\zeta_{k}-\tilde{\zeta}^{j}_{k|k}),$
(26) $\displaystyle p(\zeta_{k+1}|\mathcal{Y}_{0:n})$
$\displaystyle\approx\sum_{j=1}^{N}\omega_{k+1|n}^{j}\delta(\zeta_{k+1}-\tilde{\zeta}^{j}_{k+1|n}),$
(27)

which can be supplied by a forward–backward particle smoother [3],

  * •

initial state estimate $\zeta^{0}_{k}$,

The EM algorithm can be used to evaluate the MLE of $\zeta_{k}$ of the
smoothed density $p(\zeta_{k}|\mathcal{Y}_{0:n})$. The E-step is performed by
evaluating the expectation

$\displaystyle Q^{S}_{\zeta}(\zeta_{k},\zeta_{k}^{i})$
$\displaystyle=\mathbb{E\,}_{\zeta_{k+1}|\zeta_{k}^{i},\mathcal{Y}_{0:n}}\\{\log
p(\zeta_{k}|\zeta_{k+1},\mathcal{Y}_{0:n})|\zeta_{k}^{i},\mathcal{Y}_{0:n}\\},$
$\displaystyle=\int\log
p(\zeta_{k}|\zeta_{k+1},\mathcal{Y}_{0:n})p(\zeta_{k+1}|\zeta_{k}^{i},\mathcal{Y}_{0:n})\,d\zeta_{k+1}.$
(28)

Where in (28), using Bayes rule 3 times and Lemma 1, the density

$\displaystyle p(\zeta_{k}|\zeta_{k+1},\mathcal{Y}_{0:n})$
$\displaystyle=\frac{p(\mathcal{Y}_{k+1:n}|\zeta_{k+1},\zeta_{k},\mathcal{Y}_{0:k})p(\zeta_{k}|\zeta_{k+1},\mathcal{Y}_{0:k})}{p(\mathcal{Y}_{k+1:n}|\zeta_{k+1},\mathcal{Y}_{0:k})},$
$\displaystyle=\frac{p(\mathcal{Y}_{k+1:n}|\zeta_{k+1})p(\zeta_{k}|\zeta_{k+1},\mathcal{Y}_{0:k})}{p(\mathcal{Y}_{k+1:n}|\zeta_{k+1})},$
$\displaystyle=p(\zeta_{k}|\zeta_{k+1},\mathcal{Y}_{0:k}),$
$\displaystyle=\frac{p(\zeta_{k+1}|\zeta_{k})p(\zeta_{k}|\mathcal{Y}_{0:k})}{p(\zeta_{k+1}|\mathcal{Y}_{0:k})},$
$\displaystyle=\frac{p(\zeta_{k+1}|\zeta_{k})p(y_{k}|\zeta_{k})p(\zeta_{k}|\mathcal{Y}_{0:k-1})}{p(\zeta_{k+1}|\mathcal{Y}_{0:k})p(y_{k}|\mathcal{Y}_{0:k-1})}.$

The term $p(y_{k}|\mathcal{Y}_{0:k-1})$ is independent from $\zeta_{k}$ and
leads to an additive constant in the expectation (28). Thus,

$\displaystyle\log p(\zeta_{k}|\zeta_{k+1},\mathcal{Y}_{0:n})$
$\displaystyle\propto\log p(\zeta_{k+1}|\zeta_{k})+\log p(y_{k}|\zeta_{k})$
$\displaystyle+\log p(\zeta_{k}|\mathcal{Y}_{0:k-1}),$ $\displaystyle=\log
W_{k}(\zeta_{k+1}-f(\zeta_{k}))$ $\displaystyle\hskip-42.67912pt+\log
V_{k}(y_{k}-g(\zeta_{k}))+\log p(\zeta_{k}|\mathcal{Y}_{0:k-1}).$ (29)

The other density in (28), using the same conditional independence result and
Bayes rule,

$\displaystyle p(\zeta_{k+1}|\zeta_{k}^{i},\mathcal{Y}_{0:n})$
$\displaystyle=\frac{p(\zeta_{k}^{i}|\zeta_{k+1},\mathcal{Y}_{0:n})p(\zeta_{k+1}|\mathcal{Y}_{0:n})}{p(\zeta_{k}^{i}|\mathcal{Y}_{0:n})},$
$\displaystyle=\frac{p(\zeta_{k}^{i}|\zeta_{k+1},\mathcal{Y}_{0:k})p(\zeta_{k+1}|\mathcal{Y}_{0:n})}{p(\zeta_{k}^{i}|\mathcal{Y}_{0:n})},$
$\displaystyle\hskip-28.45274pt=\frac{p(\zeta_{k+1}|\zeta_{k}^{i})p(\zeta_{k}^{i}|\mathcal{Y}_{0:k})p(\zeta_{k+1}|\mathcal{Y}_{0:n})}{p(\zeta_{k+1}|\mathcal{Y}_{0:k})p(\zeta_{k}^{i}|\mathcal{Y}_{0:n})},$
(30)

where the terms $p(\zeta_{k}^{i}|\mathcal{Y}_{0:k})$ and
$p(\zeta_{k}^{i}|\mathcal{Y}_{0:n})$ are independent from $\zeta_{k+1}$ and
yield a positive multiplicative constant in the expectation (28). Hence,

$\displaystyle p(\zeta_{k+1}|\zeta_{k}^{i},\mathcal{Y}_{0:n})$
$\displaystyle\propto\frac{p(\zeta_{k+1}|\zeta_{k}^{i})p(\zeta_{k+1}|\mathcal{Y}_{0:n})}{p(\zeta_{k+1}|\mathcal{Y}_{0:k})},$
$\displaystyle\hskip-28.45274pt=\frac{W_{k}(\zeta_{k+1}-f(\zeta_{k}^{i}))p(\zeta_{k+1}|\mathcal{Y}_{0:n})}{p(\zeta_{k+1}|\mathcal{Y}_{0:k})}.$
(31)

Substituting (29) and (31) in (28), after dropping the additive and the
multiplicative constants, results in

$\displaystyle Q^{S}_{\zeta}(\zeta_{k},\zeta_{k}^{i})$ $\displaystyle=\log
V_{k}(y_{k}-g(\zeta_{k}))+\log p(\zeta_{k}|\mathcal{Y}_{0:k-1})+$
$\displaystyle\int\Big{\\{}\big{[}\log
W_{k}(\zeta_{k+1}-f(\zeta_{k}))\big{]}\times$
$\displaystyle\frac{W_{k}(\zeta_{k+1}-f(\zeta_{k}^{i}))p(\zeta_{k+1}|\mathcal{Y}_{0:n})}{p(\zeta_{k+1}|\mathcal{Y}_{0:k})}\Big{\\}}\,d\zeta_{k+1}.$
(32)

Notice that for the construction of (32), the densities
$p(\zeta_{k}|\mathcal{Y}_{0:k-1})$, $p(\zeta_{k+1}|\mathcal{Y}_{0:n})$ and
$p(\zeta_{k+1}|\mathcal{Y}_{0:k})$ are required. The particle approximation
(26) can be used to approximate $p(\zeta_{k+1}|\mathcal{Y}_{0:k})$,

$\displaystyle p(\zeta_{k+1}|\mathcal{Y}_{0:k})$ $\displaystyle=\int
p(\zeta_{k+1}|\zeta_{k})p(\zeta_{k}|\mathcal{Y}_{0:k})\,d\zeta_{k}$
$\displaystyle\approx\sum_{j=1}^{N}\omega_{k|k}^{j}W_{k}(\zeta_{k+1}-f(\tilde{\zeta}_{k|k}^{j})).$
(33)

Similarly, the density $p(\zeta_{k}|\mathcal{Y}_{0:k-1})$ can be approximated
using the particle approximation (25),

$\displaystyle\log p(\zeta_{k}|\mathcal{Y}_{0:k-1})$
$\displaystyle\approx\log\sum_{j=1}^{N}\omega_{k-1|k-1}^{j}W_{k-1}(\zeta_{k}-f(\tilde{\zeta}_{k-1|k-1}^{j})).$
(34)

Finally, by using the particle approximations (34) and (33) in (32), and (27)
to approximate the integral term in (32), after neglecting additive and
multiplicative constants, we get $\hat{Q}^{S.N}_{\zeta}$, the particle
approximation of $Q^{S}_{\zeta}$,

$\displaystyle\hat{Q}^{S.N}_{\zeta}(\zeta_{k},\zeta_{k}^{i})$
$\displaystyle=\log V_{k}(y_{k}-g(\zeta_{k}))+$
$\displaystyle\hskip-42.67912pt\log\sum_{j=1}^{N}\omega_{k-1|k-1}^{j}W_{k-1}(\zeta_{k}-f(\tilde{\zeta}_{k-1|k-1}^{j}))+$
$\displaystyle\hskip-42.67912pt\sum_{t=1}^{N}\Big{[}\log
W_{k}(\tilde{\zeta}_{k+1|n}^{t}-f(\zeta_{k}))\Big{]}\times$
$\displaystyle\hskip-42.67912pt\frac{W_{k}(\tilde{\zeta}_{k+1|n}^{t}-f(\zeta_{k}^{i}))}{\sum_{d=1}^{N}\omega_{k|k}^{d}W_{k}(\tilde{\zeta}_{k+1|n}^{t}-f(\tilde{\zeta}_{k|k}^{d}))}\omega_{k+1|n}^{t}.$
(35)

The EM algorithm for ML state smoothing,

Algorithm 5 EM State Smoother (EMSS) Algorithm

1:Given the sets of particles $\\{\tilde{\zeta}^{j}_{k-1|k-1}\\}_{j=1}^{N}$,
$\\{\tilde{\zeta}^{j}_{k|k}\\}_{j=1}^{N}$,
$\\{\tilde{\zeta}^{j}_{k+1|n}\\}_{j=1}^{N}$, and their corresponding
normalized importance weights, set $i~{}\leftarrow~{}0$ and let the initial
guess be $\zeta_{k}^{i}$;

2:Evaluate:  $\displaystyle\hat{Q}^{S.N}_{\zeta}(\zeta_{k},\zeta_{k}^{i})$
$\displaystyle=\log V_{k}(y_{k}-g(\zeta_{k}))+$
$\displaystyle\hskip-42.67912pt\log\sum_{j=1}^{N}\omega_{k-1|k-1}^{j}W_{k-1}(\zeta_{k}-f(\tilde{\zeta}_{k-1|k-1}^{j}))+$
$\displaystyle\hskip-42.67912pt\sum_{t=1}^{N}\Big{[}\log
W_{k}(\tilde{\zeta}_{k+1|n}^{t}-f(\zeta_{k}))\Big{]}\times$
$\displaystyle\hskip-42.67912pt\frac{W_{k}(\tilde{\zeta}_{k+1|n}^{t}-f(\zeta_{k}^{i}))}{\sum_{d=1}^{N}\omega_{k|k}^{d}W_{k}(\tilde{\zeta}_{k+1|n}^{t}-f(\tilde{\zeta}_{k|k}^{d}))}\omega_{k+1|n}^{t}.$
(36)

3:Evaluate $\zeta_{k}^{i+1}$:
$\zeta_{k}^{i+1}=\text{arg}\max_{\zeta_{k}\in\mathbb{R}^{r_{\zeta}}}\hat{Q}^{S,N}_{\zeta}(\zeta_{k},\zeta_{k}^{i})$
(37)

4:If $\\{\zeta_{k}^{i}\\}_{i\geq 0}$ satisfies a convergence criterion: set
$\zeta_{k}^{*}=\zeta_{k}^{i}$ and terminate. Else: $i\leftarrow i+1$ and go to
step 2.

##  6 Conclusion

This paper contributes a novel approach to incorporating the EM algorithm,
along with a particle filter, into the problem of maximum likelihood state
estimation. Algorithms for maximum likelihood state filtering, prediction, and
smoothing are presented, together with the desirable convergence properties
they inherit from the EM algorithm. Additionally, for the wide range of
nonlinear state-space systems possessing linear measurement equations and
affine Gaussian noises, these algorithms reduce to derivative-free
optimization through fixed-point iteration. The EMSF algorithms is tested on
two examples: one using a linear third-order Gaussian system for comparison
with the Kalman filter and to demonstrate; the concordance of the estimates,
the computational simplicity of the Kalman filter and the computational
feasibility of the EMSF combined with its attendant particle filter; the other
to estimate the MLE of a state density which exhibits multi-modality and
strong skewness. Further work is being pursued to understand the role of MLE
state estimates in constrained control, such as Model Predictive Control where
linearity and Gaussian assumptions can contravene the problem formulation.

## References

  * [1] David N DeJong, Roman Liesenfeld, Guilherme V Moura, Jean-François Richard, and Hariharan Dharmarajan.  Efficient likelihood evaluation of state-space representations.  Review of Economic Studies, 80(2):538–567, 2013.
  * [2] Arthur P Dempster, Nan M Laird, and Donald B Rubin.  Maximum likelihood from incomplete data via the EM algorithm.  Journal of the Royal Statistical Society: Series B (Methodological), 39(1):1–22, 1977.
  * [3] Arnaud Doucet, Simon Godsill, and Christophe Andrieu.  On sequential Monte Carlo sampling methods for Bayesian filtering.  Statistics and Computing, 10(3):197–208, 2000.
  * [4] Nikolas Kantas, Arnaud Doucet, Sumeetpal S Singh, Jan Maciejowski, Nicolas Chopin, et al.  On particle methods for parameter estimation in state-space models.  Statistical Science, 30(3):328–351, 2015.
  * [5] Kenneth Lange.  A gradient algorithm locally equivalent to the EM algorithm.  Journal of the Royal Statistical Society: Series B (Methodological), 57(2):425–437, 1995.
  * [6] Sylvain Le Corff and Gersende Fort.  Convergence of a particle-based approximation of the block online expectation maximization algorithm.  ACM Transactions on Modeling and Computer Simulation (TOMACS), 23(1):1–22, 2013.
  * [7] Erich L Lehmann and George Casella.  Theory of Point Estimation.  Springer Science & Business Media, 2006.
  * [8] Baibing Li, Cunjia Liu, and Wen-Hua Chen.  An auxiliary particle filtering algorithm with inequality constraints.  IEEE Transactions on Automatic Control, 62(9):4639–4646, 2016.
  * [9] Mengheng Li and Marcel Scharth.  Leverage, Asymmetry, and Heavy Tails in the High-Dimensional Factor Stochastic Volatility Model.  Journal of Business & Economic Statistics, pages 1–17, 2020.
  * [10] Andreas Lindholm and Fredrik Lindsten.  Learning dynamical systems with particle stochastic approximation EM.  arXiv preprint arXiv:1806.09548, 2018.
  * [11] Sheheryar Malik and Michael K Pitt.  Particle filters for continuous likelihood evaluation and maximisation.  Journal of Econometrics, 165(2):190–209, 2011.
  * [12] Xiao-Li Meng and Donald B Rubin.  Maximum likelihood estimation via the ECM algorithm: A general framework.  Biometrika, 80(2):267–278, 1993.
  * [13] Lyudmila Mihaylova, Avishy Y Carmi, François Septier, Amadou Gning, Sze Kim Pang, and Simon Godsill.  Overview of Bayesian sequential Monte Carlo methods for group and extended object tracking.  Digital Signal Processing, 25:1–16, 2014.
  * [14] Todd K Moon.  The expectation-maximization algorithm.  IEEE Signal Processing Magazine, 13(6):47–60, 1996.
  * [15] Michael K Pitt and Neil Shephard.  Filtering via simulation: Auxiliary particle filters.  Journal of the American Statistical Association, 94(446):590–599, 1999.
  * [16] George Poyiadjis, Arnaud Doucet, and Sumeetpal S Singh.  Particle approximations of the score and observed information matrix in state space models with application to parameter estimation.  Biometrika, 98(1):65–80, 2011.
  * [17] Herbert E Rauch, F Tung, and Charlotte T Striebel.  Maximum likelihood estimates of linear dynamic systems.  AIAA Journal, 3(8):1445–1450, 1965.
  * [18] Thomas B Schön, Adrian Wills, and Brett Ninness.  System identification of nonlinear state-space models.  Automatica, 47(1):39–49, 2011.
  * [19] François Septier, Sze Kim Pang, Simon Godsill, and Avishy Carmi.  Tracking of coordinated groups using marginalised MCMC-based particle algorithm.  In 2009 IEEE Aerospace Conference, pages 1–11. IEEE, 2009.
  * [20] Emilian R Vankov, Michele Guindani, Katherine B Ensor, et al.  Filtering and estimation for a class of stochastic volatility models with intractable likelihoods.  Bayesian Analysis, 14(1):29–52, 2019.
  * [21] Greg CG Wei and Martin A Tanner.  A Monte Carlo implementation of the EM algorithm and the poor man’s data augmentation algorithms.  Journal of the American Statistical Association, 85(411):699–704, 1990.
  * [22] CF Jeff Wu et al.  On the convergence properties of the EM algorithm.  The Annals of Statistics, 11(1):95–103, 1983.