\[ \newcommand{\Prob}[1]{\mathbb{P}\left( #1 \right)} \newcommand{\Expect}[1]{\mathbb{E}\left[ #1 \right]} \newcommand{\Var}[1]{\mathrm{Var}\left[ #1 \right]} \newcommand{\SampleVar}[1]{\widehat{\mathrm{Var}}\left[ #1 \right]} \newcommand{\Cov}[1]{\mathrm{Cov}\left[ #1 \right]} \newcommand{\TrueRegFunc}{\mu} \newcommand{\EstRegFunc}{\widehat{\TrueRegFunc}} \DeclareMathOperator{\tr}{tr} \DeclareMathOperator*{\argmin}{argmin} \DeclareMathOperator*{\argmax}{argmax} \DeclareMathOperator{\det}{det} \DeclareMathOperator{\dim}{dim} \newcommand{\TrueNoise}{\epsilon} \newcommand{\EstNoise}{\widehat{\TrueNoise}} \newcommand{\Signal}{S} \newcommand{\SignalNoise}{N} \newcommand{\AutoCov}{\gamma} \newcommand{\NoiseAutoCov}{\xi} \newcommand{\LinearExpect}[1]{\mathbb{L}\left[ #1 \right]} \newcommand{\Periodicity}{\tau} \newcommand{\Loadings}{\mathbf{b}} \newcommand{\FactorNoise}{\epsilon} \newcommand{\Uniqueness}{\mathbf{\psi}} \newcommand{\AltLoadings}{\mathbf{c}} \newcommand{\AltFactorNoise}{\zeta} \newcommand{\AltUniqueness}{\mathbf{\rho}} \newcommand{\DynMatrix}{\mathbf{a}} \newcommand{\DynNoise}{\eta} \newcommand{\DynVar}{\xi} \newcommand{\KalmanGain}{\mathbf{K}} \newcommand{\indep}{\perp} \newcommand{\Neighborhood}{\partial} \]

Previously

• There is a latent variable which evolves linearly in time (plus noise), and what we observe is linear in the latent (plus noise)
• Optimal linear prediction just is optimal prediction
• We could just use the Wiener/Kolmogorov/Krige technology with really big covariance matrices…
• … But we can do everything recursively (Kalman filter)

As you know

• Markov processes are ones where conditioning on the whole past is equivalent to conditioning on the current observation: \[ \Prob{X(t+1)|X(1:t)} = \Prob{X(t+1)|X(t)} \]
• This makes \(X(t)\) the state of the process at time \(t\), the variable which fixes the distribution of future events
• The Markov property makes it easy to:
  • Do probability calculations
  • Do prediction
  • Do inference
• The Markov property commits us to \(X(t+1) \indep X(1:t-1) | X(t)\)

Markov Property vs. Measurement Noise

• What if we don’t observe the state perfectly?
• \(S(t)\) is a Markov process
• We only observe \(X(t) = S(t) + \epsilon(t)\), where \(\epsilon(t)\) is IID noise
• Will \(X(t)\) also be Markov?

Why We Want the Probability of an Observation Sequence

\[ \Prob{X(1:n)=x(1:n)} = \sum_{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))}} \]

• When we predict, 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)}} \]
• \(\therefore\) Predicting comes down to find the probabilities of sequences and dividing
• \(\Prob{X(1:n)=x(1:n);\theta}=\) the likelihood of the data with parameter value \(\theta\) so we want this for inference as well
• But our equation is a sum over state histories, and that’s bad
  • Even 2 states means \(2^n\) possible values of \(s(1:n)\)
  • Never mind continuous states

Back to the Likelihood

• We’ve gotten predictions along the way with the forward algorithm
• The likelihood comes by multiplying: \[ \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)}} \]

Actually Doing the Calculations

• We can do everything exactly when \(S(t)\) and \(X(t)\) are both discrete
  • The original “forward algorithm”
• We can do everything exactly when \(S(t)\) and \(X(t)\) are both Gaussian, and everything’s linear
  • That’s the Kalman filter
• Otherwise we need numerical approximations

The Particle Filter: Code

