The Metropolis-Hastings algorithm by example

John Kerl

December 21, 2012

Abstract

The following are notes from a talk given to the University of Arizona Department of Mathematics Graduate Probability Seminar on February 14, 2008.

• Consider flipping a fair coin 100 times. The probability space has size 2^{100}, which is on the order of 10^{69}. Yet we can easily sample from this space: each event occurs with equal probability (1/2)^{100}, and all we have to do is conduct 100 independent Bernoulli experiments with p = 1/2.

• If the coin isn’t fair, not all events occur with equal probability. Yet it’s still easy to sample from this space: we conduct 100 independent Bernoulli experiments with 0 ≤ p ≤ 1.

• Now suppose the coin flips aren’t independent — they have some peer pressure. This system is called a one-dimensional Ising model and it’s non-trivial to sample from this space.

I will describe how the Metropolis-Hastings algorithm allows us to conduct experiments using such a model.

This paper is under construction.

Contents

1 Introduction  4

2 Probability distributions  4
2.1 IID fair flips . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2.2 IID flips . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.3 Independent flips . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.4 Dependent flips . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

3 The one-dimensional Ising model  7

4 Independence and CDF sampling  9
4.1 Sample proportions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
4.2 Sampling with the independence method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
4.3 Sampling with the CDF method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

5 Sampling with Metropolis-Hastings  13
5.1 The algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
5.2 Walk-through . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
5.3 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
5.4 Heuristics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
5.5 Random variables for 1D Ising . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

6 Design choices and caveats for Metropolis-Hastings  17
6.1 Number of trials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
6.2 Initial state . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
6.3 Sequential vs. random sweeps . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
6.4 Thermalization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
6.5 Separation of states . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
6.6 Downsampling to avoid correlation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
6.7 ΔE computations for the independent model . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
6.8 ΔE computations for the 1D Ising model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

7 Markov chains  21

7.1 The die-tipping experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
7.2 Stochastic processes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
7.3 Conditional probability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
7.4 Markov chains . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
7.5 Die-tipping revisited . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
7.6 Markov-chain properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

8 Markov chain Monte Carlo  27
8.1 Metropolis choice of Markov matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
8.2 Sampling theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

9 Generalizations  28
9.1 2D Ising . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
9.2 Non-example: percolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
9.3 Numerical integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

References  29

Index  30

1 Introduction

Coin flips are used as a motivating example to describe why one would want to use the Metropolis-Hastings algorithm. The algorithm is presented, illustrated by example, and then proved correct. See [Kerl] for probability terminology and notation used in this paper.

2 Probability distributions

Consider flipping n coins. Encode a heads-up flip as +1 and a tails-up flip as −1. Then the result of the ith flip is ω_{i} ∈ {+1, −1}. I’ll discuss four scenarios, with increasing levels of complexity. Write ordered n-tuples

ω = (ω_{1}, …, ω_{n}).

Notice that there are N = 2^{n} outcomes in the probability space

Ω = {(ω_{1}, …, ω_{n}) : ω_{i} = ±1}.

For example, with n = 3 coins, there are N = 8 outcomes: HHH, HHT, HTH, HTT, THH, THT, TTH, and TTT.

2.1 IID fair flips

In the first case I consider, the flips are independent and identically distributed. Furthermore, each coin is fair, which is to say

P_{i}(+1) = 1/2  and  P_{i}(−1) = 1/2.

Due to independence, the probabilities factor:

P (ω) = P (ω_{1}, …, ω_{n})

= P_{1}(ω_{1}) ⋯ P_{n}(ω_{n})

= ∏_{i=1}^{n} \frac{1}{2} = (\frac{1}{2})^{n} .

For example, with n = 3, the probability of heads, tails, heads is

P (+1, −1, +1) = 1/8.

The other 7 possible outcomes also occur with probability 1/8. Here’s the PMF:

(ω_{1}, ω_{2}, ω_{3})  P (ω_{1}, ω_{2}, ω_{3})
( 1 1 1 )  0.125
( 1 1 -1 )  0.125
( 1 -1 1 )  0.125
( 1 -1 -1 )  0.125
( -1 1 1 )  0.125
( -1 1 -1 )  0.125
( -1 -1 1 )  0.125
( -1 -1 -1 )  0.125

2.2 IID flips

In the second case I consider, the n flips are still IID but the coin isn’t necessarily fair. Each lands heads-up with probability p, where 0 ≤ p ≤ 1. Due to independence, the probabilities still factor. Let

η_{i} = \begin{cases}0, & ω_{i} = +1 \\ 1, & ω_{i} = −1. \end{cases}  (2.1)

That is,

ω_{i} = 1 − 2η_{i} and η_{i} = \frac{1 − ω_{i}}{2} . (2.2)

Then

P (ω) = P (ω_{1}, …, ω_{n})

= P_{1}(ω_{1}) ⋯ P_{n}(ω_{n})

= ∏_{i=1}^{n} p^{1−η_{i}}(1 − p)^{η_{i}}.

With n = 3, the probability of heads, tails, heads is

P (1, −1, 1) = (p^{1}(1 − p)^{0}) (p^{0}(1 − p)^{1}) (p^{1}(1 − p)^{0})
= p^{2}(1 − p).

The PMF is as follows, with p = 0.9:

(ω_{1}, ω_{2}, ω_{3})  P (ω_{1}, ω_{2}, ω_{3})
( 1 1 1 )  0.729
( 1 1 -1 )  0.081
( 1 -1 1 )  0.081
( 1 -1 -1 )  0.009
( -1 1 1 )  0.081
( -1 1 -1 )  0.009
( -1 -1 1 )  0.009
( -1 -1 -1 )  0.001

2.3 Independent flips

In the third case I consider, the coin flips are still independent but they are not identically distributed: each has its own probability p_{i} of heads. That is,

P_{i}(ω_{i} = +1) = p_{i}.

Then

\begin{aligned}P (ω) & = P (ω_{1}, …, ω_{n}) \\ & = P_{1}(ω_{1}) ⋯ P_{n}(ω_{n}) \\ & = ∏_{i=1}^{n} p_{i}^{1−η_{i}}(1 − p_{i})^{η_{i}}. \end{aligned}  (2.3)

With n = 3, the probability of heads, tails, heads is

P (1, −1, 1) = p_{1}(1 − p_{2})p_{3}.

With p_{1} = 0.6, p_{2} = 0.7, and p_{3} = 0.8, the PMF looks like this:

