\[ \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

The Gibbs Sampler

Reprise from Lecture 12:

• 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 et al. (2019) 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\)

Changing dependence strength

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 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 1991; Geyer and Thompson 1992)
• 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 and Thompson (1992) 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 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 2003)
  • Or maybe weirder shapes (Levina and Bickel 2006)

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

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 (1976)
• \(Z\) is hard to calculate even from the Gibbs form
  • E.g., read this 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, 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 (1975), pp. 27–28, calls this “coding”
  • \(=\) the “concliques” of Kaplan et al. (2019)
  • Adapts to other geometries (see Bartlett again)
  • Efficiency compared to pseudolikelihood is unclear