--- title: Nonlinear State-Space Models date: 22 September 2026 (Lecture 9) author: 36-740 output: slidy_presentation: math_method: engine: mathjax url: https://bactra.org/mathjax/tex-svg.js bibliography: locusts.bib --- ```{r, include=FALSE} set.seed(2026-09-22) # General set-up options library(knitr) opts_chunk$set(size="small", background="white", highlight=FALSE, cache=TRUE, autodep=TRUE, tidy=TRUE, tidy.opts=list(comment=FALSE), warning=FALSE, message=FALSE, echo=FALSE) ``` \[ \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} \] ```{r, include=FALSE} library(knitr) opts_chunk$set(size="small", background="white", cache=TRUE, autodep=TRUE, tidy=TRUE, tidy.opts=list(comment=FALSE), warning=FALSE, message=FALSE, echo=FALSE) ``` ## 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? ## A Concrete Example {.smaller} - $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} \equiv q > 0.5$ - We observe not $S(t)$ but $X(t) = S(t) + \epsilon(t)$ - Noise $\epsilon(t) \sim \mathcal{N}(0,1)$ and is IID - Is this $X(t)$ a Markov process? ```{r} # From lecture 18: 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) ``` ## The Current State "Screens Off" the Past from the Future {.smaller} $\Prob{X(t+1)=x(t+1)|S(t)=s, X(1:t)=x(1:t)}$ \begin{eqnarray} & = & \sum_{r}{\Prob{X(t+1)=x(t+1)|S(t+1)=r, S(t)=s, X(1:t)=x(1:t)}}\\ \nonumber & & {\Prob{S(t+1)=r|S(t)=s,X(1:t)=x(1:t)}}\\ & = & \sum_{r}{\Prob{X(t+1)|S(t+1)=r}\Prob{S(t+1)=r|S(t)=s}} \end{eqnarray} ## Relating $X(t)$ to Its Past through the State {.smaller} $\Prob{X(t+1)=x(t+1)|X(1:t)=x(1:t)}$ \begin{eqnarray} & = & \sum_{s}{\Prob{X(t+1)=x(t+1)|S(t)=s,X(1:t)=x(1:t)}\Prob{S(t)=s|X(1:t)=x(1:t)}}\\ & = & \sum_{s}{\left(\sum_{r}{\Prob{X(t+1)|S(t+1)=r}\Prob{S(t+1)=r|S(t)=s}}\right) \Prob{S(t)=s|X(1:t)=x(1:t)}} \end{eqnarray} $\Rightarrow$ The only way $X(t-h)$ can matter for $X(t+1)$ is by changing the conditional distribution of $S(t)$ ## Relating $X(t)$ to $S(t)$ {.smaller} - Any value of $X(t)$ _could_ have come from _either_ value of $S(t)$ ```{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), col="purple", lty="dotted") curve(dnorm(x,-1,1), add=TRUE, col="blue") curve(dnorm(x,+1,1), add=TRUE, col="red") legend("topright", legend=c("X(t)", "X(t)|S(t)=-1", "X(t)|S(t)=+1"), lty=c("dotted", "solid", "solid"), col=c("purple", "blue", "red")) ``` - $X(t) \gg 0$ is _unlikely_ if $S(t)=-1$ but not impossible - $X(t) \approx 0$ is about equally like for either value of $S(t)$ ## History Matters for $X(t)$ {.smaller} - If $X(t) = `r signif(max(sim$x), 2)`$ then it's really likely that $S(t)=+1$ and conditioning on earlier $X$s can't really change that - If $X(t) = 0.1$ then $S(t)$ is almost equally likely to be $+1$ or $-1$ - But not if $X(t-1)=-1$ or $X(t-1)=1.2$ - $\therefore$ Earlier values of $X$ give information about $S(t)$ when $X(t)$ isn't too extreme, _and so_ they give information about $X(t+1)$ - In a Markov process, $X(t-1) \indep X(t+1)|X(t)$ - That isn't true here - $\therefore$ the observed process isn't Markov ## Hidden Markov Models / State-Space Models {.smaller} - **Hidden** or **latent** (= unobserved) **state** $S(t)$, which follows a Markov process - Transition probabilities: \[ q(s,r) \equiv \Prob{S(t+1)=r|S(t)=s} \] + Need a distribution for initial state as well + A.k.a. transition model, process model, state model, plant model - **Observable** or **manifest** $X(t)$ is a noisy function of the current state: \[ X(t) = f(S(t), \epsilon(t)) \] for some function $f$ and IID noise $\epsilon$ + Function $f$ needn't be additive and might be very nonlinear - Leads to **observation probabilities** or **observation model** \[ g(s,x) \equiv \Prob{X(t)=x|S(t)=s} \] (for continuous variables, replace PMFs by pdfs as needed) ## Some Basic Calculations for HMMs {.smaller} - Probability of a state sequence from Markov property: \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} - Probability of an observed sequence, given a state sequence: \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} - Joint probability of states and 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} - Probability of an observation sequence: \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)}}\\ & = & \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))}} \end{eqnarray} ## 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 ## A Recursion Trick (a.k.a. "The Forward Algorithm") [1] {.smaller} - If we knew $S(1)$, it'd be easy to predict $X(1)$: \[ \Prob{X(1)=x|S(1)=s} = g(s, x) \] - Leads to the _marginal_ distribution of $X(1)$, what we should predict for the first observation: \[ \Prob{X(1)=x} = \sum_{s}{\Prob{X(1)=x|S(1)=s}\Prob{S(1)=s}} = \sum_{s}{g(s, x) \Prob{S(1)=s}} \] - Now estimate the state at time $1$ _after_ seeing that $X(1)$ is actually $x(1)$: \begin{eqnarray} \Prob{S(1)=s|X(1)=x(1)} & = & \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^{\prime}} g(s^{\prime}, x(1))}} \equiv F_1(s) \end{eqnarray} ## A Recursion Trick (a.k.a. "The Forward Algorithm") [2] {.smaller} - We have the conditional distribution of $S(1)$, extrapolate it forward: \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} - The extrapolated state lets us predict $X(2)$: \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} - After we see $X(2)$, use both it _and_ $X(1)$ to update $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} ## A Recursion Trick (a.k.a. "The Forward Algorithm") [3] {.smaller} #### 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(t+1)) \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 - 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: Words {.smaller} - The key thing is to get $F_t(s) = \Prob{S(t)=s|X(1:t)=x(1:t)}$ - In the **Monte Carlo filter**, or **particle filter**, we approximate $F_t(s)$ by a large number of random samples from the right distribution #### The Particle Filter 1. At $t=1$, make $m$ independent draws $S^*_1, S^*_2, \ldots S^*_m$ from $\Prob{S(1)}$ + These are the **particles** 2. For each $S^*_i$, calculate $g(S^*_i, x(t)) \equiv w_{i}$ + $\Prob{X(t)=x(t)|X(1:t-1)=x(1:t-1)} \approx m^{-1} \sum_{i=1}^{m}{w_i}$ 3. Resample the particles: draw $m$ of them, with replacement, from $S^*$, with probabilities $\propto w_i$; call these $\tilde{S}_1, \tilde{S}_2, \ldots \tilde{S}_m$ 4. Set $\hat{F}_t$ to the sample distribution of the $\tilde{S}$ + The distribution of the resampled particles approximates $\Prob{S(t)|X(1:t)=x(1:t)}$ 5. Increment $t$ by 1, and **evolve** the particles forward in time: for each $\tilde{S}_i$, draw $S^*_i$ from $q(\tilde{S}_i, S^*_i)$ + Distribution of the new $S^*$ approximates $\Prob{S(t+1)|X(1:t)=x(1:t)}$. 6. Go to (2) ## The Particle Filter: Code ```{r, echo=TRUE, tidy=TRUE, tidy.opts=list(comment=FALSE)} # 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) } ``` ## The Particle Filter: For Our Little Demo ```{r, echo=TRUE, tidy=FALSE} 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 ```{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)) pf.means <- apply(sim.pf, 2, mean) lines(1:100, pf.limits[1,], col="blue", lty="dashed") lines(1:100, pf.limits[2,], col="blue", lty="dashed") lines(1:100, pf.means, col="blue") legend("bottomright", legend=c("Hidden state S(t)", "Observable X(t)", "Particle filter (mean)", "Particle filter (90% Interval)"), pch=c(NA, 1, NA), lty=c("solid", "blank", "solid"), "dashed", col=c("grey", "black", "blue", "blue"), cex=0.5) ``` ## 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 [@Big-Laplace-Gauss-filter-paper] ## 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-dirty-pictures + 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 [@Geometry-from-a-time-series] ## State-Space Reconstruction (cont'd.) ```{r} alpha <- runif(1, min=-pi, max=pi) omega <- rnorm(1) observable <- function(x) { cos(omega*x+alpha) } logistic.map <- function(s, r=4) { r*s*(1-s) } logistic.iteration <- function(n, init=runif(1), r=4) { s.vec <- vector(length=n) s <- init for (t in 1:n) { # record the current state s.vec[t] <- s # update the state s <- logistic.map(s, r=r) } return(s.vec) } ``` - $S(t) \in [0,1]$, $S(t+1) = 4 S(t) (1-S(t))$ + **Logistic map** - $X(t) = \cos{(`r omega` S(t) + `r alpha`)}$ + Randomly chosen frequency and phase shift ```{r, fig.cap="Vectors $(X(t), X(t+1), X(t+2))$, color-coded by the value of the initial condition $S(t)$."} s.seq <- seq(from=0, to=1, length.out=1000) s.mat <- t(sapply(s.seq, logistic.iteration, n=3)) x.mat <- observable(s.mat) library(scatterplot3d) scatterplot3d(x.mat, color=terrain.colors(length(s.seq)), pch=16, xlab="X(t)", ylab="X(t+1)", zlab="X(t+2)") ``` ## State-Space Reconstruction (cont'd.) - Fact from differential topology [@Whitney-embedding]: > **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-et-al-embedology]) - Application to this problem [@Takens-embedding]: > **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 ```{r} s.traj <- logistic.iteration(n=300) x.traj <- observable(s.traj) ``` ```{r} plot(s.traj, xlab="t", ylab="S(t)", pch=16, type="o") ``` ```{r} plot(x.traj, xlab="t", ylab="X(t)", pch=16, type="o", col=terrain.colors(1000)[ceiling(1000*s.traj)]) ``` ```{r} time.delay.vector <- function(x, k) { n <- length(x) delay.column <- function(i) { x[(1+i):(n-(k-1)+i)] } time.delay.matrix <- sapply((k-1):0, delay.column) return(time.delay.matrix) } ``` ```{r} y.traj <- time.delay.vector(x.traj, k=3) ``` ```{r} scatterplot3d(y.traj, xlab="X(t)", ylab="X(t-1)", zlab="X(t-2)", color=terrain.colors(1000)[ceiling(1000*head(s.traj, -2))], pch=16) ``` ## State-Space Reconstruction (cont'd.) - There is a whole alternative tradition of state-space modeling based on these theorems [@Kantz-Schreiber-2nd] - 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-conclique-Gibbs-sampler for a clever idea for speeding this up by simultaneous updating ## Backup: The Gibbs Sampler: History {.smaller} ```{r} # Do some calculations related to the size of the sample space in the original # Geman and Geman paper on the Gibbs sampler here, rather than in-line # They considered images which were 64 pixels * 64 pixels with 5 grey levels # per pixel, so the size of the state space was 5^(64*64), which R thinks # is infinity; calculate the logs (base 10) instead log.geman2.space <- 64*64*log(5, base=10) power.of.ten <- floor(log.geman2.space) prefactor <- 10^(log.geman2.space - power.of.ten) ``` - ... was not invented by Gibbs, but by @Geman-and-Geman + As they said (p. 722), it's very similar to older ideas in computational physics, going back to @Metropolis-et-al-Monte-Carlo - J. W. Gibbs was an American physicist who pioneered "statistical mechanics" --- building probability models for physical systems - Gibbs's models were based on the idea that the probability of the system being in state $x$ would be $\propto e^{- k U(x)/T}$, where $U(x)$ is the energy of the system in state $x$, $T$ is temperature, and $k$ is a universal constant + Normalizing, $\mathrm{Pr}(X=x) = \frac{e^{-k U(x)/T}}{\sum_{x^{\prime}}{e^{-k U(x^{\prime})/T}}} = Z^{-1}(T) e^{-k U(x)/T}$ + Computing the **partition function** $Z(T) = \sum_{x}{e^{-k U(x)/T}}$ is generally hard, because there are a very large number of states * Physicists have spent whole careers getting short-cut formulas for $Z(T)$ for particular models - Distributions like this came to be called "Gibbs distributions" or "Gibbs measures" + Statisticians would call them exponential families [@Mandelbrot-sufficiency-and-estimation-in-thermo] - In the 1970s, people realized that every Gibbs distribution corresponds to a Markov random field, and vice versa, every Markov random field corresponds to a Gibbs distribution [@Griffeath-on-random-fields] - @Geman-and-Geman were making statistical models of images, and modeled them as Gibbs distributions, where the "energy" favored things like "color changes slowly from one pixel to the next", but then they wanted / needed to generate images from their model + Directly sampling from $\mathrm{Pr}(X=x)$ would require calculating $Z(T)$, which was intractable even for their $64 \times 64$ images; with 5 level grey-scale, that's $5^{64 \times 64} = `r signif(prefactor, 2)` \times 10^{`r power.of.ten`}$ states! - Geman and Geman used the correspondence between Gibbs distributions and Markov random fields to get the conditional distribution of each pixel given the neighboring pixels + Each conditional distribution is easy to simulate from + Sweeping through all the conditional distributions a few times gives something close to the joint distribution ## 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) \] ## Backup: Iterated Conditional Modes (cont'd.) {.smaller} \[ \Prob{S(r)=s|X(r)=x, S(\Neighborhood(r)) = q} \propto g(x|s) p(s|q) \] #### The algorithm - Start with a guess at $\hat{S}(r)$ for each $r$ + Maybe $\hat{S}(r) \leftarrow \argmax_{s}{ f(x(r)|s) }$ - Go over the sites $r$ in arbitrary order: + $\hat{S}(r) \leftarrow \argmax_{s}{ f(x(r)|s) p(s| \hat{S}(\Neighborhood(r))) }$ * This is the (posterior) mode for $S(r)$, conditional on $X(r)$ and the guess about the rest of the Markov field, - Iterate to convergence + Each update never decreases $\Prob{S|X}$ for the whole fields, so this must converge [@Besag-dirty-pictures, p. 265] - Naturally, lots of variants - Everyone cites @Besag-dirty-pictures, but it was actually introduced a few years earlier by multiple authors [@Besag-dirty-pictures, p. 265] ## References