---
title: Markov Random Fields
author: 36-467/36-667
date: 17 November 2020 (Lecture 20)
output: slidy_presentation
bibliography: locusts.bib
---

```{r, include=FALSE}
# General set-up options
library(knitr)
opts_chunk$set(size="small", background="white",
               cache=TRUE, autodep=TRUE,
               tidy=FALSE, warning=FALSE, message=FALSE,
               echo=FALSE)
```

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

## In our previous episodes

- Markov processes in time: $\Prob{X(t+1)|X(0), X(1), \ldots X(t)} = \Prob{X(t+1)|X(t)}$
    + Markov chains: $X(t)$ is discrete-valued
- Asymptotic behavior of Markov chains
- Likelihood-based inference for Markov models
- What about space?  What about space and time?

## Markov Random Fields

- $X(r) =$ value of the random field at site $r$
- $X$ or $X(\cdot) =$ value of the field at all sites, the **configuration**
- $X(-r) =$ value of the random field everywhere other than $r$
- $X(\Neighbors(r)) =$ value of the random field at the neighbors of $r$
    + the **neighborhood configuration**
- $X$ is a **Markov random field** when
\[
\Prob{X(r)|X(-r)} = \Prob{X(r)|X(\Neighbors(r))}
\]
- $p_r(x, y) \equiv \Prob{X(r)=x|X(\Neighbors(r))=y}$ is called the **local characteristic** for site $r$
    + People also write $q_{r}(x,y)$
	+ Can vary with parameters, $p_{r}(x,y;\theta)$

## What Does a Markov Random Field Look Like?

- Working out joint $\Prob{X}$ from $\Prob{X(r)|X(\Neighbors(r))}$ is not easy!
    + No natural order to use in factorizing
- Approach I: Heroic mathematics
    + Truly heroic in some cases...
- Approach II: Gibbs-Markov equivalence
    + Upshot: For the right functions $V_r$, $V_{rq}$, the **potentials**,
\[
\log{\Prob{X=x}} = \text{(constant)} + \sum_{r}{V_r(x(r))} + \sum_{r}{\sum_{q\in\Neighbors(r)}{V_{rq}(x(r), x(q))}}
\]
    + Getting the constant is hard
    + See backup
- Approach III: Gibbs sampler
    + Same Gibbs, different idea

## An Example: The Ising Model {.smaller}

- Space is 2D square lattice
- Each $X(r)$ is binary, say $+1$ or $-1$
- Potentials:
\[
V_{rq}(x(r), x(q)) = \theta x(r) x(q)
\]
    + "Homogeneous potentials"
- Add up potentials across neighbors:
\[
U_r(x|y) = \sum_{q \in \Neighbors(r)}{\theta x y(q)}
\]
- Local characteristics:
\[
p_r(x,y) = \Prob{X(r) = x|X(\Neighbors(r))=y} = \frac{e^{U_r(x|y)}}{e^{U_r(-1|y)} + e^{U_r(+1|y)}}
\]
    + "Autologistic" model

## The Gibbs Sampler