(ω_{1}, ω_{2}, ω_{3})  P (ω_{1}, ω_{2}, ω_{3})
( 1 1 1 )  0.336
( 1 1 -1 )  0.084
( 1 -1 1 )  0.144
( 1 -1 -1 )  0.036
( -1 1 1 )  0.224
( -1 1 -1 )  0.056
( -1 -1 1 )  0.096
( -1 -1 -1 )  0.024

2.4 Dependent flips

In the fourth and final case I consider, the coin flips are no longer independent — perhaps we are flipping magnetized coins, where the result of a flip is influenced by one or more of its neighbors. The PMF is a function with potentially 2^{n} different input/output values. There are all sorts of PMFs one could write down to define a probability model for dependent coin flips. For example:

(ω_{1}, ω_{2}, ω_{3})  P (ω_{1}, ω_{2}, ω_{3})
( 1 1 1 )  0.4
( 1 1 -1 )  0.0
( 1 -1 1 )  0.1
( 1 -1 -1 )  0.1
( -1 1 1 )  0.0
( -1 1 -1 )  0.1
( -1 -1 1 )  0.1
( -1 -1 -1 )  0.2

However, it’s not clear what kind of physical situation could lead to such an ad-hoc PMF.

There is one particular family of models — the Ising family — which contains enough flexibility to describe plenty of dependent-flip situations. Furthermore, the Ising family has some physical plausibility.

3 The one-dimensional Ising model

I’ll refer to spins rather than flips: heads and tails will be replaced with spin-up and spin-down, respectively. A sequence such as HHTH will now be + + −+, and I will call that a spin configuration.

In the one-dimensional Ising model, there are n spins and N = 2^{n} possible states (spin configurations), and the individual spins in a configuration are not necessarily independent.

Let ω = (ω_{1}, …, ω_{n}) as before. We define an energy function on the spins ω_{i} by

E(ω) = ∑_{i=1}^{n} ∑_{j=1}^{n} S_{ij}ω_{i}ω_{j} + ∑_{i=1}^{n} h_{i}ω_{i}.

The S_{ij}’s are coupling coefficients which specify the dependence of one spin on another. The h_{i}’s are the magnetic field term, which breaks the symmetry.

We get a probability model using the statistical-mechanical notation that energy is the logarithm of probability. The probability of each configuration ω is a Boltzmann distribution with inverse temperature β = 1/kT (where k is Boltzmann’s constant):

P (ω) ∝ e^{−βE(ω)}.  (3.1)

To normalize this, so that ∑_{ω} P (ω) = 1, we divide by

Z = ∑_{k=1}^{N} e^{−βE(ω^{(k)})}

where ω^{(k)} runs over all possible states. This denominator has a special name: it is called a partition function. Thus,

P (ω^{(j)}) = \frac{e^{−βE(ω^{(j)})}}{∑_{k=1}^{N} e^{−βE(ω^{(k)})}} .

Some handy notation: Write the probability numerator as

P^{∗}(ω^{(j)}) = e^{−βE(ω^{(j)})}.  (3.2)

Then

P (ω^{(j)}) = \frac{P^{∗}(ω^{(j)})}{Z}.  (3.3)

We will see (if I get to it) that the diagonal terms S_{ii} do not affect spin probabilities. Here are some examples of S matrices, with n = 3.

Non-interacting:  Mean field:  Aligning nearest-neighbor:
S = ( \begin{matrix}0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{matrix} )  S = ( \begin{matrix}0 & 1 & 1 \\ 1 & 0 & 1 \\ 1 & 1 & 0 \end{matrix} )  S = ( \begin{matrix}0 & 1 & 0 \\ 0 & 0 & 1 \\ 0 & 0 & 0 \end{matrix} )

Anti-aligning nearest-neighbor:  Aligning nearest-neighbor  Anti-aligning nearest-neighbor
with periodic boundary conditions: with periodic boundary conditions:
S = ( \begin{matrix}0 & −1 & 0 \\ 0 & 0 & −1 \\ 0 & 0 & 0 \end{matrix} )  S = ( \begin{matrix}0 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \end{matrix} )  S = ( \begin{matrix}0 & −1 & 0 \\ 0 & 0 & −1 \\ −1 & 0 & 0 \end{matrix} )

Here are the PMFs for the first three of those, with h = 0:

(ω_{1}, ω_{2}, ω_{3})  P (ω_{1}, ω_{2}, ω_{3}) P (ω_{1}, ω_{2}, ω_{3})  P (ω_{1}, ω_{2}, ω_{3})
Non-interacting Mean field Aligning nearest neighbor
( 1 1 1 )  0.125  0.4994973  0.3879017
( 1 1 -1 )  0.125  0.0001676  0.0524968
( 1 -1 1 )  0.125  0.0001676  0.0071047
( 1 -1 -1 )  0.125  0.0001676  0.0524968
( -1 1 1 )  0.125  0.0001676  0.0524968
( -1 1 -1 )  0.125  0.0001676  0.0071047
( -1 -1 1 )  0.125  0.0001676  0.0524968
( -1 -1 -1 )  0.125  0.4994973  0.3879017

4 Independence and CDF sampling

In this section we discuss two sampling methods which are simpler than Metropolis-Hastings: independence sampling, which works for the independent cases, and CDF sampling which works for all four examples of section 2 but doesn’t scale well as n increases.

4.1 Sample proportions

Suppose we flip four fair coins. Getting all four heads, or all four tails, is no huge surprise — but getting two heads and two tails is more likely. On the other hand, if we flip a thousand fair coins, all heads or all tails are very unlikely and we expect about 500 heads in general. It is quite possible to quantify these expectations precisely, in terms of sample mean and variance, using the binomial theorem as described in [Kerl]: IID fair flips lend themselves to neat theoretical analysis. For purposes of this paper, though, the goal is to describe experimental methods which are useful for probability models which are not theoretically neat. It is more important to look at sample proportions, which allow us to tabulate our experimental results.

Suppose we conduct (in a manner which subsequent sections will make precise) M experiments of an n-coin flip, using any of the four models of section 2, with outcomes

ω^{(1)}, …, ω^{(M)}.

Then, for j = 1, …, N = 2^{n}, define Q_{i} to be the fraction of ω’s which landed in state i. I also sometimes write, for brevity,

P = (P_{1}, …, P_{N}) and Q = (Q_{1}, …, Q_{N}).

For example, if I flip N = 3 fair coins M = 10 times, I might get the following:

