F.Beukers, April 15, 2004
Notebook 9 - The Ehrenfest urn model

Background

First we recall some theory. Let M be a square matrix. Suppose that it has a real eigenvalue (Subscript[λ, 0])which is larger or equal to the absolute values of all other eigenvalues. We call (Subscript[λ, 0]) dominant if all other n - 1 eigenvalues have absolute value strictly less than (Subscript[λ, 0]). Suppose that (Subscript[λ, 0]) is dominant and let (Subscript[v, eigen])be the corresponding eigenvector (up to a scalar vector). Then, for almost all starting vectors (Subscript[v, 0]) we have that (((Subscript[λ, 0]))^(-n)) (M^n)((Subscript[v, 0])) goes to a scalar multiple of (Subscript[v, eigen]) as n goes to ∞. 

Suppose in addition that M is a probability matrix corresponding to a Markov chain. Then M has the eigenvalue 1 and all other eigenvalues have absolute value ≤1. If 1 is a dominant eigenvalue, the behaviour of the Markov chain is simple, for almost all choices of (Subscript[v, 0]) the vector (M^n)((Subscript[v, 0 ])) goes to the eigenvector (Subscript[v, eigen]) as n goes to ∞. Here, (Subscript[v, eigen]) is chosen so that the sum of the components is 1, i.e. it represents the final probability distribution.

There is a nice criterion to see if 1 is a dominant eigenvalue. Draw the graph corresponding to the Markov chain and check if it is irreducible and a-periodic. If that is the case, 1 is dominant. However, there are interesting Markov chains that do not have a dominant eigenvalue.

Exercise 1. Draw the transition graph of the coin game between Abe and Bert of the previous Mathematica Notebook. Check its irreducibility and a-periodicity.

The Ehrenfest Urn Model

Another Markov chain where 1 is not dominant is the Eherenfest urn model, named after the famous Austrian physicist Paul Eherenfest (1880-1933).  Suppose we have a gas container consisting of two compartments which are connected by a hole.

Suppose our gas consists of N molecules, which are initially in the right hand compartment. As the molecules bounce around, a number of them may enter the left hand compartment. Experience tells us that for large N the molecules eventually spread evenly over both compartments. The famous physicist Ehrenfest has proposed a simple model to mimick this behaviour.

Time is divided into time steps. At each time step we randomly select a molecule and take it to the other compartment. We suppose that the probability for a molecule to be chosen is 1/N. In this way we get a Markov chain whose states may be characterised by the number of molecules in the left hand compartment. So we have N + 1 states.

Exercise 2. In the case when N = 4, draw a transition graph of the Markov chain and determine the probabilty matrix M. (Check with the staff whether you have found the correct matrix, since you are going to work with it).
Determine the eigenvalues of M. Note that  1 is not a dominant eigenvalue.

Exercise 3. Let (Subscript[v, 0]) be the intial vector representing the state with no molecules in the left compartment. Determine (M^n)((Subscript[v, 0])) for a couple of large consecutive values of n and convince yourself that there is no convergence. Give an explanation for this phenomenon.

A Revised Model

To make the Eherenfest model more realistic, we give it an extra twist. Every time we have selected a molecule, we toss a fair coin to decide if we move it to the other compartment.

Exercise 4. Write down the transition graph and transition matrix M for the new model in the case N = 4. Determine the eigenvalues. Is 1 a dominant eigenvalue?

Exercise 5. Check that the new transition matrix is actually the average of the old transition matrix and the identity matrix. Can you explain the relation between the old eigenvalues and the new ones?

Exercise 6. Repeat Exercise 3 with the new matrix and convince yourself numerically that we do have convergence.

Now that we have a model which seems to work, let us look at its prediction for (somewhat) large numbers of molecules. First we construct a function giving the transition matrix for n molecules

urn[n_] := 
      Module[{m, k}, 
m = 
          Table[0, {n + 1}, {n + 
                1}]; 
Do[
m[[k + 1, 
                k]] = 
            1 - (k - 1)/n; 
m[[k, k + 1]] = 
            k/n, 
{k, 1, 
            n}
]; 
m = (m + 
                IdentityMatrix[n + 1])/2
]

Exercise 7. Compute the eigenvalues of the transition matrix for a number of values of N below 20. Make a guess what eigenvalues you will get for any value of N.
Can you also make a guess about the shape of the eigenvector corresponding to the eigenvalue 1? To get an idea, compare the eigenvectors for N=3,4,5,6.

Exercise 8. Perform the following steps for N = 50, 100, 200. Construct the initial state vector corresponding to zero molecules  in the left compartment. Apply the Markov matrix to this state as many times as you think is needed to reach equilibrium. Make a ListPlot of the resulting limit vector.