--- title: Hidden Markov Models and State Estimation author: 36-467/36-667 date: 19 November 2020 (Lecture 23) output: slidy_presentation bibliography: locusts.bib --- \[ \newcommand{\Prob}[1]{\mathbb{P}\left( #1 \right)} \] ```{r, include=FALSE} set.seed(2018-12-05) 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) ``` ## Housekeeping - Class is canceled for Tuesday + I'll post notes for what would have been the optional lecture - Homework 12 will be due on 3 December (as scheduled) ## Previously - 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) \perp 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 - $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 $\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)}\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 $\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$ Given $X(t)$, the only way $X(t-1)$ can matter for $X(t+1)$ is by changing the conditional distribution of $S(t)$ ## Relating $X(t)$ to $S(t)$ - 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)$ - 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 - 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) \perp X(t+1)|X(t)$ - That isn't true here - $\therefore$ the observed process isn't Markov ## Hidden Markov Models / State-Space Models - **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 + "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$? ## 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 all our simulation-based-inference tricks ## What About Space? What About Time? - Hidden Markov random field: $S(r)$ is a Markov random field, $X(r) = f(S(r), \epsilon(r))$ + $X(r)$ isn't Markov + Simulation: Gibbs for the $S$ field, then draw from $X(r)|S(r)$ independently + Likelihood: There are recursive approaches, though rather more complicated than the forward algorithm - 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))$ - OR: dynamic factor models + $S(t)$ is a low-dimensional vector + $S(t)$ evolves according to a Markov process (often assumed linear) + $X(r,t) = f_r(S(t), \epsilon(r,t))$, so $X(r,t) \perp X(q,t) | S(t)$ * This measurement model is also often assumed to be linear ## 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 - 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: Applying Indirect Inference to Our Demo - We need: 1. Simulator for the generative model 2. Estimator for the auxiliary model 3. Code to optimize the generative model so auxiliaries fit to simulations match the auxiliary fit to the data - We'll use AR(2) as the auxiliary + One more auxiliary parameter than we really need to estimate $q = \Prob{S(t+1)=+1|S(t)=+1} = \Prob{S(t+1)=-1|S(t)=-1}$ ## Backup: Finding the Auxiliary Parameter Estimates ```{r, echo=TRUE} # 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)) } ``` ## Backup: Can We Go From Auxiliary Parameters to Generative Parameters? ```{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") ``` ## Backup: Optimizing the Generative Model so Simulation-Auxiliaries Match Data-Auxiliaries ```{r, echo=TRUE} 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) } ``` - Run: ```{r, echo=TRUE} (toy.ii <- toy.indirect.inference(x=sim$x, starting.q=0.5, s=500)) ``` - Remember we set the true $q$ to $0.75$ - This isn't a fluke: ```{r, echo=TRUE} summary(replicate(30, toy.indirect.inference(x=sim$x, starting.q=0.5, s=500))) ``` + All the variation here is from simulation to simulation