Reprise from [Lecture 12](http://www.stat.cmu.edu/~cshalizi/dst/20/lectures/12/lecture-12.html):

- 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(\Neighbors(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

## The Gibbs Sampler (actually running)

- The Ising model from before, run on a $150\times 150$ lattice, using the
  `mrf2d` package
    + 4 independent runs, all with $\theta=1$

```{r, echo=FALSE}
library(mrf2d)
# After the example in help(rmrf2d)
ising.theta <- array(rep(c(0,-1,-1,0),2), dim=c(2,2,2))
par(pty="s", mfrow=c(2,2))
image(rmrf2d(c(150,150), mrfi(1), ising.theta))
image(rmrf2d(c(150,150), mrfi(1), ising.theta))
image(rmrf2d(c(150,150), mrfi(1), ising.theta))
image(rmrf2d(c(150,150), mrfi(1), ising.theta))
```

## Changing dependence strength


```{r, echo=FALSE}
par(pty="s", mfrow=c(2,2))
image(rmrf2d(c(150,150), mrfi(1), 0.1*ising.theta), main=expression(theta==0.1))
image(rmrf2d(c(150,150), mrfi(1), 0.5*ising.theta), main=expression(theta==0.5))
image(rmrf2d(c(150,150), mrfi(1), 2*ising.theta), main=expression(theta==2))
image(rmrf2d(c(150,150), mrfi(1), 10*ising.theta), main=expression(theta==10))
```


## Inference: Basics

- If we have multiple, independent copies of the field
    + Can find $p_r(x,y)$ for each $r$
	+ Use conditional density estimation or even regression
- If $p_r(x,y)$ is the same for all $r$
    + Could do conditional density estimation or even regression
	+ Usual error statistics aren't valid because of dependence
- Can also do simulation-based inference
    + [Conditional auto-regressions](http://www.stat.cmu.edu/~cshalizi/dst/20/lectures/20/lecture-12.html) make good auxiliaries for indirect inference

## Inference: Likelihood (1)

- Full likelihood inference is _hard_
    + Need the joint distribution $\Prob{X=x}$
	+ Which means needing the normalizing factor as a function of the parameters
	+ Often need to do Monte Carlo just to get the likelihood [@Geyer-on-MCMCMLE;@Geyer-Thompson-MCMCMLE]
- What are our alternatives which stay close to likelihood?

## Inference: Likelihood vs. Matching Sufficient Statistics

- Summary statistics are **sufficient** when we can write the likelihood in terms of them and not the full data
- For Markov chains, the sufficient statistics were the counts of pairs of successive states ($N_{ij}$)
- For MRFs, the sufficient statistics are counts of clique configurations
- E.g., in Ising model, counts of $(-1,-1)$ for (self, neighbor), of $(-1,+1)$, of $(+1,-1)$, of $(+1,+1)$
- We can get these easily from the data
- Handy fact about exponential families: the MLE is the value of $\theta$ where the _expected_ value of the sufficient statistics equals the observed value
- Pick an initial value of $\theta$; simulate multiple times to get $\approx$ expected sufficient statistics; adjust to align these with the observed values (**stochastic approximation**)


	
## Inference: Likelihood (3)

- **Pseudo-likelihood**: Product of all the local characteristics $p_r$
\[
\prod_{r}{p_r(x(r)|x(\Neighbors(r)))} = \prod_{r}{\Prob{X(r)=x(r)|X(\Neighbors(r))=x(\Neighbors(r))}}
\]
    + $\neq$ the joint probability
    + Also called a **composite likelihood** (or pseudolikelihood with no hyphen...)
	+ Maximizing this is generally a consistent estimator, but not as statistically efficient as likelihood
        * See @Geyer-Thompson-MCMCMLE for head-to-head comparisons on the Ising model; may be unsurprising that pseudo-likelihood works better when interactions are weaker


## Inference: Uncertainty

- If you can get the likelihood, use the Hessian of the log-likelihood
- [Sandwich-covariance methods](http://www.stat.cmu.edu/~cshalizi/dst/20/lectures/14/lecture-14.html) work as usual (if we can calculate gradients and second derivatives)
- Parametric bootstrap works as usual
- Nonparametric bootstrap for stationary fields
    + Use rectangular blocks [@Lahiri-resampling-for-dependent]
	+ Or maybe weirder shapes [@Levina-Bickel-resampling-random-fields]

## Adding Time Back In

- Now have $X(r,t)$
    + Write $X(\cdot, t)$ for the whole configuration at time $t$
- Two choices for how we think about the Markov property:

1. $X(r,t)$ should be conditioned on $X(\Neighbors(r), t-1)$ and $X(r,t-1)$
    + So $X(r,t)$ and $X(q,t)$ are independent given $X(\cdot, t-1)$
2. Condition on $X(\Neighbors(r), t-1)$, $X(r,t-1)$ _and_ $X(\Neighbors(r), t)$
    + Neighbors at time $t$ are still dependent given $X(\cdot, t-1)$
- Sometimes called **recursive** vs. **simultaneous**

## Spatio-temporal Markov Random Fields

- Recursive (first option) fields are _much easier_
    + $\Prob{X(\cdot, t)|X(\cdot, t-1)} = \prod_{r}{\Prob{X(r,t)|X(\Neighbors(r), t-1), X(r,t-1)}}$
	+ Can actually calculate the likelihood
	+ Similar to a vector-valued Markov process, but with more structure to the transition probabilities (because they have to products of local characteristics)
	+ Can estimate $\Prob{X(r,t)|X(\Neighbors(r), t-1), X(r,t-1)}$ the same way as in a Markov chain
- Simultaneous (second option) fields are not as nice
    + Basically, back to spatial models with an extra coordinate
	+ Pseudo-likelihood often the best bet
    + Simultaneous dependence may be necessary if the time between measurements is slow compared to the time-scale of the actual dynamics

## Cellular Automata {.smaller}

- $X(r,t)$ discrete (usually finite), $r$ on a regular lattice, $t$ discrete
- Recursive (option 1) spatio-temporal Markov model
    + Sometimes conditionally deterministic
- These are _very_ expressive models
    + Some examples: [I](http://www.stat.cmu.edu/~cshalizi/462/lectures/10/10.pdf), [II](http://www.stat.cmu.edu/~cshalizi/462/lectures/11/11.pdf)
- A fun (?) exercise: code up the spatial SIRS model:
    + Each cell is in one of the 3 states S, I, R
    + An S can get infected from a _neighboring_ I, with some probability each time step
    + An I turns in to an R at some rate
    + Rs turn back in to Ss at some other rate
    + Tune the parameters right and you see traveling waves of infection, sometimes [quite pretty](http://psoup.math.wisc.edu/archive/recipe2.html):
		
![](http://psoup.math.wisc.edu/archive/soup2.gif)

## Summary

- In a (spatial) Markov random field, $X(r)$ is "screened off" from the rest of the field by its neighbors
    + Conditional distribution of $X(r)$ given neighbors = local characteristic of site $r$
- Local characteristics determine the joint distribution, but actually _solving_ for the joint distribution is hard
    + Gibbs-Markov equivalence gives a general result, but rarely easy thing to calculate with
- Estimation via
    + Local conditional model (under homogeneity or multiple realizations)
	+ Likelihood (if you can calculate it)
	+ Likelihood-ish objective functions
	+ Simulation-based inference
- Recursive spatio-temporal Markov random fields are, despite the name, actually a lot simpler, and a very powerful class of models


## Backup: Gibbs-Markov Theorem

- $X=$ the whole random field, the vector of $X(r)$ for all $r$
- $X$ has a **Gibbs distribution** when
\[
\Prob{X=x} = \frac{1}{Z}\exp{\left\{-\left(\sum_{A}{V_A(x)}\right)\right\}}
\]
    + $A =$ subsets of sites
	+ $V_A=$ **potential function**, depending only on the state of sites in A
	+ $Z \equiv \sum_{x}{\exp{\left\{-\left(\sum_{A}{V_A(x)}\right)\right\}}} =$ **partition function**
- A set of sites is a **clique** if every site in $A$ is a neighbor of every other site
- A potential is a **nearest-neighbor potential** if $V_A = 0$ whenever $A$ isn't a clique

## Backup: Gibbs-Markov Theorem

- On a finite set of sites, if $\Prob{X=x}> 0$ for all configurations $x$,
then $X$ is a Markov random field if and only if $X$ has a Gibbs
distribution with a nearest-neighbor potential
- Gibbs $\Rightarrow$ Markov: direct calculation
- Markov $\Rightarrow$ Gibbs: complicated!
    + Our readings from Guttorp sketch a proof due to @Griffeath-on-random-fields
- $Z$ is hard to calculate even from the Gibbs form
    + E.g., read [this](https://www.quantamagazine.org/the-cartoon-picture-of-magnets-that-has-transformed-science-20200624/) for the saga of calculating $Z$ for the 2D and 3D Ising model
- Sampling from Gibbs distributions:
    + Use Gibbs sampler (since it's Markov)
	+ Use [MCMC](http://www.stat.cmu.edu/~cshalizi/statcomp/13/lectures/16/mcmc.pdf), since that doesn't need $Z$

## Backup: Gibbs-Markov Theorem

- Statisticians usually call this the "Hammersley-Clifford" theorem
- It was proved simultaneously and independent by H&C, Besag, Griffeath, Grimmett, ...
- Priority is a mess, and Griffeath was one of my thesis advisers...
- "Gibbs-Markov" expresses the content better, and slights everyone equally


## Backup: A Conditional Likelihood for a Fraction of the Data

- On a square lattice, color the sites alternately blue and red
- Each red site has only blue neighbors
- White sites are independent given the blue sites (by Markov property)
- Multiply the local characteristics of the red sites to get a conditional likelihood
\[
\Prob{X(\text{red})=x(\text{red})|X(\text{blue})=x(\text{blue})} = \prod_{r \in \text{red}}{\Prob{X(r)=x(r)|X(\Neighbors(r)) = x(\Neighbors(r))}}
\]

- @Bartlett-spatial-pattern, pp. 27--28, calls this "coding"
	+ $=$ the "concliques" of @Kaplan-conclique-Gibbs-sampler
    + Adapts to other geometries (see Bartlett again)
	+ Efficiency compared to pseudolikelihood is unclear




## References