particle.filter <- function(x, rinitial, rtransition, dobs, m) {
    n <- length(x)
    particles <- replicate(m, rinitial())
    particle.history <- matrix(0, nrow = m, ncol = n)
    for (t in 1:n) {
        weights <- dobs(state = particles, observation = x[t])
        particles <- sample(particles, size = m, replace = TRUE, prob = weights)
        particle.history[, t] <- particles
        particles <- sapply(particles, rtransition)
    }
    return(particle.history)
}

The Particle Filter: For Our Little Demo

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)

The Particle Filter: Run on Our Little Demo

Issues with the Particle Filter

1. Monte Carlo error: really want \(m\) to be as large as possible
   • But it costs computing time
   • Error is \(O(1/\sqrt{m})\) so diminishing returns
2. “Particle collapse”: if one \(w_i \gg w_j\) for all the other \(j\), then (almost) all new particles are copies of \(S^*_i\) and the approximation is bad
   • Lots of hacks around this
3. What if we don’t know transitions \(q\) or observation probabilities \(g\)?

Alternatives to the particle filter

• Extended Kalman filter: make a first-order Taylor expansion around your last estimate of the state, to get a (local) linear model, run that through the Kalman filter
  • Many variants with increasingly funny names (“unscented”, etc.)
• Laplace-Gauss Filter: Use Laplace approximation of the observation model to get a mean and variance, pretend that’s Gaussian, run it through the dynamic model (Koyama et al. 2010)

Parameter Estimation

• Be explicit about parameters in the notation: it’s \(q(s,s^{\prime};\theta)\) and \(g(s,x;\theta)\); how do we find \(\theta\)?
• We could use the likelihood: use forward algorithm to get \(\prod_{t}{\Prob{X(t)=x(t)|X(1:t-1) = x(1:t-1);\theta}}\), and maximize
  • Or, as usual, maximize the log
  • It’s a complicated function of \(\theta\) but it is just a function of \(\theta\)
• Alternate between estimate states and estimating parameters
  • Start with a guess about \(q\) and \(g\), estimate \(S(1:t)\)
  • Using \(S(1:t)\), estimate \(q\) and \(g\)
  • Iterate to convergence
• Calculating the likelihood is still computationally demanding, but simulating is easy \(\Rightarrow\) Use simulation-based-inference tricks (take 741)

What About Space?

• Hidden Markov random field: \(S(r)\) is a Markov random field, \(X(r) = f(S(r), \epsilon(r))\)
• \(\Rightarrow\) \(X(r)\) is not Markov
• Simulation in two stages:
  • First, simulate the \(S\) field
    • Usually by the Gibbs sampler (see backup)
  • Second, draw from \(X(r)|S(r)\) independently
• State estimation: usually by the “iterated conditional mode” algorithm of Besag (1986)
  • See back-up
• Likelihood: nothing as simple as the forward algorithm, often forced to do EM on the global state

What About Space-Time?

• Spatio-temporal hidden Markov random fields: \(S(r,t)\) is a spatio-temporal random field, \(X(r,t) = f(S(r,t), \epsilon(r,t))\)
• Proceed as for purely spatial random fields
• Or fall back on a dynamic factor model

State-Space Reconstruction / “Geometry from a Time Series”

• Suppose time-evolution and observation are both deterministic: \[\begin{eqnarray} X(t) & = & a(S(t))\\ S(t+1) &= & b(S(t)) \end{eqnarray}\] for some smooth functions \(a\), \(b\)
  • Continuous-time variants are very close
• Suppose \(S(t)\) lives on some \(d\)-dimensional manifold
• Consequence: \(X(t)\) sequences also live in some \(d\)-dimensional manifold: \[ (X(t), X(t+1), X(t+2), \ldots X(t+k)) = (b(S(t)), b(a(S(t))), \ldots (b \circ a^{k})(S(t))) \]
• \(\therefore\) If \(k\) is big enough, identical observed trajectories \(\Rightarrow\) identical latent states (Packard et al. 1980)

State-Space Reconstruction (cont’d.)

• \(S(t) \in [0,1]\), \(S(t+1) = 4 S(t) (1-S(t))\)
  • Logistic map
• \(X(t) = \cos{(0.188462 S(t) + 0.2437264)}\)
  • Randomly chosen frequency and phase shift

