Coupling+of+Markov+chains+-+Perfect+sampling


 * Tutorial 01: Coupling of Markov chains - Perfect sampling **

toc 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 ) 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.

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 //t// coup, 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 //t// coup 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 //π//(//t// coup ) of //t// coup 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).



=Algorithms=

Forward coupling
We consider below the forward coupling for the generic //N// sites problem depicted above. The jump probabilities are math p_{k\to k\pm 1}=p_{k\to k}=\frac 13 math excepted at the boundaries where math 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 math This choice implements reflecting boundary conditions. Since the rates are symmetric math p_{k\to k'}=p_{k'\to k} math detailed balance is satisfied with the uniform distribution probability math \pi_k=\frac 1N math To obtain a histogram of the site occupation at coupling time, one may run the following Python program: code format="python" 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 code 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: code format="python" 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 code 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:

code format="python" 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

code

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: math \delta^{(n)}_{ij} = \sum_{k} |p^{(n)}_{ik} - p^{(n)}_{jk}| math where //p// (//n//) is the //n//-fold iterated transfer matrix. Weakly ergodic means that math \forall i,j\,,\quad \lim_{n \to \infty}\delta^{(n)}_{ij} \to 0 math The use of coupling arguments to prove weak ergodicity seems to go back to W. Doeblin, in 1937 These arguments were thus used for decades before Propp and Wilson were able to show its practical (and revolutionary) use...

url}?f=print|print this page