matthiasm@8: function pi = mc_stat_distrib(P)
matthiasm@8: % MC_STAT_DISTRIB Compute stationary distribution of a Markov chain
matthiasm@8: % function pi = mc_stat_distrib(P)
matthiasm@8: % 
matthiasm@8: % Each row of P should sum to one; pi is a column vector
matthiasm@8: 
matthiasm@8: % Kevin Murphy, 16 Feb 2003
matthiasm@8: 
matthiasm@8: % The stationary distribution pi satisfies pi P = pi
matthiasm@8: % subject to sum_i pi(i) = 1,  0 <= pi(i) <= 1
matthiasm@8: % Hence
matthiasm@8: % (P'  0n   (pi  = (pi 
matthiasm@8: %  1n  0)    1)     1)
matthiasm@8: % or P2 pi2 = pi2.
matthiasm@8: % Naively we can solve this using (P2 - I(n+1)) pi2 = 0(n+1)
matthiasm@8: % or P3 pi2 = 0(n+1), i.e., pi2 = P3 \ zeros(n+1,1)
matthiasm@8: % but this is singular (because of the sum-to-one constraint).
matthiasm@8: % Hence we replace the last row of P' with 1s instead of appending ones to create P2, 
matthiasm@8: % and similarly for pi.
matthiasm@8: 
matthiasm@8: n = length(P);
matthiasm@8: P4 = P'-eye(n);
matthiasm@8: P4(end,:) = 1;
matthiasm@8: pi = P4 \ [zeros(n-1,1);1];
matthiasm@8: 
matthiasm@8: