Module 3 : Finite Mathematics

304 : Markov Processes

O B J E C T I V E

We will construct transition matrices and Markov chains, automate the transition process, solve for equilibrium vectors, and see what happens visually as an initial vector transitions to new states, and ultimately converges to an equilibrium point.

S E T U P

In this project we will use the following command packages. Type and execute this line before begining the project below. If you re-enter the worksheet for this project, be sure to re-execute this statement before jumping to any point in the worksheet.

> restart; with(linalg): with(plots): with(plottools):

Warning, the protected names norm and trace have been redefined and unprotected

Warning, the name changecoords has been redefined

___________________________________________________________________________________

A. Making The Transition

___________________________________________________________________________________

First, we will simply construct a transition matrix for a Markov Process and later use it to create a Markov Chain. Suppose we begin with the situation where all of the students in a class are earning grades of A, B, or C and the teacher does not believe in giving grades of D or F. On the most recent exam, 28% earned As, 32% earned Bs, and 50% earned Cs. We can express this break down with a current state vector.

\

> current_state := vector( [ .28, .32, .5]);

Note that the entries in the vector add to 1 since these three states represent the entire class.

Lets further assume that in this class, a person who earns an A on one exam has a 48% of earning an A on the next exam, a 37% of earning a B on the next exam, and a 15% chance of earning a C on the next exam. We can indicate these probabilities by a transition vector from state A to states A, B, and C.

Vector a represents the probabilities of an A student earning an A, B, or C on the next exam.

> a := vector( [.48, .37, .15]);

Note that the sum of the entries is 1 because the student must earn an A, B, or C. There is a 100% probability that the grade will be one of these because no Ds or Fs are given.

In a similar way, lets assume that the probability that a student who earns a B on an exam will get an A, B, or C are 30%, 45%, and 25% respectively, and the probability of a C student earning an A, B, or C on the next exam are 10%, 30%, 60% respectively. Thus, the transition vectors for B and C students are as follows.

Vector b represents the probabilities of a B student earning an A, B, or C on the next exam, and vector c holds the probabilities of a C student earning an A, B, or C. Again note that the sum of the entries of each is 1.

> b := vector( [.30, .45, .25]);

> c := vector( [.10, .30, .60]);

All of the information about the transition of grades from one exam to the next is encoded in these three vectors. We can now form a transition matrix with these transition vectors as the rows. Alternatively, you can also make a direct definition of the transition matrix in this way.

> A := matrix( [a,b,c]);

> A := matrix( [ [.48, .37, .15], [.30, .45, .25], [.10, .30, .60] ]);

The direct method of defining a transition matrix is to list the transition vectors in order.

And now the big question : What will be the proportion of grades on the next exam?

> current_state := vector( [ .28, .32, .5]);
next_state := evalm(current_state &* A);

We multiply the vector times the matrix, and simplify to find the proportion of each grade on the next test.

What about the next exam, and the next? This process of transitioning to each new result create a Markov chain. Here is a sequence of 12 transitions in the Markov Chain for this scenario.

> for k from 1 to 12 do print(k), evalm( current_state &* A^k); od;

___________________________________________________________________________________

B. The Equilibrium Vector

___________________________________________________________________________________

As we saw above, the transition states do not necessarily fluctuate wildly. In fact, they often converge to particular values. The ultimate transition, is called the equilibrium vector. It is, in a sense, the limit of the transitions.

Using pencil and paper you can solve a matrix to find the equilibrium vector using this process. Write the matrix equation vA = v, which can be transformed into vA - v = 0. When you simplify this expression you get a number of linear equations. When you attempt to solve this system you will find that the system is dependent and can not be solved completely as it is. However, using the fact that the solution must be a probability vector, that is, the sum of the components is 1, its usually possible to solve the system and find a unique solution. We follow a similar path using Maple to find an equilibrium solution.

SLOVE A 2 X 2

Here is an example of a transition matrix for two states and an indefinite vector which we will solve to find the equilibrium vector.

> A := matrix( [[ .8,.2],[.35,.65]]);
v := vector( [s,t]);

Now we will find the equilibrium vector by simplifying, converting to linear equations, and solving.

> evalm( v &* A - v); <-------- Simplify the expression (vA - v)
convert(%, set); <-------- Convert it into a set of linear expression
subs( {t = 1-s}, %); <-------- Eliminate t since s + t = 1
solve( %[1] = 0); <-------- Solve the resulting equation
equivect := [ %, 1-%]; <-------- Form the solution

> solve_markov := proc( A )
local s,t,x;
x := vector( [s,t]);
convert(evalm( x &* A - x), set);
solve( subs( {t = 1-s}, %)[1] = 0);
[ %, 1 -%];
end:

AUTOMATICALLY SOLVING 2 X 2 's

Here is a custom function which will solve a two state system automatically. Enter the following commands and execute. Nothing will happen, but the definition for this new function will have been created. You can then access it.

> A := matrix( [[ .8,.2],[.35,.65]]);
solve_markov( A );

SOLVE A 3 X 3

Here is a similar process for solving a three state system to what we used above for the two state system.

> A := matrix([[.90,.07,.03],[.25, .37,.38],[.6,.1,.3]]);

> v := vector([a,b,c]);

> evalm( v &* A - v);
convert(%, set);
subs( {c = 1-a-b}, %);
solve( { %[1], %[2], c=1-a-b}, {a,b,c});

___________________________________________________________________________________

C. Geometric View of Markov Convergence

___________________________________________________________________________________

Its quite illuminating to see what happens when the current state transitions to various new states along with the equilibrium vector.

2 DIMENSIONAL MARKOV CHAINS

Here is a special procedure which will show all of that.

> markov_plot := proc( A, v)
local r, u, vp; r:= .03;
vp := convert(v, list); u := convert(solve_markov( A ), list);
display( {disk( vp, 1.5*r, color=yellow),
textplot( [vp [1]+ 2*r, vp[2],`Initial Proportion`],align={ABOVE,RIGHT}),
seq( disk( convert( evalm(v &* A^k ), list), r,color=blue), k=1..5),
disk( u ,1.5*r, color=red),
textplot([u[1]+ 2*r, u[2], `Equilibrium Vector`], align={ABOVE,RIGHT}),
plot( 1-x, x = 0..1, color = green, linestyle = 2) }, scaling = constrained );
end:

When you execute the command block above, nothing will happen. However, this new command is defined and ready to use. Lets put it to use at once.

Note that the initial points are different, but the equilibrium points are the same. This is due to the fact that the equilibrium vector depends only on the matrix and not the initial vector.

> A := matrix( [[ .15,.85],[.44,.54]]);
v := vector( [.05,.95]);
markov_plot( A, v);
v := vector( [.83,.17]);
markov_plot( A, v);

3 DIMENSIONAL MARKOV CHAINS

We can also visualize a three dimensional probability vector and its Markov transitions.

> A := matrix( [[.15,.65,.2],[.4,.35,.25],[.8, 0,.2]]);
v := vector( [ .31, .32, .41]);

> display( pointplot3d( convert( v, list), color = red, symbol = box ,axes = normal),
seq( pointplot3d( convert( evalm(v &* A^k), list) , color = blue, symbol = box,
axes = normal) , k =2..8),
plot3d( 1-x-y, x =0..1, y =0..1, view =0..1, style = patchnogrid, orientation=[35,65]) );

If you click on the graph and drag the mouse you can rotate the graph to inspect it from any view point. Note that all of the points lie on the plane x + y + z = 1.

>