MATH 384
13 February 2024

Example

Working by hand, we found the following matrix P. Now let’s confirm that it has the desired eigenvector with eigenvalue 1.

p={{1/3,1/6,1/9},{1/3,1/2,2/9},{1/3,1/3,2/3}}
MatrixForm[p]

Eigensystem[p]

Finding Markov chain transition matrices

Here is a module that constructs a transition matrix given a desired steady-state probability distribution v and a proposal matrix Q. This implements Algorithm 6.5 in the text.

(* inputs: v is the steady-state vector, q is the proposal transition matrix *)
buildTransitionMatrix[v_,q_]:=Module[
(* list of local variables within the module *)
{n=Length[v] (* dimension of the vector *),
p (* the transition matrix *),
tot (* used for adding up each column*)},

(* initialize p to a n-by-n zero matrix *)
p=Table[0,n,n];

(* loop over each column *) 
Do[
(* set total to zero *)
tot=0;

(* loop over each row *) 
Do[
(* we only set off-diagonal entries in this loop*)
If[row!=col,
(* compare steady-state frequencies and set p[[row,col]] accordingly *)
If[v[[col]]<v[[row]],
p[[row,col]]=q[[row,col]],
p[[row,col]]=((v[[row]])/(v[[col]]))q[[row,col]]
];
(* keep track of the sum of the entries in this column *) 
tot+=p[[row,col]]
]
,{row,1,n}];

(* the last entry in the column to set is the diagonal entry *) 
p[[col,col]]=1-tot
,{col,1,n}];

(* all done! return matrix p *) 
p
]

Example 1: Let’s recompute the matrix from class with v=((1/6),(1/3),(1/2)):

v={1/6,1/3,1/2};
q=Table[1/3,3,3];
p=buildTransitionMatrix[v,q];
MatrixForm[p]

Example 2: Find a matrix with steady-state vector v=((3/10),(1/2),(1/5)):

v={3/10,1/2,1/5};
q=Table[1/3,3,3];
p=buildTransitionMatrix[v,q];
MatrixForm[p]

Eigensystem[p]

Example 3: Find another matrix with steady-state vector v=((3/10),(1/2),(1/5)):

v={3/10,1/2,1/5};
q={{2/3,1/6,1/6},{1/6,2/3,1/6},{1/6,1/6,2/3}}
p=buildTransitionMatrix[v,q];
MatrixForm[p]

Eigensystem[p]

Example 3: How about another matrix with steady-state vector v=((3/10),(1/2),(1/5)):

v={3/10,1/2,1/5};
q={{1/2,1/4,1/4},{1/4,1/2,1/4},{1/4,1/4,1/2}}
p=buildTransitionMatrix[v,q];
MatrixForm[p]

Eigensystem[p]

Simulating a Markov chain without a transition matrix

We will generate a Markov chain of a specified length, given a desired steady-state distribution and a proposal transition function. We will not store a matrix in memory. This implements Algorithm 6.8 in the text.

For this example, our states will be the integers 1,2,3,4,5,6,7,8,9,10.

Our proposal transitions will be as follows: from any state k, we propose a transition to k-1, k, or k+1, each with probability (1/3). However, we must handle the endpoints 1 and 10 specially. If k=1, then we propose to stay put with probability (2/3) and move to 2 with probability (1/3). If k=10, then we propose to stay put with probability (2/3) and move to 9 with probability (1/3). (See Figure 6.4 in the text for a diagram.)

First, store the steady state distribution as a function. Specifically, the function v(k) will give the relative frequency of state k.

Clear[v]
v[k_]:=k

Plot[v[k],{k,1,10}]

Next, define parameters and initialize a list to store the sequence of states.

length=1000000; (* number of steps to simulate *)
start=4; (* starting state *)
chain=Table[0,length]; (* initialize a list to store the sequences of states *)
chain[[1]]=start; (* store the starting state in the list *)

Now we can simulate the Markov chain:

(* repeat for m from 2 to length *)
Do[
(* first, select a proposed transition *)
curr=chain[[m-1]];
prop=RandomChoice[{curr-1,curr,curr+1}]; (* choose with equal probability *) 
If[prop<1,prop=1]; (* don't go off the left end *) 
If[prop>10,prop=10]; (* don't go off the right end *) 

(* second, decide whether or not to move to the proposed state *)
ratio=((v[prop])/(v[curr])); (* this is the acceptance ratio *) 
r=RandomReal[];  (* generate a real number uniformly between 0 and 1 *) 
(* if v[prop]<v[curr], move with probability given by ratio; otherwise stay put *) 
If[r<ratio,chain[[m]]=prop,chain[[m]]=curr] 

,{m,2,length}]

chain

Histogram[chain]

The histogram matches the desired long-term frequencies!