# Coupling of Markov chains - Perfect sampling

Tutorial 01: Coupling of Markov chains - Perfect sampling

The big problem in a Monte Carlo simulation is to determine the convergence time, or at least to estimate it correctly. In many applications, this is a complex, controversial problem, often a nightmare. Here we study the subject of perfect simulations that completely solve this problem, whenever they are applicable.

# Forward Coupling of Markov chains

## Random maps

Consider a complete Monte Carlo simulation (see below figure, taken from[1] ) where a random map prescribes the motion (up, down, straight) of a particle on a lattice of N=5 sites, for every position of the particle at each time step.

 A simulation of a particle moving on the 1D lattice of 5 sites, using random maps.
Transition probabilities are as follows: on site 1 (site 5), arrows go straight with probability 2/3, and up (down) with probability 1/3. On sites 2,3,4, arrows go into the three possible directions with probability 1/3 each.

## Computation of coupling time (example)

Use an (electronic) pencil to determine the coupling time tcoup, after which the Monte Carlo simulation (as realized in the above figure) no longer depends on the initial condition.

## Probablity distribution of coupling time

The coupling time tcoup is not a fixed number, but a random variable, which depends on the way we have drawn the arrows in the figure. Show (without any calculation) that the distribution π(tcoup) of tcoup decays at least exponentially for large coupling times.

## Exponential decay of coupling-time distribution

Use the reasoning of the previous question to show that the distribution decays exactly exponentially (without any calculation).

## Bias in forward coupling

Show (use the algorithm given below) that the position at which the coupling takes place (from 1 to 5) is not uniformly distributed. Comment.

# Backward Coupling of Markov chains

## Coupling from the past: perfect sampling

Below, please find an infinite Monte Carlo calculation (not all arrows have been drawn). Use pencil and eraser to find the position at t = 0 for any of the initial conditions at t=-. Show that this position is uniformly distributed (if we average over the ways of drawing the arrows).

 A simulation on a one-dimensional lattice that has run for more than one billion years.

# Algorithms

### Forward coupling

We consider below the forward coupling for the generic N sites problem depicted above. The jump probabilities are
$p_{k\to k\pm 1}=p_{k\to k}=\frac 13$
excepted at the boundaries where
$p_{1\to 1}=p_{N\to N}=\frac 23 \quad\text{and}\quad p_{1\to 2}=p_{N\to N-1} = \frac 13$
This choice implements reflecting boundary conditions. Since the rates are symmetric
$p_{k\to k'}=p_{k'\to k}$
detailed balance is satisfied with the uniform distribution probability
$\pi_k=\frac 1N$
To obtain a histogram of the site occupation at coupling time, one may run the following Python program:
import random, pylab

N = 8
pos = []
for statistic in range(5000):
positions = set(range(N))
while True:
arrows = [random.randint(-1,1) for i in range(N)]
if arrows[0] == -1: arrows[0] = 0
if arrows[N-1] == 1: arrows[N - 1] = 0
positions = set([b + arrows[b] for b in positions.copy()])
if len(positions) == 1: break
pos.append(positions.pop());
pylab.title('Forward coupling: distribution of the coupled configuration')
pylab.hist(pos, bins = N, range = (-0.5, N - 0.5), normed = True)
pylab.xlabel('Position')
pylab.ylabel('Histogram')
pylab.show()
One checks that the histogram is not flat. It does not yield the uniform steady state πk = 1/N :

### Exact sampling from backward coupling

We consider the same system, now using the backward sampling method: starting from t = 0, we go backwards in time, constructing the graph starting from decreasing initial times t = -1, t = -2, ... , and waiting until the graph at t = 0 contains only one branch.
Note that, in contrast to the forward coupling method, the coupling may have occurred at any negative time t < 0, and not at the observation time t = 0.
To histogram the site occupation at t = 0, one may run the following Python program:
import random, pylab

N = 8
pos = []
times_tot = []
for statistic in range(5000):
all_arrows = {}
time_tot = 0
while True:
time_tot -= 1
arrows = [random.randint(-1, 1) for i in range(N)]
if arrows[0] == -1: arrows[0] = 0
if arrows[N-1] == 1: arrows[N-1] = 0
all_arrows[time_tot] = arrows
positions = set(range(N))
for t in range(time_tot, 0):
positions = set([b + all_arrows[t][b] for b in positions.copy()])
if len(positions) == 1: break
times_tot.append(-time_tot) # here we store the coupling times
pos.append(positions.pop())
pylab.title('Backward coupling for a particule diffusion on [1,N]: histogram of the position at t=0')
pylab.hist(pos, bins = N, range = (-0.5,N - 0.5), normed = True)
pylab.savefig('backward_coupling_position.png')
pylab.xlabel('Position')
pylab.ylabel('Histogram')
pylab.show()
This coupling-from-the-past method for sampling is perfect: configurations at time t = 0 are now exactly distributed according to the uniform steady state distribution:

### Relation to exponential convergence

To see the relation between the coupling time and the exponential convergence seen in the transfer matrix approach, we can record the distribution of the coupling times in above algorithm. Here are the lines you have to add at the end of the previous algorithm so as to perform the analysis of this distribution:

pylab.title('Backward coupling: distribution of coupling times')
pylab.hist(times_tot, bins = 150, range=(1,150), normed = True)
pylab.savefig('backward_coupling_times.png')
pylab.show()

pylab.title('Backward coupling: cumulative distribution of coupling times')
pylab.hist(times_tot, bins = 150, normed = True, histtype = 'step', cumulative = -1)
pylab.savefig('backward_coupling_times_cumulative.png')
pylab.show()

pylab.title('Backward coupling: cumulative distribution of coupling times in log scale')
pylab.hist(times_tot, bins = 150, range=(1,150), normed = True, histtype = 'step', cumulative = -1)
pylab.savefig('backward_coupling_times_cumulative_log-scale.png')
pylab.semilogy()
pylab.axis([-1, 151, 0.001, 1.1])
pylab.show()


The following histogram shows the distribution of coupling times. Integrating this distribution over time, we can obtain the number of states that have perfectly coupled after a given time t. When we then plot the ratio uncoupled states over the total number of states, we see that it decays exponentially. To reconcile the transfer matrix picture with what we learned here about perfect sampling, we can say, that the exponentially decaying error (that we've seen, e.g., in the transfer matrix approach to the pebble game) stems from an exponentially decaying number of imperfect samples.

Compare numerically the distributions of the backward and of the forward coupling times. Explain.

# Historical remark: how coupling arose in mathematics

In mathematics, one defines the concept of weak ergodicity as follows:
$\delta^{(n)}_{ij} = \sum_{k} |p^{(n)}_{ik} - p^{(n)}_{jk}|$
where p(n) is the n-fold iterated transfer matrix. Weakly ergodic means that
$\forall i,j\,,\quad \lim_{n \to \infty}\delta^{(n)}_{ij} \to 0$
The use of coupling arguments to prove weak ergodicity seems to go back to W. Doeblin, in 1937[2] These arguments were thus used for decades before Propp and Wilson[3] were able to show its practical (and revolutionary) use...