i  ω^{(i)}
1 ( 1 1 1 )
2 ( -1 1 -1 )
3 ( -1 1 1 )
4 ( -1 1 -1 )
5 ( -1 -1 1 )
6 ( 1 -1 1 )
7 ( 1 -1 1 )
8 ( -1 -1 1 )
9 ( 1 -1 -1 )
10 ( 1 1 1 )

Then I can tabulate the sample proportions, which I can compare against the PMF:

ω  P (ω) Q(ω)
( 1 1 1 ) 0.125 0.2
( 1 1 -1 ) 0.125 0.0
( 1 -1 1 ) 0.125 0.2
( 1 -1 -1 ) 0.125 0.1
( -1 1 1 ) 0.125 0.1
( -1 1 -1 ) 0.125 0.2
( -1 -1 1 ) 0.125 0.2
( -1 -1 -1 ) 0.125 0.0

If I re-run this experiment with M = 10000 trials, I might get the following sample proportions:

ω  P (ω) Q(ω)
( 1 1 1 ) 0.125 0.1247
( 1 1 -1 ) 0.125 0.1191
( 1 -1 1 ) 0.125 0.1250
( 1 -1 -1 ) 0.125 0.1281
( -1 1 1 ) 0.125 0.1237
( -1 1 -1 ) 0.125 0.1307
( -1 -1 1 ) 0.125 0.1231
( -1 -1 -1 ) 0.125 0.1256

The law of large numbers tells us that each of the eight Q_{i}’s approach the respective P_{i}’s as the number M of trials increases; the central limit theorem tells us about the distribution of the Q_{i}’s.

The next question is how to actually conduct such an experiment in a computer program.

4.2 Sampling with the independence method

For independent flips, e.g. any of the first three models in section 2, it’s easy to produce samples. Recall that p_{j} is the probability that the jth flip is heads, which we encode as +1. Thus, the algorithm is straightforward.

Initialize the histogram bins: For i = 1, …, N (number of possible outcomes):

counts_{i} = 0.

Generate states and collect data: For i = 1, …, M (number of trials):

For j = 1, …, n (number of coins): U = uniform random on [0, 1). If U < p_{j} then ω_{j} = +1 else ω_{j} = −1.

Find the index k of the state ω = (ω_{1}, …, ω_{n}) which was just generated. Increment counts_{k} by 1.

Scale the histogram bins down to sample proportions: For i = 1, …, N (number of possible outcomes):

Divide counts_{i} by M.

The output of such an experiment was shown in section 4.1.

The main limitation of this algorithm is obvious: it only works for independent-flip models.

4.3 Sampling with the CDF method

For all four models of section 2 — independent or not — it is possible to sample using a technique which is almost as simple as the independence sampling of the previous section.

The key concept is the cumulative distribution function. That is, we write down the N possible outcomes in some order ω^{(1)} through ω^{(N)}, then tabulate at the jth slot the probability of obtaining states ω^{(1)} through ω^{(j)}. That is, the CDF, as a function of j from 1 to N, is

C(j) = P ({ω^{(1)}, …, ω^{(j)}}).

E.g.

