---
title: Hidden Markov Models and State Estimation
author: 36-467/36-667
date: 4 December 2018
bibliography: locusts.bib
---

\[
\newcommand{\Prob}[1]{\mathbb{P}\left( #1 \right)}
\]

```{r, include=FALSE}
set.seed(2018-12-3)
library(knitr)
opts_chunk$set(size="small", background="white",
               cache=TRUE, autodep=TRUE,
               tidy=FALSE, warning=FALSE, message=FALSE,
               echo=TRUE)
```


The last few lectures have focused on Markov processes, where conditioning on
the whole past is equivalent to conditioning on the most recent value:
\[
\Prob{X(t+1)|X(1:t)} = \Prob{X(t+1)|X(t)}
\]
When this is true, we say that $X(t)$ is the **state** of the process at time
$t$, the variable which determines the whole distribution of future
observations.  When this assumption holds, we can easily do likelihood-based
inference and prediction.  But the Markov property commits us to
$X(t+1)$ being independent of all earlier $X$'s given $X(t)$.

The most natural route from Markov models to hidden Markov models is to ask
what happens if we don't observe the state perfectly.  I have been using $X(t)$
to always stand for the time series we observe, so I am going to introduce a
new stochastic process, $S(t)$, which will be a Markov process, and say that
what we observe is $X(t)$ plus independent noise.  Will this $X(t)$ be Markov?

To be concrete, consider the following set-up.  $S(t)$ is a Markov chain with
two states, $-1$ and $+1$, and $\Prob{S(t+1)=+1|S(t)=+1} =
\Prob{S(t+1)=-1|S(t)=-1} q > 0.5$.  This chain thus alternates between two
states, but tends to stick in the same state for a while, to be "persistent".
What we observe is $X(t) = S(t) + \epsilon(t)$, where $\epsilon(t)$ are IID
Gaussian noise variables with mean 0 and variance $1$.  (Nothing particularly
turns on the choice of Gaussian noise or variance 1, etc.)  The figure below
illustrates both a typical realization of the observables $X(t)$, and of the
chain $S(t)$ that generated it.

```{r}
# From lecture 20: generic, abstract Markov simulator
rmarkov <- function(n, rinitial, rtransition) {
  x <- vector(length=n)
  x[1] <- rinitial()
  for (t in 2:n) {
    x[t] <- rtransition(x[t-1])
  }
  return(x)
}

# Simulate a chain that's +1/-1, with persistence in state
  # Pr(X(t+1)=+1|X(t)=+1) = Pr(X(t+1)=-1|X(t)=-1) = q, assumed > 0.5
  # Exercise: Prove that the invariant distribution is (1/2, 1/2) for all q
rpmchain <- function(n, q) {
    # Take one step of the chain
    transition.with.persistence <- function (x) {
        # Stay in the same state with probability q
        stay <- sample(c(TRUE, FALSE), size=1, prob=c(q, 1-q))
        if (stay) { return(x) }
        # Or flip to the opposite state
        if (!stay) { return(-x) }
    }        
    pm.seq <- rmarkov(n, rinitial=function() { sample(c(-1, +1),
                                                          size=1) },
                      rtransition=transition.with.persistence)
    return(pm.seq)
}

# Simulate a +1/-1 chain, but observe it through noise
noisy.rpmchain <- function(n, q, sd=1) {
    pm.seq <- rpmchain(n, q)
    x <- rnorm(n=n, mean=pm.seq, sd=sd)
    return(list(x=x, s=pm.seq))
}
```

```{r}
sim <- noisy.rpmchain(n=100, q=0.75)
plot(sim$x, xlab="t", ylab="X(t)", xlim=c(0, 120))
lines(1:100, sim$s, col="grey")
legend("right", legend=c("Hidden state S(t)", "Observable X(t)"),
       pch=c(NA, 1), lty=c("solid", "blank"), col=c("grey", "black"),
       cex=0.5)
```

Now any value of $X(t)$ _could_ have come from either value of $S(t)$.  (See
next figure.)  Large positive values of $X(t)$ are _really unlikely_ if
$S(t)=-1$, but not impossible, and values of $X(t)$ near zero are about equally
likely under either state.


```{r}
curve(0.5*dnorm(x,-1,1)+0.5*dnorm(x,1,1), xlab="x", ylab="p(x)",
      main="Marginal and Conditional Distributions of X(t)", xlim=c(-5,5), ylim=c(0, 0.5))
curve(dnorm(x,-1,1), lty="dashed", add=TRUE)
curve(dnorm(x,+1,1), lty="dotted", add=TRUE)
legend("topright", legend=c("X(t)", "X(t)|S(t)=-1", "X(t)|S(t)=+1"),
       lty=c("solid", "dashed", "dotted"))
```

This has consequences when we ask about the distribution of $X(t+1)$.  If $X(t)
= `r signif(max(sim$x), 2)`$, then it's very likely that $S(t)=+1$ --- so
likely that earlier history won't change much of anything.  But if $X(t) = 0$,
well, it will make a difference whether $X(t-1)$ was positive or negative.

This simple observation is enough, in itself, to tell us that $X(t)$ can't be a
Markov process.  In any Markov process, $X(t-1)$ is irrelevant to $X(t+1)$
given $X(t)$.  In particular, if $X$ is numerical and we regress $X(t+1)$ on
both $X(t)$ and $X(t-1)$, the true coefficient on $X(t-1)$ is _exactly_ zero.
This is _not_ the case here (Exercise 1).

Quite generally, if $S(t)$ is a Markov process, and $X(t)$ is a function of
$S(t)$, plus, possibly, some extra noise, then $X(t)$ is _not_
Markov[^functionofmarkov].

[^functionofmarkov]: This observation suggests the interesting mathematical question of whether _any_ functions of a Markov process are themselves Markovian, and, if so, which.  This has been studied by a variety of authors, often interested in "lumping" the states of large simulation models.  @Rosenblatt-on-markov-processes and @Rogers-Pitman-markov-functions are particularly notable contributions to this problem; @Nilsson-Jacobi-Gornerup-lumpability gives a necessary-and-sufficient algebraic condition.  On the related problem of characterizing the processes which are functions of Markov chains (whether or not the result is Markovian), see @Heller-functions-of-Markov-chains and @Erickson-functions-of-markov-chains.


# Definition of HMMs

Examples like these lead to a general notion of a **hidden Markov model**, or
**state-space model**.  In these models, there is a **latent** or **hidden**
state $S(t)$, which follows a Markov process.  We'll write
$\Prob{S(t+1)=r|S(t)=s} = q(r,s)$.  As in Markov models, the transitions need
to be complemented with a distribution for the initial state.  Because these
states are hidden, we don't actually get to observe them.  The **observable**
or **manifest** variable $X(t)$ is a random or noisy function of $S(t)$.  That
is, there's some sequence of IID noise variables, say $\epsilon(t)$, and $X(t)
= f(S(t), \epsilon(t))$ for some function $f$.  As a consequence, given $S(t)$,
$X(t)$ is independent of $X(t^\prime), S(t^\prime)$ for all $t^\prime \neq t$.
An _example_ of this would be when $X(t) = a S(t) + \epsilon(t)$, but the
function can be nonlinear and the noise can be non-additive.  In fact, the
exact function $f$ usually doesn't matter very much; what matters are the
conditional observation probabilities
$\Prob{X(t)=x|S(t)=s} = g(s,x)$.  

Some comments are in order on this definition.

- Both $S(t)$ and $X(t)$ can be numerical or qualitative variables, or indeed
more complicated.
    + Some people make a distinction between "hidden Markov models", where $S(t)$ is discrete, and "state-space models", where $S(t)$ is continuous.  Some of these people would further limit "state-space models" to situations where $X(t)$ is also continuous.  Since the theory and methods are almost the same, however, I will use the terms interchangeably.
    + Obviously, when variables are continuous, probability mass functions need to get replaced by probability density functions, and sums by integrals.
- Both $S(t)$ and $X(t)$ can have multiple dimensions.  I won't indicate this
in the notation.
- Any type of state variable can, potentially, be combined with any type of
observable.  For instance, you can have discrete $S(t)$ and
continuous-function-valued $X(t)$.  (These are sometimes used in speech
recognition, where the discrete states are the qualitative linguistic sounds,
"phonemes", people are saying, and the function-valued observations are the
actual noises they make over brief intervals of time.)
- HMMs can be extra _useful_ when the state is somehow simpler than the
  observable --- say a low-dimensional or even binary state explaining a
  high-dimensional vector of observables --- but that's not required by the
  definition.
     + In particular, if $S(t)$ has $r$ dimensions, but we observe $k < r$ of
	 them, the variables we _observe_ will not be Markov, even if they're
	 observed perfectly.

# Some Basic Calculations for HMMs

Because the states $S(t)$ form a Markov process, we can easily calculate
the probability of a state sequence:
\begin{eqnarray}
\Prob{S(1:n) = s(1:n)} & = & \Prob{S(1)=s(1)} \prod_{t=2}^{n}{S(t)=s(t)|S(1:t-1)=s(1:t-1)}\\
& = &  \Prob{S(1)=s(1)} \prod_{t=2}^{n}{S(t)=s(t)|S(t-1)=s(t-1)}\\
& = & \Prob{S(1)=s(1)} \prod_{t=2}^{n}{q(s(t-1), s(t))}
\end{eqnarray}
The first line is basic probability; the second uses the Markov property.

It's equally easy to calculate the probability of a sequence of observations,
given a sequence of states:
\begin{eqnarray}
\Prob{X(1:n)=x(1:n)|S(1:n)=s(1:n)} & = & \prod_{t=1}^{n}{\Prob{X(t)=x(t)|S(t)=s(t)}}\\
& = & \prod_{t=1}^{n}{g(s(t), x(t))}
\end{eqnarray}
This is because $X(t)$ is independent of everything else, given $S(t)$.

Putting these together gives us the joint probability of a sequence of states and a sequence of observations:
\begin{eqnarray}
\Prob{X(1:n)=x(1:n), S(1:n)=s(1:n)} & = & \Prob{S(1:n)=s(1:n)}\Prob{X(1:n)=x(1:n)|S(1:n)=s(1:n)}\\
& = & \Prob{S(1)=s(1)} \prod_{t=2}^{n}{q(s(t-1), s(t))}\prod_{t=1}^{n}{g(s(t), x(t))}
\end{eqnarray}

This, in turn, in principle gives us the probability of a sequence of observations:
\begin{eqnarray}
\Prob{X(1:n)=x(1:n)} & = & \sum_{s(1:n)}{\Prob{X(1:n)=x(1:n), S(1:n)=s(1:n)}}\\
\end{eqnarray}

This probability is a very valuable thing to know. When we make predictions, we
want
\[
\Prob{X(t)=x|X(1:t-1)=x(1:t-1)} = \frac{\Prob{X(1:t)=x(1:t)}}{\Prob{X(1:t-1)=x(1:t-1)}}
\]
so if we can calculate the probability of a sequence, and we can divide,
then we can predict.  If $q$ and $g$ contain some unknown parameters $\theta$,
and we want to do inference on them, then we really should be
writing things like
\[
\Prob{X(1:n)=x(1:n); \theta}
\]
and so calculating sequence probabilities will let us do likelihood-based
inference.

Unfortunately, the sum-over-state-histories expression, $
\sum_{s(1:n)}{\Prob{X(1:n)=x(1:n), S(1:n)=s(1:n)}}$, is just not suitable for
direct evaluation.  If there are even two states, there are $2^n$ possible
values of $s(1:n)$, and that's clearly not a sum you want to do if $n$ is at
all large.

# The Recursive Trick for Likelihood, Prediction, and State Estimation

#### a.k.a. "The Forward Algorithm"

There is a neat trick which lets us avoid the exponential explosion in
evaluating the likelihood, or in making predictions.  The same trick will also
let us estimate the hidden state at a given time.  The trick is to do all of
these things **recursively**, based on having solved the same problem with less
of the data.

Let's start with the base case.  If we knew $S(1)$, it'd be easy to find
the distribution of $X(1)$:
\[
\Prob{X(1)=x(1)|S(1)=s} = g(s, x(1))
\]
It's therefore easy to find the marginal distribution of $X(1)$, i.e.,
our predictive distribution for $X(1)$:
\[
\Prob{X(1)=x(1)} = \sum_{s}{\Prob{X(1)=x(1)|S(1)=s}\Prob{S(1)=s}} = \sum_{s}{g(s, x(1)) \Prob{S(1)=s}}
\]
Now we can use the definition of conditional probability to estimate the state at time 1:
\begin{eqnarray}
\Prob{S(1)=s|X(1)=x} & = & \frac{\Prob{S(1)=s, X(1)=x(1)}}{\Prob{X(1)=x(1)}}\\
& = & \frac{\Prob{S(1)=s}\Prob{X(1)=x(1)|S(1)=s}}{\Prob{X(1)=x(1)}}\\
& = & \frac{\Prob{S(1)=s} g(s, x(1))}{\sum_{s^{\prime}}{\Prob{S(1)=s} g(s, x(1))}}
\end{eqnarray}
This is, of course, just Bayes's rule.  Because this will be important
later, let's abbreviate it by $F_1(s)$.

Now that we know (the distribution of) the state at time 1, we can extrapolate
forward to time 2:
\begin{eqnarray}
\Prob{S(2)=r|X(1)=x(1)} & = & \sum_{s}{\Prob{S(2)=r, S(1)=s|X(1)=x(1)}}\\
& = & \sum_{s}{\Prob{S(2)=r|S(1)=s, X(1)=x(1)}\Prob{S(1)=s|X(1)=x(1)}}\\
& = & \sum_{s}{\Prob{S(2)=r|S(1)=s}\Prob{S(1)=s|X(1)=x(1)}}\\
& = & \sum_{s}{q(s,r) F_1(s)}
\end{eqnarray}
using the Markov property in the next-to-last line, and the definitions of $g$ and $F_1$ in the last line.

We want to extrapolate the state forward because it lets us make a prediction:
\begin{eqnarray}
\Prob{X(2)=x|X(1)=x(1)} & = & \sum_{s}{\Prob{X(2)=x, S(2)=s|X(1)=x(1)}}\\
& = & \sum_{s}{\Prob{X(2)=x|S(2)=s, X(1)=x(1)}\Prob{S(2)=s|X(1)=x(1)}}\\
& = & \sum_{s}{\Prob{X(2)=x|S(2)=s}\Prob{S(2)=s|X(1)=x(1)}}\\
& = & \sum_{s}{g(s, x)\Prob{S(2)=s|X(1)=x(1)}}
\end{eqnarray}
using the assumption that $S(t)$ screens off everything else from $X(t)$ in
the next-to-last line.


Finally, once we've seen $X(2)$, we can use that to get the distribution
of $S(2)$:
\begin{eqnarray}
F_2(s) & \equiv & \Prob{S(2)=s|X(2)=x(2), X(1)=x(1)}\\
& = & \frac{\Prob{S(2)=s, X(2)=x(2)|X(1)=x(1)}}{\Prob{X(2)=x(2)|X(1)=x(1)}}\\
& = & \frac{\Prob{X(2)=x(2)|S(2)=s, X(1)=x(1)}\Prob{S(2)=s|X(1)=x(1)}}{\Prob{X(2)=x(2)|X(1)=x(1)}}\\
& = & \frac{\Prob{X(2)=x(2)|S(2)=s}\Prob{S(2)=s|X(1)=x(1)}}{\Prob{X(2)=x(2)|X(1)=x(1)}}\\
& = & \frac{g(s, x(2))\Prob{S(2)=s|X(1)=x(1)}}{\Prob{X(2)=x(2)|X(1)=x(1)}}
\end{eqnarray}
Again, this is Bayes's rule.  Notice that the role of "prior distribution" here
isn't played by $\Prob{S(1)}$ or even $\Prob{S(1)|X(1)}$, but by
$\Prob{S(2)|X(1)}$.

Now we're ready to go on to $t=3$, and to state the general, recursive
pattern.  This is sometimes called the **forward algorithm**.

- Begin with $\Prob{S(t)=s|X(1:t)=x(1:t)} \equiv F_t(s)$
- Extrapolate the state forward:
\[
\Prob{S(t+1)=s|X(1:t)=x(1:t)} = \sum_{s^{\prime}}{q(s^{\prime}, s)\Prob{S(t)=s^{\prime}|X(1:t)=x(1:t)}} = \sum_{s^{\prime}}{q(s^{\prime}, s) F_t(s^{\prime})}
\]
- Predict the next observation:
\[
\Prob{X(t+1)=x|X(1:t)=x(1:t)} = \sum_{s}{g(s,x) \Prob{S(t+1)=s|X(1:t)=x(1:t)}}
\]
- Condition $S(t+1)$ on $X(t+1)$:
\[
\Prob{S(t+1)=s|X(1:t+1)=x(1:t+1)} = \frac{g(s,x) \Prob{S(t+1)=s|X(1:t)=x(1:t)}}{\Prob{X(t+1)=x|X(1:t)=x(1:t)}}
\]
- End with $F_{t+1}(s) \equiv \Prob{S(t+1)=s|X(1:t+1)=x(1:t+1)}$, or begin the cycle over again with $t+2$.

## Back to the Likelihood; Time Complexity of the Forward Algorithm

We began this by wanting an alternative to direct sum-over-state-histories
for the likelihood $\Prob{X(1:n)=x(1:n)}$.  We can get that from this
approach, by multiplying together predictive distributions:
\[
\Prob{X(1:n)=x(1:n)} = \Prob{X(1)=x(1)}\prod_{t=2}^{n}{\Prob{X(t)=x(t)|X(1:t-1)=x(1:t-1)}}
\]
We haven't, however, established that this approach is any faster than
summing over histories!

To see that it is, think about how time-consuming it will be to go through the
cycle above once.  To make a prediction for $X(t+1)$, we need to sum over
states $S(t+1)$, and to get the probability that $S(t+1)=s$, we need to sum
over states $S(t)$.  So if there are $k$ states, we need to so $O(k^2)$
operations in that double sum.  Each of the $n$ time-steps requires one pass
through the cycle, so we need up needing to do $O(k^2 n)$ operations, not
$O(k^n)$ as with direct summation.

In this analysis, it's important that the procedure is recursive, that we
calculate $F_{t+1}$ using just $F_t$ and $x(t+1)$ (and the model!), so that we
don't need to keep track of earlier states or their distributions, or even, in
fact, earlier values of $X$.  To the extent they matter, they've been
summarized in $F_t$, which is a sufficient statistic for predicting $X(t+1)$.

To sum up, in an HMM it's straightforward to calculate $\Prob{X(1)=x(1)}$ and
$F_1(s) = \Prob{S(1)=s|X(1)=x(1)}$.  Thereafter, we can predict by cycling
through (i) extrapolating the state distribution $F_t$ forward in time (using
the transition probabilitiess $q$), (ii) predicting the next observation (using
(i) and the observation probabilities $g$), and (iii) updating the distribution
of the state to $F_{t+1}$ (using (ii) and Bayes's rule).  This lets
us make predictions, and so calculate a likelihood, or log-likelihood.

## The Backward Algorithm

You may have noticed that we're just calculating $\Prob{S(t)|X(1:t)}$ and not
$\Prob{S(t)|X(1:n)}$, let alone $\Prob{S(1:n)|X(1:n)}$.  Partly, this is
because all we really _need_ for prediction or the likelihood is
$\Prob{S(t)|X(1:t)}$.  Partly, however, it is because it's just a lot _easier_
to calculate this than it is to calculate $\Prob{S(t)|X(1:n)}$.  Nonetheless,
there is _also_ a recursive algorithm for the latter.  To use it, one must
first have run the forward algorithm, getting a sequence $\Prob{S(t)|X(1:t)}$
for each $t\in 1:n$.  The last of these, $\Prob{S(n)|X(1:n)}$, is the last
probability we want, so the recursion starts at the _end_ of the time series
and works backwards in time (hence the name).  I will omit all the details
here, because they're neither as illuminating as in the forward algorithm,
and because you can find them in readily-available references, such
as @Andy-Fraser-on-HMMs.


## The Kalman Filter

An important special case is when $S(t)$ and $X(t)$ are both vectors,
and everything is linear-and-Gaussian:
\begin{eqnarray}
S(t+1) & = & \mathbf{a} S(t) + \eta(t)\\
X(t) & = & \mathbf{b} S(t) + \epsilon(t)\\
\eta(t) & \sim_{IID} & \mathcal{N}(0, \mathbf{\Sigma}_{S})\\
\epsilon(t) & \sim_{IID} & \mathcal{N}(0, \mathbf{\Sigma}_X)
\end{eqnarray}
In this case, it's possible to write out all the conditional probabilities in
the recursion _exactly_ --- everything is Gaussian, with mean vectors and
variance matrices that are straightforward functions of the data.  Finding
the filtering distribution then becomes known as the **Kalman filter**, after
@Kalman (and @Kalman-Bucy).  One of the clearest presentations I've seen
is in @Andy-Fraser-on-HMMs, which you should read if this is at all
interesting.  There is a base R implementation in the function `KalmanRun` 
(with various companion functions `KalmanLike`, `KalmanForecast`, etc.).

While the Kalman-filter solution is is nicely straightforward, it is heavily
tied to the linear-and-Gaussian assumptions.  There are many important
situations where that's good enough
[@McGee-Schmidt-discovery-of-kalman-filter], but it can fail in important
ways.  This has led to a number of attempts to use local linear approximations
to the Kalman filter in non-linear problems, and to alternative
closed-form approximate filters [@Big-Laplace-Gauss-filter-paper].

A different direction, however, is to go to back to the core ideas of the
forward algorithm, and to simulate them.  This is the "particle filter".


## The Particle Filter

While the forward algorithm is straightforward in principle, it still involves
a double sum over the state space.  This can be a time-consuming when the
number of states $k$ is large.  When the state-space is continuous, the sum
becomes an integral, and doing even a double integral numerically can be quite
time-consuming and/or inaccurate.  In this situation, a useful trick is to
replace complicated or time-consuming calculations with simulations,
i.e., to use Monte Carlo.

Applied here, the **Monte Carlo filter**, or **particle filter**, approximates
the distribution $F_t$ with the empirical distribution of a large
number of simulations of the HMM, conditioned on the data.  It works as follows[^particlefilternote].

[^particlefilternote]: See @Kunsch-on-particle-filters for a good, brief review of particle filters, and @Sequential-Monte-Carlo-book for a more exhaustive, though now a little old, treatment.  The actual history of who introduced which idea when is a bit more involved than I feel like going into, but covered in both these references.  @Douc-Moulines-Stoffer discusses statistical applications.

1. Make a large number of independent draws $S^*_1, S^*_2, \ldots S^*_m$ from
   the distribution $\Prob{S(1)}$.  Each random state value is called a
   **particle**, and the number of particles $m$ is, ideally, as large as
   practical.  (See more on this below.)  Set $t=1$.
2. For each particle $S^*_i$, calculate $g(S^*_i, x(t)) = w_{i}$.  Adding these
   up approximates the predictive probability,
   $\Prob{X(t)=x(t)|X(1:t-1)=x(1:t-1)} \approx \sum_{i=1}^{m}{w_i}$.
3. Now resample the particles, i.e., draw $m$ of them, with replacement, from
   $S^*$, with probabilities proportional to the $w_{i}$.  Call the new values
   $\tilde{S}_1, \tilde{S}_2, \ldots \tilde{S}_m$.  Set $\hat{F}_t$ to be the
   empirical distribution of the $\tilde{S}$.  That is, their empirical
   distribution approximates $\Prob{S(t)|X(1:t)=x(1:t)}$.
4. Increment $t$ by 1, and evolve the particles forward in time.  That is, for each $\tilde{S}_i$, pick
a $S^*_i$ according to $q(\tilde{S}_i, S^*_i)$.  The distribution of these new
$S^*$ approximates $\Prob{S(t+1)|X(1:t)=x(1:t)}$.
5. Go to (2).

This all lends itself wonderfully to code.

```{r}
# Particle filter for an HMM
# Inputs: time series of length n (x)
  # function generating draws for initial state (rinitial)
  # function generating draws for next state (rtransition)
  # function giving observational probabilities (dobs)
  # number of particles (m)
# Output: n*m array of particles, approximating Pr(S(t)|X(1:t))
# Presumes: states can be elements in a vector
  # observations can be elements in a vector
  # rinitial is a no-argument function (for compatibility with rmarkov())
  # rtransition is a one-argument function (ditto)
  # dobs is a two-argument function
particle.filter <- function(x, rinitial, rtransition, dobs, m) {
    # How long is the time series?
    n <- length(x)
    # Our initial guesses at the state, before seeing anything
    particles <- replicate(m, rinitial())
    # A storage area for our guesses
    particle.history <- matrix(0, nrow=m, ncol=n)
    for (t in 1:n) {
        # Calculate observationa likelihoods/weights for current observation
        # given state=particles
        weights <- dobs(state=particles, observation=x[t])
        # Resample the particles, with probabilities (proportional to)
        # the weights
        particles <- sample(particles, size=m, replace=TRUE, prob=weights)
        # Store the current set of guesses
        particle.history[,t] <- particles
        # Evolve the particles forward using the Markov transitions for the
        # next time step
        particles <-sapply(particles, rtransition)
    }
    return(particle.history)
}
```
I will leave it as an exercise to extend this function so that it also
returns (an approximation to) the sequence $\Prob{X(t)|X(1:t-1)}$.

Here is an implementation using the HMM introduced at the beginning,
where $S(t)$ is either $+1$ or $-1$, and $X(t)|S(t) \sim \mathcal{N}(S(t), 1)$.

```{r}
sim.pf <- particle.filter(x=sim$x,
                          rinitial=function() { sample(-1,+1, size=1) },
                          rtransition=function (x) {
                              stay <- sample(c(TRUE, FALSE), size=1, prob=c(0.75, 0.25))
                              if (stay) { return(x) }
                              if (!stay) { return(-x) }
                          },
                          dobs=function(state, observation) {
                              dnorm(x=observation, mean=state, sd=1)
                          },
                          m=100)
```

We have 100 particles, i.e., 100 guesses for $S(t)$ at each time point.  Rather
than try to visualize this directly ---- it's rather "busy" --- I'll summarize
with lines at the 5% and 95% percentiles (so a 90% probability interval).

```{r}
plot(sim$x, xlab="t", ylab="X(t)", xlim=c(0, 120))
lines(1:100, sim$s, col="grey", lwd=2)
pf.limits <- apply(sim.pf, 2, quantile, prob=c(0.05, 0.95))
lines(1:100, pf.limits[1,], col="blue")
lines(1:100, pf.limits[2,], col="blue")
legend("right", legend=c("Hidden state S(t)", "Observable X(t)", "Particle filter"),
       pch=c(NA, 1, NA), lty=c("solid", "blank", "solid"), col=c("grey", "black", "blue"),
       cex=0.5)
```

Notice that there are places where the limits become very narrow --- these
are mostly situations where $X(t)$ has a large absolute value, and do it's
very easy to tell what $S(t)$ must be.

### An Analogy for the Particle Filter

A common analogy for understanding how the particle filter works is to think of
the particles as an evolving population, with the state of each particle being
its genotype.  The probability $w_{i}$ that particle $i$ gives to the current
observation is then the "fitness" of particle $i$.  In the resampling step,
each particle has a chance to reproduce, and the expected number of offspring
of particle $i$ is $\propto w_i$.  Finally, moving the new particles according
to the Markov chain is like introducing mutations.  So the whole particle
filter process is analogous to "breeding" a population with the goal of giving
high likelihood (=fitness) to the data.[^bayesasevol]

[^bayesasevol]: In one direction, you can extend this analogy to techniques of "evolutionary optimization", such as "genetic algorithms", which try to breed good solutions to complex optimization problems [@Holland-adaptation;@MM-on-GAs].  On the other, you can extend it to all forms of Bayesian inference, which are mathematically equivalent to natural selection without any sort of variation [@CRS-dynamics-of-bayes].


### Issues with the Particle Filter

1. Using a finite number of particles (= Monte Carlo samples) introduces
   approximation error --- the distribution $\hat{F}_t$ isn't quite the same as
   $F_t$.  This error shrinks as the number of particles $m \rightarrow
   \infty$.  All else being equal, $m$ should always be as large as possible.
   + What might not be equal is that increasing $m$ costs more time to run.
     Each step of the particle filter --- moving each particle according to
     the Markov process, calculating a probability for $X(t)=x(t)$ given the
     particle state, and re-sampling the particles --- takes $O(m)$ operations.
     (You can parallelize the transitions and the likelihood calculations, but
     not the resampling.)
   + There are diminishing returns to increasing $m$ --- the approximation
     error is generally $O(1/\sqrt{m})$.
   + Exactly how large you need to make $m$ thus depends on how you
     value computing time vs. numerical accuracy, which I can't tell you.
2. The particle filter will perform badly if at some time $t$, one particle $i$
   gives much higher probability to $X(t)$ than all the others, since then $w_i
   \gg w_j$, and all (or almost all) of the resampled particles will just be
   copies of $w_i$.  This "particle collapse" is most apt to be a problem when
   the data does something which is very improbable according to the model,
   either because of an actual outlier, o because of model mis-specification.
   @Bengtsson-Bickel-Li-on-particle-filters studies this problem in some detail.
   A slightly hacky solution is to add a little bit of noise on to the resampled
   particles.  (This is equivalent to smoothing the empirical distribution.)
3. I have presented the particle filter assuming the state-transition
   probabilities $q$ and the state-observation probabilities $g$ are _known_.
   Obviously, this is often not the case for real data problems, which brings
   us to the topic of parameter estimation.

# Parameter Estimation

Up to now, I have been taking the state-transition probabilities
$q(s,s^{\prime})$ and the state-observation probabilities $g(s,x)$ as known
functions.  This is obviously unrealistic.  They will generally depend on
one or more unknown parameters $\theta$, so I should really be writing
$q(s,s^{\prime};\theta)$ and $g(s,x;\theta)$.  How might we estimate $\theta$?

- _Use the likelihood_ We can use the forward algorithm to calculate, or at least approximate, $\Prob{X(1:n)=x(1:n);\theta}$ for each $\theta$, as
$\prod_{t}{\Prob{X(t)=x(t)|X(1:t-1) = x(1:t-1);\theta}}$.  This is a function
of $\theta$, which we can maximize.  In the usual way, it will be easier,
and more numerically stable, to work with the log-likelihood,
\[
\sum_{t}{\log{\Prob{X(t)=x(t)|X(1:t-1)=x(1:t-1);\theta}}}
\]
A fully rigorous proof that maximum likelihood estimation is consistent for
general HMMs is surprisingly tricky, but doable
[@Cappe-Moulines-Ryden-inference-in-HMMs].  As a practical matter, however,
the log-likelihood is a function, and can be optimized like anything else.
    + One issue with the likelihood-based approach is that if either $X(t)$ or
	$S(t)$ is a continuous variable, we _must_ impose some parametric structure
	on $g$ and/or $q$, since otherwise maximizing the likelihood leads to bad
	answers.  Checking these assumptions is important, and potentially tricky.
- _Alternate between guessing at states and estimation_ If we _knew_ the
   sequence of states $S(t)$, it would be straightforward to estimate both $q$
   and $g$.  We could estimate $q$ using any technique we'd use to estimate
   transition probabilities for a Markov process, and we'd estimate $g$
   using any technique for conditional distributions.  (We might even use
   regression in both problems, if we could believe in additive noise.)
   In the other direction, if we knew $q$ and $g$, we can estimate $S(t)$
   for each $t$.  This suggests a whole slew of iterative-approximation
   methods, which start _either_ with a guess for the $S$'s or for $q$ and $g$,
   and then alternate between using the current estimate for one to improve
   the estimate of the other.
   + The most famous iterative-approximation approach is the
     **Baum-Welch**, **forward-backward**, **expectation-maximization**, or
     **EM** algorithm.  It uses the backward algorithm (see above) to find
     $\Prob{S(t)|X(1:n)}$, using the current estimates for $q$ and $g$.  (This
	 is the "E" step, since a probability is an expectation.)  It then does
	 weighted maximum likelihood to re-estimate $q$ and $g$, with every
	 possible sequence of states getting a weight given by its probability
	 under the backward algorithm.  (This is the "M" step, since it's
	 maximizing the likelihood.)  It then re-runs the backward algorithm to
	 re-estimate the states, and so forth to convergence.  One can
	 show that both the E-step and the M-step increase the log-likelihood, so
	 this has to converge on at least a local maximum, though perhaps not
	 a global one.  (See @CRS-ADAfaEPoV for a tutorial introduction to EM.)
   + Variants include replacing finding the whole distribution
     $\Prob{S(t)|X(1:n)}$ with drawing samples from it (perhaps with the
	 particle filter), or even just finding the most probable sequence
	 $\hat{s}(1:n)$.
   + Non-parametric variants of the iterative-approximation idea are
     straightforward to code up, but hard to prove things about.
 - _Simulation-based inference_ Even with all the tricks for the forward
   algorithm, filtering, etc., calculating the likelihood can be difficult
   for complex HMMs.  These can, however, be very straightforward to simulate.
   All of the techniques for simulation-based inference from earlier lectures
   can be brought to bear here, especially indirect inference.

### Example of Simulation-Based Inference

Let's work through an example of simulation-based inference, with our toy
problem.  (The advantage of trying it out on a toy is that we can tell whether
or not it's working.)  Specifically, we'll use indirect inference, with an
autoregressive auxiliary model.  We'll need (i) a function to simulate the
generative model, (ii) a function to estimate the auxiliary model,
and (iii) a function which optimizes the generative model so that auxiliaries
fit to its simulation match auxiliaries fit to the data.

We wrote (i) already, `noisy.rpmchain`.  We need to write (ii) and (iii).

For (ii), the estimator of the auxiliary model, we'll use an autoregressive
model.  There is one[^moreparams] parameter we need to estimate, the probability
$\Prob{S(t+1)=+1|S(t)=+1} = \Prob{S(t+1)=-1|S(t)=-1} = q$.  This suggests
an AR(1) should be enough.  To give this a bit more information, though,
we'll fit an AR(2).

[^moreparams]: If we wanted to, we could treat all the entries in the transition matrix for $S$ as unknown, and estimate them, in the same way we're doing here, but this is just a demo.  Similarly, we could treat the variance of the observational noise as unknown and estimate it, etc.

```{r}
# Estimate the auxiliary model (here, an AR(2)) for indirect inference
# Inputs: a time series (x)
# Output: Point estimate of the auxiiary model (here, a vector of length 2)
aux.estimator <- function(x) {
    ar2.fit <- ar.ols(x, aic=FALSE, order=2)
    return(ar2.fit$ar)
}

# Repeatedly simulate the generative model, fit the auxiliary, and average
# Inputs: parameter of generative model (q)
  # number of simulations to run (s)
  # length of each simulation (n)
# Output: vector of averaged auxiliary parameter estimates
aux.from.sim <- function(q, s, n=length(sim$x)) {
    estimates <- replicate(s, aux.estimator(noisy.rpmchain(n=n, q=q)$x))
    return(rowMeans(estimates))
}    
```

We should check that there's a unique mapping from the parameter(s) in the
generative model to the auxiliary parameters (so that inverting this mapping in
principle lets us recover the generative parameter).  We can do this here by
plotting how the auxiliary parameters vary as we sweep over the generative
parameter.  (Obviously this sort of visual check gets harder with more
parameters on each side.)

```{r}
plottable.aux.from.sim <- Vectorize(aux.from.sim, vectorize.args="q")
ar2.vs.q <- sapply(seq(from=0.01, to=0.99, by=0.01),
                   plottable.aux.from.sim, s=500)
```

```{r}
plot(ar2.vs.q[1,], ar2.vs.q[2,], type="l", xlab="Coefficient on X(t-1)",
     ylab="Coefficient on X(t-2)")
text(ar2.vs.q[1,1], ar2.vs.q[2,1], labels="q=0.01")
text(ar2.vs.q[1,99], ar2.vs.q[2,99], labels="q=0.99")
```




We can now follow the examples in earlier lectures to write out the
indirect-inference function.  One thing worth noting here is that we know the
parameter $q$ has to be between $0$ and $1$ (it's a probability), and,
fortunately, `optim` has a good method for optimizing a single parameter
between known limits.  With multiple parameters, it might be worthwhile to
use a package like `alabama`, which allows for very general constraints on
the parameters.

```{r}
toy.indirect.inference <- function(x, starting.q, s) {
    n <- length(x)
    # What's the auxiliary model on the actual data?
    aux.data <- aux.estimator(x)
    # Calculate the discrepancy between the auxiliary on the data and
    # the auxiliary we get from simulating at a particular parameter value
    discrepency <- function(param) {
        q <- param
        # Get the average of auxiliary estimates over s simulation runs
        aux.mean <- aux.from.sim(q=q, s=s, n=n)
        # Mean squared difference over auxiliary parameters
        return(sum((aux.data-aux.mean)^2))
    }
    # Now minimize the discrepancy
    fit <- optim(par=starting.q, fn=discrepency, method="Brent", lower=0,
                 upper=1)
    # Return just the best-fitting value
    return(fit$par)
}    
```

Now let's actually run it:

```{r}
(toy.ii <- toy.indirect.inference(x=sim$x, starting.q=0.5, s=500))
```

Recall that the true value of $q$ was $0.75$, so this is doing pretty well.  If
we want to get a standard error, we could either follow the theory in Lecture
19, or just do a bootstrap.  (Since we know this is a stationary process, that
could either be a parametric or a non-stationary block bootstrap.)

Note that now that we could estimate the parameter of the generative
model without also _having_ to estimate the unobserved state.  Having estimated
$q$, however, we can now estimate the states $S(t)$, say with the particle
filter.

## Some Important Special Cases

### Mixture Models

A sequence of IID variables is, technically, a Markov process.
So if the sequence $S(1), \ldots S(n)$ is IID, we _could_ apply the techniques
above to estimate the latent variables $S(t)$.  But most of the work would be
unnecessary.  Because $S(t)$ is independent of both $S(t-1)$ and $S(t+1)$,
it breaks up into a collection of independent problems, of working out
the distribution of $S(t)$ given $X(t)$ (and _only_ given $X(t)$).
When the latent variable $S(t)$ is discrete, we call this problem
a "mixture model" or a "latent class model".  When $S(t)$ is continuous, the
problem has multiple names; "factor model" is common when $X(t)$ is linearly
related to $S(t)$.

@CRS-ADAfaEPoV goes over these situations in some detail.

### Deterministic Latent Dynamics

If there is actually a _deterministic_ relationship between $S(t)$ and $S(t+1)$,
say $S(t+1) = \phi(S(t))$ for some function $\phi$, that, too, is technically
a Markov process.  What makes this case special is that giving the initial
value $S(1)$ fixes the whole rest of the $S(t)$ sequence.  Thus there isn't
the exponential growth in the number of possible trajectories we see with
other Markov models.  In fact, it makes sense to treat the initial state
$S(1)$ as an unknown but _fixed_ parameter, say $s(1)$.  The likelihood is then
\[
\prod_{t=1}^{n}{g(s(t), x(t); \theta)}
\]
with the understanding that $s(t) = \phi^t(s(1))$.  Of course, if the dynamics
in the function $\phi$ are very sensitive to small changes in the initial state,
optimizing this likelihood can be very difficult, and so it might make more
sense to still fall back on techniques like simulation-based inference.


   



# Exercises

To think through / practice on, not to hand in.

1. Think back to the first example: $S(t)$ is a Markov chain with states $+1$
   and $-1$, $\Prob{S(t+1)=+1|S(t)=+1} = \Prob{S(t+1)=-1|S(t)=-1}=q$, and $X(t)
   \sim \mathcal{N}(S(t), 1)$.  Assume that $S(1)$ is drawn from the invariant
   distribution, which is $(1/2, 1/2)$, so that the process is stationary.
   Find (a) the variance of $X(t)$; (b) all the covariances between
   $X(t-1)$, $X(t)$ and $X(t+1)$, and (c) the optimal coefficients for
   a linear prediction of $X(t+1)$ from $X(t)$ and $X(t-1)$.  (You will need
   to invert a $2\times 2$ matrix.)


# References
