% MATLAB Demo (Positive Markov matrices)
% File: markov
%
% Let's experiment with the positive Markov matrix
%
%
%      [ 0.84  0.04  0.04 ]                          
%  A = [ 0.12  0.52  0.12 ]  and an initial vector u0. 
%      [ 0.04  0.44  0.84 ]                          
%
% Each entry in a positive Markov matrix is positive,
% and the entries in each column add up to 1.
% 
% Positive Markov matrices have the additional property
% that one of the eigenvalues is 1, and the other 
% eigenvalues have magnitude less than 1.
% (It is possible for some of the smaller eigenvalues 
% to be complex!)
%
% The objective is to investigate the problem of finding
% the steady-state vector associated with a positive
% Markov matrix.
%
>> A = [0.84 0.04 0.04; 0.12 0.52 0.12; 0.04 0.44 0.84]
A =
    0.8400    0.0400    0.0400
    0.1200    0.5200    0.1200
    0.0400    0.4400    0.8400
%
%%% u0 is an initial probability vector.
%%% Its components are probabilities that add up to 1.
%
>> u0 = [0.5; 0.3; 0.2]
u0 =
    0.5000
    0.3000
    0.2000
%
%%% Let's examine powers of the matrix A
%
>> A5 = A^5
A5 =
    0.4621    0.1345    0.1345
    0.1980    0.2082    0.1980
    0.3399    0.6573    0.6676

>> A10 = A^10
A10 =
    0.2859    0.1785    0.1785
    0.2000    0.2001    0.2000
    0.5141    0.6214    0.6215
%
%%% Let's examine A * u0, A^5 * u0, and A^10 * u0.
%
>> u1 = A * u0
ans =
    0.4400
    0.2400
    0.3200

>> u5 = A5 * u0
ans =
    0.2983
    0.2010
    0.5007

>> u10 = A10 * u0
ans =
    0.2322
    0.2000
    0.5678
%
%%% Note that u1, u5, u10 are also probability vectors.
%%% Note A5 and A10 are also positive, Markov matrices.
%
>> sum(A5)
ans =
    1.0000    1.0000    1.0000

>> sum(A10)
ans =
    1.0000    1.0000    1.0000
%
%%% After k=100 iterations, A^k converges to a limit. 
%
>> A100 = A^100
A100 =
    0.2000    0.2000    0.2000
    0.2000    0.2000    0.2000
    0.6000    0.6000    0.6000
%
%%% Each column of the limit matrix is the same.
%%% Each column is the eigenvector associated with the
%%% eigenvalue lambda = 1.
%%% (The eigenvector is scaled so that its components add to 1.)
%
>> [v, d] = eig(A)
v =
   -0.7071    0.3015    0.0000
    0.0000    0.3015   -0.7071
    0.7071    0.9045    0.7071
d =
    0.8000         0         0
         0    1.0000         0
         0         0    0.4000

%
%%% The eigenvector (associated with lambda = 1) is in the
%%% second column of v.
%
>> x1 = v(:,2)
x1 =
    0.3015
    0.3015
    0.9045
%
%%% Scale the eigenvector so that the components add to 1.
%
>> uinfty = x1 / sum(x1)
uinfty =
    0.2000
    0.2000
    0.6000
%
%%% We can also find the steady-state vector for a positive Markov
%%% matrix in a more direct manner:
%%% Find the special solution to (A - 1*I) x = "0" since it is the
%%% eigenvector corresponding to lambda = 1.
%%% Scale this eigenvector so that its components add to 1.
%
>> I = eye(3)
I =
     1     0     0
     0     1     0
     0     0     1

>> nullbasis(A-eye(3))
ans =
    0.3333
    0.3333
    1.0000

>> ans / sum(ans)
ans =
    0.2000
    0.2000
    0.6000

>> quit