C(1) = P (ω^{(1)}
C(2) = P (ω^{(1)} + P (ω^{(2)}
C(3) = P (ω^{(1)} + P (ω^{(2)} + P (ω^{(3)}
⋮
C(N) = P (ω^{(1)} + … + P (ω^{(N)}.

For example, using the IID-fair model, we have

ω  P (ω) CDF
( 1 1 1 ) 0.125 0.125
( 1 1 -1 ) 0.125 0.250
( 1 -1 1 ) 0.125 0.375
( 1 -1 -1 ) 0.125 0.500
( -1 1 1 ) 0.125 0.625
( -1 1 -1 ) 0.125 0.750
( -1 -1 1 ) 0.125 0.875
( -1 -1 -1 ) 0.125 1.000

Likewise, using the aligning nearest-neighbor Ising model from section 3, with h = 0.1, we have

ω  P (ω)  CDF
( 1 1 1 ) 0.3879017 0.3879017
( 1 1 -1 ) 0.0524968 0.4403985
( 1 -1 1 ) 0.0071047 0.4475032
( 1 -1 -1 ) 0.0524968 0.5000000
( -1 1 1 ) 0.0524968 0.5524968
( -1 1 -1 ) 0.0071047 0.5596015
( -1 -1 1 ) 0.0524968 0.6120983
( -1 -1 -1 ) 0.3879017 1.0000000

The algorithm is as follows.

Turn the PMF into a CDF: Sum = 0. For i = 1, …, N (number of possible states):

Sum := sum + p_{i}.

CDF_{i} = sum.

Initialize the histogram bins: For i = 1, …, N (number of possible outcomes):

counts_{i} = 0.

Generate states and collect data: For i = 1, …, M (number of trials):

U = uniform random on [0, 1). For j = 1, …, N:

If U < CDF_{j}:

Increment counts_{j} by 1; go the next iteration of the i loop.

Scale the histogram bins down to sample proportions: For i = 1, …, N (number of possible outcomes):

Divide counts_{i} by M.

Example results are as follows, with M = 10000 trials:

ω  P (ω)  Q(ω) Difference
( 1 1 1 ) 0.3879017 0.3832000 0.0047017
( 1 1 -1 ) 0.0524968 0.0480000 0.0044968
( 1 -1 1 ) 0.0071047 0.0070000 0.0001047
( 1 -1 -1 ) 0.0524968 0.0506000 0.0018968
( -1 1 1 ) 0.0524968 0.0557000 −0.0032032
( -1 1 -1 ) 0.0071047 0.0063000 0.0008047
( -1 -1 1 ) 0.0524968 0.0522000 0.0002968
( -1 -1 -1 ) 0.3879017 0.3970000 −0.0090983

The limitation of this CDF-sampling algorithm is a bit subtle: it requires computing probabilities for all N states. In the coin-flip example, N = 2^{n}. For, say, n = 100, 2^{n} ≈ 10^{69} which is (as I recall) on the order of the number of protons in the universe. One cannot even write down all possible states. Yet 100 coins isn’t a lot so it seems like we should be able to do something.

There is even worse news. Computing the partition function Z requires a sum over all N states, which is also infeasible. Not only can we not list the probability of all states, we can’t even compute the (normalized) probability of even one state, since (from equations 3.2 and 3.3) P (ω) = P^{∗}(ω)/Z. It would be nice if we had something which only depended on the probability numerator P^{∗}(ω) = e^{−βE(ω)}, which is easy to compute.

We will see in the following sections that the Metropolis-Hastings algorithm overcomes both these problems. First we will see what that algorithm is and what it does; then, we’ll see why it works.

5 Sampling with Metropolis-Hastings

Here we present the essential components of the Metropolis-Hastings algorithm, in pseudocode and a worked example. Important, non-negligible practical considerations are deferred until section 6; correctness is proved in section 8.

5.1 The algorithm

The goal is to select M samples ω^{(1)}, …, ω^{(M)} from the PMF of the 1D Ising model with specified S, h, and β. The states have probability numerator P^{∗}(ω) = exp(−βE(ω)) as given by equation 3.2. Here, I display spins as up (filled) or down (hollow):

• There is an inverse temperature β: For physical systems, this has some meaning and, presumably, would be varied under some controlled pattern. For the independent models, β cancels and has no effect (as we’ll see in section 6.7) so you can give it any non-zero value, e.g. β = 1.

• Pick an initial configuration. Typically, there are three choices: (1) Start with all spins down, i.e. ω = (−1, …, −1). (2) Start with all spins up, i.e. ω = (+1, …, +1). (3) Start with ω selected from a uniform probability distribution on Ω. (See section 6.2 for more information on choice of initial state.)

• Select a site j and decide whether to flip ω_{j} to −ω_{j}. (See section 6.3 for more information on site selection.)

PSfrag replacements P (change) = min{1, e^{−ΔH}}

• This decision is made using the Metropolis prescription, namely:

– Compute the change in energy ΔE = E(ω′) − E(ω) which would be obtained if ω were sent to ω′ by flipping ω_{j}.

– One may compute ΔE by separately computing E(ω′) and E(ω) and subtracting the two. However, since the only change is at the site j, one may do some ad-hoc algebra (sections 6.7 and 6.8) to derive an expression for ΔE which is less computationally expensive.

– Accept the change with probability

min{1, e^{−ΔE}}.

Otherwise, the change is said to be rejected.

This is called a Metropolis step.

• Looping through all n sites from j = 1 to j = n, performing a Metropolis step at each site j, is called a Metropolis sweep.

• Doing M sweeps is a complete experiment. (See section 6.1 for more information on choice of M.)

• If one realizes a random variable X(ω) at each of M sweeps, averaging X over the M sweeps, one obtains an approximation \overline{X} for the expectation E[X]. (See section 5.5 for more information on choice of random variable.)

The above may be rephrased in pseudocode as follows:

Pick the initial state ω: For i = 1, …, n (number of spins):

ω_{i} = 0 (spins-down) OR For i = 1, …, n (number of spins):

ω_{i} = 1 (spins-up) OR For i = 1, …, n (number of spins):

ω_{i} = ±1 with probabilities 1/2, 1/2 (uniform random initial state).

Generate states and collect data: For i = 1, …, M (number of trials):

One Metropolis sweep:

For j = 1, …, n (number of spins): One Metropolis step: ΔE = E(ω′) − E(ω) where ω′ has opposite spin at site j. Let V = min{1, e^{−βΔE}}. Generate U uniformly distributed on [0, 1). If U < V:

Accept the new state: ω := ω′. Else:

Reject the new state: keep ω. Find the index k of the state ω = (ω_{1}, …, ω_{n}) which was just generated. Increment counts_{k} by 1. (Do this whether the new state was accepted or not.)

Scale the histogram bins down to sample proportions: For k = 1, …, N (number of possible outcomes):

Divide counts_{k} by M.

Notice that a step of this algorithm touches one of the n sites; a sweep touches all n sites. We collect data only once per sweep, since only then have all sites been given a chance to change.

This is the basic algorithm. Important caveats, though, are presented in section 6.

5.2 Walk-through

Here are two sweeps with n = 3, β = 1, aligning nearest-neighbor model, initial state (1, 1, 1), and h = 0.1. These are the energies and probability numerators:

State number  ω  E(ω) P^{∗}(ω)
1  ( 1 1 1 ) -1.7 5.4739474
2  ( 1 1 -1 ) 0.1 0.9048374
3  ( 1 -1 1 ) 2.1 0.1224564
4  ( 1 -1 -1 ) -0.1 1.1051709
5  ( -1 1 1 ) 0.1 0.9048374
6  ( -1 1 -1 ) 1.9 0.1495686
7  ( -1 -1 1 ) -0.1 1.1051709
8  ( -1 -1 -1 ) -2.3 9.9741825

Initial state:       ( 1  1  1 )
Sweep 1:
omega       j  state#  Delta E  U         V         Accept/reject New E
------------ --- ------  -------  --------- --------- ------------- -----
( 1  1  1 )   1   1      1.8      0.9731406 0.1652989 Reject        -1.7
( 1  1  1 )   2   1      3.8      0.3204335 0.0223708 Reject        -1.7
( 1  1 -1 )   3   2      1.8      0.1590330 0.1652989 Accept         0.1
State after sweep 1: ( 1  1 -1 )
Sweep 2:
omega       j  state#  Delta E  U         V         Accept/reject New E
------------ --- ------  -------  --------- --------- ------------- -----
( 1  1 -1 )   1   2      1.8      0.6781967 0.1652989 Reject         0.1
( 1 -1 -1 )   2   4     -0.2      0.9644570 1.0000000 Accept        -0.1
( 1 -1 -1 )   3   4      2.2      0.2357908 0.1108032 Reject        -0.1
State after sweep 2: ( 1 -1 -1 )

5.3 Numerical results

Example results are as follows, n = 3, aligning nearest-neighbor model, h = 0.1, and M = 10000 trials:

ω  P (ω)  Q(ω) Difference
( 1 1 1 ) 0.3879017 0.3875000 0.0004017
( 1 1 -1 ) 0.0524968 0.0514000 0.0010968
( 1 -1 1 ) 0.0071047 0.0066000 0.0005047
( 1 -1 -1 ) 0.0524968 0.0523000 0.0001968
( -1 1 1 ) 0.0524968 0.0530000 −0.0005032
( -1 1 -1 ) 0.0071047 0.0070000 0.0001047
( -1 -1 1 ) 0.0524968 0.0523000 0.0001968
( -1 -1 -1 ) 0.3879017 0.3899000 −0.0019983

5.4 Heuristics

Recall from section 5.1 that the Metropolis transition probability is

min{1.0, exp(−βΔ_{E})}.

So, the Metropolis algorithm certainly transitions to states with non-negative exp(−βΔE), and randomly transitions to states with negative exp(−βΔE). The certain transition occurs when

exp(−βΔE) ≥ 1
−βΔE ≥ log(1) = 0
βΔE ≤ 0
ΔE ≤ 0.

Thus, we may also say that the Metropolis algorithm certainly transitions to lower-energy or same-energy states, and randomly transitions to higher-energy states.

Certain transition exp(−βΔE) ≥ 1 ΔE ≤ 0

Random transition exp(−βΔE) < 1 ΔE > 0

5.5 Random variables for 1D Ising

Here, the RVs are counts Q_{i} as described in section 4.1. This is just for didactic purposes. Often people look at sample means (averages) of some quantity, taken at each sample. One can also define energy, total magnetization, correlations between specified pairs of sites, etc.

6 Design choices and caveats for Metropolis-Hastings

For clarity, I presented the basic algorithm in section 5.1. Here are additional practical considerations which must be addressed in any serious use of Metropolis-Hastings. Correctness of the algorithm is proved in section 8.

6.1 Number of trials

Choice of the number of iterations M: How to do this wisely is a story in itself. For this paper, I suggest trying a few runs with each M, for varying M’s. For example, do half a dozen runs with M = 10000, another half-dozen runs with M = 20000, etc. Then choose an M which reduces the standard deviation of the sample mean to your desired number of decimal places.

6.2 Initial state

Choice of the initial state is as described in section 5.1: typically lowest-energy/highest-energy (e.g. here the all-down/all-up states), or random configuration. Results should be independent of initial state. Failure to achieve this may be due to inadequate thermalization (section 6.4).

6.3 Sequential vs. random sweeps

The decision here is whether to sweep through all sites from (say) left to right, or to move through them randomly. There is some debate in the literature on this topic (for which I lack citations here).

6.4 Thermalization

One should first run L Metropolis sweeps of the system, discarding any data collected, before running the M sweeps in which data are accumulated. The L sweeps are called the thermalization phase; the M sweeps are called the accumulation phase.

xxx clarify with numerical/graphical examples (describe in terms of P plots in section xxx): The underlying concern is the convergence of the Metropolis probability distribution to the stable distribution of its implicit Markov chain.

Put up E plots?

According to Tom Kennedy, there is no general method to determine whether the system has thermalized. xxx discuss E-plot idea (if it works out ...); manual approaches.

6.5 Separation of states

xxx highly annp — the distribution is supported on all-spins-up and all-spins-down. The Metro one-at-a-time method will always give rejects. Generalize, and cite some further reading.

6.6 Downsampling to avoid correlation

Downsampling and batched means are discussed thoroughly in appendix B of [Kerl2010]. The discussion there also includes information about Berg’s autocorrelation detection.

6.7 ΔE computations for the independent model

As noted in section 5.1, it is usually possible to compute ΔE at less computational cost than E(ω′) − E(ω). That is, computing E(ω) requires visiting all n sites of a spin configuration, as does computing E(ω′), but since the Metropolis algorithm only modifies one spin, we can skip most of the redundant values which will only cancel out anyway.

When you use a Metropolis-Hastings on a new problem, you will likely be doing some calculations of the type shown here.

xxx note β shouldn’t affect probabilities — it shouldn’t matter whether we’re flipping hot coins or cold ones. then do the algebra and show why this intuition is right.

xxx independence from temperature.

Using equations 3.1 and 2.3,

P (ω) ∝ e^{−βE(ω)}
xxx
E(ω) = −\frac{1}{β} log(P (ω))

= −\frac{1}{β} log (∏_{i=1}^{n} p_{i}^{1−η_{i}}(1 − p_{i})^{η_{i}})

= −\frac{1}{β} ∑_{i=1}^{n} ((1 − η_{i}) log(p_{i}) + η_{i} log(1 − p_{i})) .

xxx flip only at jth slot: ω′_{j} = −ω_{j} and so η′_{j} = 1 − η_{j}:

ΔE = −\frac{1}{β} (η_{j} log(p_{j}) + (1 − η_{j}) log(1 − p_{j}) − (1 − η_{j}) log(p_{j}) − η_{j} log(1 − p_{j}))

= −\frac{1}{β} ((1 − 2η_{j}) log(1 − p_{j}) − (1 − 2η_{j}) log(p_{j}))

= −\frac{1}{β} (1 − 2η_{j}) log (\frac{1 − p_{j}}{p_{j}} )

= \frac{ω_{j}}{β} log (\frac{p_{j}}{1 − p_{j}} ) .

xxx (or perhaps just state that the β’s cancel):

exp (−βΔE) = exp (−ω_{j} log (\frac{p_{j}}{1 − p_{j}} )) = (\frac{1 − p_{j}}{p_{j}} )^{ω_{j}} .

xxx certain-transition case is

(\frac{1 − p_{j}}{p_{j}} )^{ω_{j}} ≥ 1.

Thus, when ω_{j} = 1, there is a certain transition when p_{j} ≤ 1/2; when ω_{j} = 1, there is a certain transition when p_{j} ≥ 1/2. Since p was P (ω_{j} = 1) — the probability of heads at the jth slot — the Metropolis algorithm certainly sends individual spins to their more likely states, and randomly sends them to their less likely states.

6.8 ΔE computations for the 1D Ising model

As discussed at the beginning of the previous section, and as mentioned in section 5.1, we can compute ΔE at less computational cost than E(ω′) − E(ω).

The current state is ω and the energy is

E = − ∑_{ij} S_{ij}ω_{i}ω_{j} + ∑_{i} h_{i}ω_{i}.

The candidate new state ω′ has opposite spin ω_{k}. For example, with n = 4, k = 2, and

ω = (+1, +1, +1, +1),

we have (numbering states from 0 up to n − 1):

ω′ = (+1, +1, −1, +1).

Then the new energy is

E′ = − ∑_{ij} S_{ij}ω′_{i}ω′_{j} + ∑_{i} h_{i}ω′_{i}

and the energy change is

ΔE = E′ − E

= − ∑_{ij} S_{ij}ω′_{i}ω′_{j} + ∑_{ij} S_{ij}ω_{i}ω_{j} + ∑_{i} h_{i}(ω′_{i} − ω_{i}).

Now, in my computer program I’m going to be computing this energy difference very many times, so I want to omit as much redundant computation as possible. For i ≠ k we have ω′_{i} = ω_{i}, and for i = k we have ω′_{k} = −ω_{k}. So the last term is

h_{k}(ω′_{k} − ω_{k}) = −2h_{k}ω_{k}.

The sums split into four pieces: (1) i ≠ k and j ≠ k, (2) i = k and j ≠ k, (3) i ≠ k and j = k, and (4) i = k and j = k.

The first piece is

− ∑_{i,j≠k} S_{ij}ω′_{i}ω′_{j} + ∑_{i,j≠k} S_{ij}ω_{i}ω_{j} = − ∑_{i,j≠k} S_{ij}ω_{i}ω_{j} + ∑_{i,j≠k} S_{ij}ω_{i}ω_{j} = 0.

The second piece is

− ∑_{j≠k} S_{kj}ω′_{k}ω′_{j} + ∑_{j≠k} S_{kj}ω_{k}ω_{j} = ∑_{j≠k} S_{kj}ω_{k}ω_{j} + ∑_{j≠k} S_{kj}ω_{k}ω_{j} = 2 ∑_{j≠k} S_{kj}ω_{k}ω_{j}.

The third piece is

− ∑_{i≠k} S_{ik}ω′_{i}ω′_{k} + ∑_{i≠k} S_{ik}ω_{i}ω_{k} = 2 ∑_{i≠k} S_{ik}ω_{i}ω_{k}.

The fourth piece is

−S_{kk}ω′_{k}ω′_{k} + S_{kk}ω_{k}ω_{k} = −S_{kk}ω_{k}ω_{k} + S_{kk}ω_{k}ω_{k} = 0.

Thus the energy difference is

2 (∑_{j≠k} S_{kj}ω_{k}ω_{j} + ∑_{i≠k} S_{ik}ω_{i}ω_{k} − h_{k}ω_{k}) .

7 Markov chains

In order to prove the correctness of Metropolis-Hastings, we first need to understand Markov chains. We get Markov chains by putting together two concepts: stochastic processes and conditional probability. I preface discussion of Markov chains with the die-tipping example which helps motivate the discussion.

Throughout, I assume elementary probability terminology as in [Kerl, Lawler, GS]. Namely: sample space; probability measure; probability mass function (PMF); probability density function (PDF); random variable.

7.1 The die-tipping experiments

Most of this document uses coin-flipping and the 1D Ising model as a running example. Here, though, I describe two even simpler experiments (experiment A and experiment B), which involve tipping a single die. As simple as these are, they encapsulate all the central points I want to make about Markov chains. I’ll return to coin flips and the Ising model in section xxx.

First recall that the pips on opposite faces of a die add to 7:

Experiment A: I set the die on the table with the one-face up. Then I close my eyes momentarily and tip it, so the one-face is now on the side. I then close my eyes and tip it again, and so on, until I’ve tipped it 20 times.

If I repeat this experiment — placing the side and doing the 20 tips — many times, what’s the probability distribution for the initial state, before tipping? I always start with 1 on top, so we have

P^{(0)} = (1, 0, 0, 0, 0, 0).

(The six-tuple notation is the probability of 1-face up, 2-face up, and so on. The superscript on the P indicates how many tips have occurred.)

What are the possibilities for the new top face after the first tip? The 1-face was tipped over, and the 6-face was on the bottom so it can’t be on top. But 2, 3, 4, and 5 are equally likely. So the probability distribution is

P^{(1)} = (0, \frac{1}{4}, \frac{1}{4}, \frac{1}{4}, \frac{1}{4}, 0) .

Note in particular that the up-face after each tip is not independent of the up-face before that tip: given that the current up-face is, say, 3, there are transition probabilities for what the next up-face will be.

What about P^{(2)}, P^{(3)}, and so on? It’s an elementary, but somewhat tedious, exercise in cases to show that

P^{(2)} = (\frac{1}{4}, \frac{1}{8}, \frac{1}{8}, \frac{1}{8}, \frac{1}{8}, \frac{1}{4}) .

One can continue this pattern, but one of my goals is to make an advertisement for the computational advantages of Markov matrices. We’ll find out in section 7.5 how to easily compute P^{(k)}.

Even before computing that, though, I certainly expect that after many tips, the die won’t remember that the starting position was with the 1-face up. So, P^{(k)} ought to begin to look uniform.

Experiment B: The only difference here is that instead of starting with the 1-face up, I instead first roll the die. Then the initial distribution is

P^{(0)} = (\frac{1}{6}, \frac{1}{6}, \frac{1}{6}, \frac{1}{6}, \frac{1}{6}, \frac{1}{6}) .

Again, you can go through the cases and find that

P^{(1)} = (\frac{1}{6}, \frac{1}{6}, \frac{1}{6}, \frac{1}{6}, \frac{1}{6}, \frac{1}{6}) .

and so on. Again, Markov matrices in section 7.5 will make this easier.

In experiment A, the probability distribution started out supported only on 1-up but (we expect) begins to approach uniform; in experiment B, the probability distribution is already uniform and stays that way.

Key points:

• We have a choice of initial state for the die-tipping experiment: for example, one-up, or roll first, or others we could have used.

• We have a probability distribution P^{(0)} for the initial state.

• We have probability distributions P^{(k)} for the subsequent states. These distributions evolve in k.

• We have transition probabilities constraining which face can appear on each step of the experiment. These connect the probability distributions P^{(k)} and P^{(k+1)}. In section 7.5 we will describe the evolution and the transition probabilities in terms of Markov chains.

• Some distributions evolve, or converge, toward a limiting distribution; this is called mixing or thermalization. The limiting distribution (here, uniform) doesn’t evolve. This limiting distribution is called a stationary distribution for the experiment.

• For this experiment, it looks various distributions evolve toward the uniform distribution. We’ll see in section xxx how, and when, this happens.

7.2 Stochastic processes

Definition 7.1. Given a probability space Ω and a sequence P^{(k)} of probability measures for k = 0, 1, 2, …, a discrete-time stochastic process is a sequence of random variables^{1} X^{(k)} from Ω to a state space S.

For concreteness, think of Ω as all N = 2^{n} heads/tails configuration of n coin flips, or all N = 6 results of the landing of a single die, where we allow the probability distributions P^{(k)} to vary in time.

The state space might be another copy of Ω, or ℝ, etc. For example, X^{(k)} : Ω → ℝ might count the number of heads in a set of n coin flips. Or, W^{(k)} : Ω → Ω might be the identity function.

(I should point out that many, if not most, authors define a stochastic process with one “long Ω”, e.g. the probablity space is all sequences of n-coin-flip experiments.)

We categorize stochastic processes by their inputs and outputs:

1For the measure-theoretically inclined reader, X^{(k)} is a measurable function from (Ω, F, P^{(k)}) to (S, G) where F and G are sigma algebras on Ω and S, respectively.

• If k = 0, 1, 2, …, as above, then we have a discrete-time stochastic process, consisting of countably many random variables X^{(k)}. If k ∈ [0, +∞), then we have a continuous-time stochastic process, consisting of uncountably many random variables X^{(k)}.

• If each X^{(k)} has only finitely or countably infinitely many possible output values, then we have a discrete-valued stochastic process, and the distribution of X^{(k)} is a PMF. If each X^{(k)} is a continuous random variable, then we have a continuous-valued stochastic process, and the distribution of X^{(k)} is a PDF.

For this paper, I am considering discrete-time, discrete-valued stochastic processes; Brownian motion is an example of a continuous-time, continuous-valued stochastic process.

7.3 Conditional probability

Forgetting for just a moment about stochastic process, consider two discrete random variables X and Y with joint PMF

f_{X,Y} (x, y) = P (X = x, Y = y).

We can sum over x’s or y’s to obtain the marginal distributions of X and Y:

P (X = x_{i}) = ∑_{j} P (Y = y_{j}, X = x_{i})  and  P (Y = y_{j}) = ∑_{i} P (Y = y_{j}, X = x_{i}).

Recall the definition of conditional probability:

P (Y = y_{j}) = \frac{P (Y = y_{j}, X = x_{i})}{P (X = x_{i})} (if P (X = x_{i}) ≠ 0, otherwise 0).

Multiply both sides by the denominator and sum over i:

P (Y = y_{j}, X = x_{i}) = P (Y = y_{j} | X = x_{i})P (X = x_{i}).

Then we have

P (Y = y_{j}) = ∑_{i} P (Y = y_{j} | X = x_{i})P (X = x_{i}).  (7.1)

For example, take X to be the starting up-face of the die in experiment A or B of section 7.1, and take Y to be the up-face after the first tip. Then the intuition is that equation 7.1 is the following. How likely is it that Y will be in state j? Add up all the likelihoods of all previous states i, weighted by the probability that state i transitions to state j.

In fact, if X and Y take finitely many values, say n each, we can write down an even nicer equation. Take n = 2 for example. Then equation 7.1, for j = 1, 2 simultaneously, is just a matrix product:

(P (Y = y_{1}), P (Y = y_{2})) = (P (X = x_{1}), P (X = x_{2})) ( \begin{matrix}P (Y = y_{1} | X = x_{1}) & P (Y = y_{1} | X = x_{1}) \\ P (Y = y_{1} | X = x_{1}) & P (Y = y_{1} | X = x_{1}) \end{matrix} )

For shorthand, write

P_{i} = P (X = x_{i}),  P′_{j} = P (Y = y_{j}),  and  M_{ij} = P (Y = y_{j} | X = x_{i}).

Then equation 7.1 is

P′_{j} = ∑_{i=1}^{n} P_{i}M_{ij}.  (7.2)

or, more compactly,

P′ = PM.  (7.3)

This matrix M of transition probabilities is called a stochastic matrix, which means the following.

Definition 7.2. A stochastic matrix has all entries between 0 and 1, and all row sums equal to 1.

The reason we choose the row-sum condition is that, given X = x_{i}, it is certain that Y has to go to some y_{j}:

∑_{j=1}^{n} M_{ij} = ∑_{j=1}^{n} P (Y = y_{j} | X = x_{i}) = 1.

7.4 Markov chains

Remember that a (discrete-time) stochastic process is just a sequence of random variables on the same probability space Ω. Those random variables can be related or unrelated. We are particularly interested in the case when the (k + 1)st random variable depends only on the previous one.

Definition 7.3. A Markov chain is a stochastic process X^{(k)} such that

P (X_{n+1} = x | X_{n} = x_{n}, …, X_{1} = x_{1}) = P (X_{n+1} = x | X_{n} = x_{n}).

The intuition is that a Markov process has (at most) only one level of memory. The die-tipping example of section 7.1 is a Markov process. If the up-face at the 19th step is a 2, then the up-face at the 20th step can be 1, 3, 4, or 6. But it doesn’t matter at all how the 19th face got to be a 2.

(I said “at most” since a sequence of independent random variables is also a Markov process: for example, replace the die-tipping experiment with a die-rolling experiment.)

(An example of a non-Markov process is a modified die-tipping experiment wherein I tip to a new face that avoids the previous one. E.g. if the 10th face is a 3 and the 9th face was a 2, I would choose the 11th face with equal likelihood from 1, 2, and 6.)

To apply conditional probability to Markov processes, simply replace X and Y of the previous section with sequence of W^{(k)} which are identity functions from Ω → Ω. The W^{(k)} is simply the random variable which is the state (of the coin flips, die-tips, etc.) at the kth time step. Also replace the P’s of the previous section with a sequence of P^{(k)}’s. Then we have probablity measures

P^{(k+1)}({ω^{(j)}})  and  P^{(k)}({ω^{(i)}}).

on the sample space Ω, we have transition probabilities

M_{ij}^{(k,k+1)} = P (W_{k+1} = ω^{(j)} | W^{(k)} = ω^{(i)}),

and the PMFs at each successive time step are related by matrix multiplication: just as in equation 7.3:

P_{j}^{(k+1)} = ∑_{i} P_{i}^{(k)}M_{ij}^{(k,k+1)}

i.e.

P^{(k+1)} = P^{(k)}M^{(k,k+1)}.  (7.4)

If additionally (as is the case in this paper) the transition probabilities are the same for each k, then we have a time-homogeneous Markov chain. Then the successive PMFs are related by

P^{(k+1)} = P^{(k)}M  (7.5)

and so

P^{k} = P^{(0)}M^{k}.  (7.6)

Thus, it will be interesting below [xxx] to see what happens to powers M^{k} of the Markov matrix M.

For our purposes, think of a Markov chain as a fixed matrix M along with with the sequence P^{(k)} of PMFs, remembering that there’s only the choice of P^{(0)} after which point subsequent PMFs are obtained by matrix multiplication.

When implementing such things in software (the electronics don’t have the error-correcting intuition your brain does), it’s especially important to get the i’s and j’s straight:

• M_{ij} = P (W_{k+1} = ω^{(j)} | W^{(k)} = ω^{(i)}).

• M_{ij} = P (ω^{(i)} ↦ ω^{(j)}).

• When looking for the transition probability from state i to state j, look in row i and column j of M.

7.5 Die-tipping revisited

Now that we know about the evolution of PMFs using Markov matrices, it’s much easier to compute the PMFs for the die-tipping experiments, as promised in section 7.1. In both experiments, the transition probabilities are

M = ( \begin{matrix}0 & \frac{1}{4} & \frac{1}{4} & \frac{1}{4} & \frac{1}{4} & 0 \\ \frac{1}{4} & 0 & \frac{1}{4} & \frac{1}{4} & 0 & \frac{1}{4} \\ \frac{1}{4} & \frac{1}{4} & 0 & 0 & \frac{1}{4} & \frac{1}{4} \\ \frac{1}{4} & \frac{1}{4} & 0 & 0 & \frac{1}{4} & \frac{1}{4} \\ \frac{1}{4} & 0 & \frac{1}{4} & \frac{1}{4} & 0 & \frac{1}{4} \\ 0 & \frac{1}{4} & \frac{1}{4} & \frac{1}{4} & \frac{1}{4} & 0 \end{matrix} )

E.g. if the 1-face is up, afte the die-tip the up-face can’t be 1 or 6 but it can be any of the other four with equal likelihood.

In experiment A, the initial PMF was

P^{(0)} = (1, 0, 0, 0, 0, 0).

Then we can just turn the crank and compute

k  P^{(k)}
0 1.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
1 0.0000000 0.2500000 0.2500000 0.2500000 0.2500000 0.0000000
2 0.2500000 0.1250000 0.1250000 0.1250000 0.1250000 0.2500000
3 0.1250000 0.1875000 0.1875000 0.1875000 0.1875000 0.1250000
4 0.1875000 0.1562500 0.1562500 0.1562500 0.1562500 0.1875000
5 0.1562500 0.1718750 0.1718750 0.1718750 0.1718750 0.1562500
6 0.1718750 0.1640625 0.1640625 0.1640625 0.1640625 0.1718750
7 0.1640625 0.1679688 0.1679688 0.1679688 0.1679688 0.1640625
⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮
19 0.1666660 0.1666670 0.1666670 0.1666670 0.1666670 0.1666660
20 0.1666670 0.1666665 0.1666665 0.1666665 0.1666665 0.1666670

xxx comment, and plots. Likewise, for experiment B,

k  P^{(k)}
0 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667
1 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667
2 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667
⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮

xxx compute higher powers of M.

7.6 Markov-chain properties

These are discussed thoroughly in chapter 4 of [Kerl2010]:

• irreducibility

• aperiodicity

• ergodicity

• reversibility

• stationarity

Numerical examples from handwritten notes. xxx plot them with and without transposing — informative! xxx note that this is ideal for detecting thermalization, but for large N, writing down the entire PMF is precisely what we cannot do. Hence the art of thermalization detection.

8 Markov chain Monte Carlo

8.1 Metropolis choice of Markov matrix

Given P, find M. Various possibilities; “good” convergence rate. Def the particular M. Prove it satisfies DB by construction. Prove irr. Prove aper. Not true for uniform P_{0}! Need something else there …. Note high correlation of X^{(k)}’s.

8.2 Sampling theory

Thm: Sample mean over chain converges to true mean, when M is finite, irr, and aper. Prove this. c/\sqrt{n} convergence rate.

9 Generalizations

xxx to do. Citations.

9.1 2D Ising 9.2 Non-example: percolation 9.3 Numerical integration

xxx continuous here …. xxx emph large number of variables; Riemann mesh impossible.

References

[GS] Grimmett, G. and Stirzaker, D. Probability and Random Processes, 3rd ed. Oxford, 2001.

[Huang] Huang, K. Introduction to Statistical Physics. CRC Press, 2001.

[Kerl] Kerl, J. Probability notes.

http://math.arizona.edu/˜ kerl/doc/prb.pdf

[Kerl2010] Kerl, J. Critical behavior for the model of random spatial permutations (doctoral dissertation) http://math.arizona.edu/˜ kerl/rcm/kerl-dis-gc-ssp.pdf

[Lawler] Lawler, G. Introduction to Stochastic Processes (2nd ed.). Chapman & Hall/CRC, 2006.

[Mackay] MacKay, D.J.C. Introduction to Monte Carlo Methods.

http://www.cs.toronto.edu/˜ mackay/p0.html#BayesMC.html

[NR] Press, W. et al. Numerical Recipes (2nd ed.). Cambridge, 1992.

[Wik] Articles on Ising Model, Statistical Ensemble, and Statistical Mechanics.

http://en.wikipedia.org.

30

Index

B Boltzmann distribution . . . . . . . . . . . . . . . . . . . . . . . . . . .7 Boltzmann’s constant . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

C central limit theorem . . . . . . . . . . . . . . . . . . . . . . . . . . . .10 conditional probability . . . . . . . . . . . . . . . . . . . . . . 21, 23 continuous-time stochastic process . . . . . . . . . . . . . . 23 continuous-valued stochastic process . . . . . . . . . . . . .23 converge . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 coupling coefficients . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 cumulative distribution function . . . . . . . . . . . . . . . . .11

D die-tipping example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 discrete-time stochastic process . . . . . . . . . . . . . . 22, 23 discrete-valued stochastic process . . . . . . . . . . . . . . . .23

E energy function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

F flips . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

G given . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .21

I initial state . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 inverse temperature . . . . . . . . . . . . . . . . . . . . . . . . . . 7, 13 Ising model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

J joint PMF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

L law of large numbers . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 limiting distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . .22

M magnetic field term . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .7 marginal distributions . . . . . . . . . . . . . . . . . . . . . . . . . . .23 Markov chain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 Markov chains . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .21, 22 mixing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

P partition function . . . . . . . . . . . . . . . . . . . . . . . . . . . . .7, 12 PDF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 PMF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .21 probability density function . . . . . . . . . . . . . . . . . . . . . 21 probability distribution . . . . . . . . . . . . . . . . . . . . . . . . . 22 probability mass function . . . . . . . . . . . . . . . . . . . . . . . 21 probability measure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 probability numerator . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

R

random variable . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .21

S

sample space . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 spin configuration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 spins . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 state space . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 stationary distribution . . . . . . . . . . . . . . . . . . . . . . . . . . 22 step . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .13 stochastic matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .24 stochastic process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 stochastic processes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 subsequent states . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 sweep . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

T

thermalization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 time-homogeneous Markov chain . . . . . . . . . . . . . . . . 25 transition probabilities . . . . . . . . . . . . . . . . . . . 21, 22, 24 transition probability . . . . . . . . . . . . . . . . . . . . . . . . 15, 23