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 Fraser (2008).

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 (1960) (and Kalman and Bucy (1961)). One of the clearest presentations I’ve seen is in Fraser (2008), 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 and Schmidt 1985), 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 (Koyama et al. 2010).

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².

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.

# 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)\).

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).

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.

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⁴ 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).

# 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.)

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)

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.

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:

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

## [1] 0.7069305

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