Vectors \((X(t), X(t+1), X(t+2))\), color-coded by the value of the initial condition \(S(t)\).

State-Space Reconstruction (cont’d.)

• Fact from differential topology (Whitney 1936):

Whitney embedding theorem: Suppose \(\dim{\mathcal{M}} = d\), and each coordinate of \(h: \mathcal{M} \mapsto \mathbb{R}^{k}\) is \(\mathcal{C}^1\). Then once \(k \geq 2d+1\), the typical \(h\) is a diffeomorphism, a differentiable mapping with a differentiable inverse

• (typical or generic is a somewhat complicated topological notion which we can pretty much gloss as “almost every” (Sauer, Yorke, and Casdagli 1991))

• Application to this problem (Takens 1981):

Takens embedding theorem: Suppose \(S(t)\) is on a \(d\)-dimensional manifold, we observe \(X(t) = b(S(t))\), and for the time-delay vector \(Y(t) = (X(t), X(t-\tau), \ldots X(t-(k-1)\tau))\). Then for typical choices of \(b\) and \(\tau\), \(Y(t)\) is diffeomorphic to \(S(t)\), once \(k \geq 2d+1\).

• Implication: Somehow find the right \(k\), and then do prediction in the time-delay-vector space, which is as good as we could observe the state space

State-Space Reconstruction (cont’d.)

• A single trajectory of the logistic map, from a random initial condition

State-Space Reconstruction (cont’d.)

• There is a whole alternative tradition of state-space modeling based on these theorems (Kantz and Schreiber 2004)
• Can actually tolerate moderate amounts of measurement noise, and use these ideas to do de-noising

Summing Up

• If we observe a Markov process with noise, the observables aren’t Markov
• We often postulate such hidden Markov models, because they let us fit non-Markov processes but still have some of the nice structure of Markov models
• Estimating the distribution of the hidden state \(S(t)\) given the history of observations \(X(1:t)\) is crucial
  • The forward algorithm lets us do this recursively
  • The particle filter lets us do this recursively by simulation
• Inference then proceeds either by likelihood (found using the forward algorithm) or maybe by simulation-based inference
  • Take 741
• The spatial and spatio-temporal variants are easy to formulate, but harder to work with

Many more details in the handout on the website

Backup: The Gibbs Sampler

• Assumes space \(r\) is discrete, but state \(X(r)\) can be continuous
• Start with some initial value of \(X(r)\) for all \(r\)
• Fix an order on the sites \(r\), and then:
  • \(X(r) \sim \Prob{X(r)|X(\Neighborhood(r))}\)
  • \(=\) draw a new value for \(X(r)\) from the conditional distribution / local characteristic and replace the old value of \(X(r)\)
  • N.B., later sites get conditioned on the updated value for \(X(r)\)
• Cycle through all sites at least once, and generally multiple times
• Each cycle is one step of a big Markov chain \(\Rightarrow\) converges on an invariant joint distribution; that invariant distribution is what we want
• See Kaplan et al. (2019) for a clever idea for speeding this up by simultaneous updating

Backup: Iterated Conditional Modes

• \(S(r)\) is a Markov random field, meaning:
  • \(\Neighborhood(r)\) are the sites adjacent to \(r\)
  • \(S(-r)\) is the field everywhere except at \(r\)
  • Markov property: \(\Prob{S(r)|S(-r)} = \Prob{S(r) | S(\Neighborhood(r))}\)
    • The neighbors screen off the rest of the field
• Notation: \(\Prob{S(r)=s|\Neighborhood(r) = q)} = p(s | q)\)
• \(X(r)\) is a function of \(S(r)\) and \(\TrueNoise(r)\), for IID noise \(\TrueNoise\)
• Notation: \(\Prob{X(r)=x|S(r)=s} = g(x|s)\)
• Use Markov to simplify inverse probability: \[ \Prob{S(r)|X, S(-r)} = \Prob{S(r)|X(r), S(\Neighborhood(r))} \]
• Use Bayes’s rule: \[ \Prob{S(r)=s|X(r)=x, S(\Neighborhood(r)) = q} \propto g(x|s) p(s|q) \]