Mercurial > hg > mauch-mirex-2010
annotate _FullBNT/KPMstats/mc_stat_distrib.m @ 9:4ea6619cb3f5 tip
removed log files
author | matthiasm |
---|---|
date | Fri, 11 Apr 2014 15:55:11 +0100 |
parents | b5b38998ef3b |
children |
rev | line source |
---|---|
matthiasm@8 | 1 function pi = mc_stat_distrib(P) |
matthiasm@8 | 2 % MC_STAT_DISTRIB Compute stationary distribution of a Markov chain |
matthiasm@8 | 3 % function pi = mc_stat_distrib(P) |
matthiasm@8 | 4 % |
matthiasm@8 | 5 % Each row of P should sum to one; pi is a column vector |
matthiasm@8 | 6 |
matthiasm@8 | 7 % Kevin Murphy, 16 Feb 2003 |
matthiasm@8 | 8 |
matthiasm@8 | 9 % The stationary distribution pi satisfies pi P = pi |
matthiasm@8 | 10 % subject to sum_i pi(i) = 1, 0 <= pi(i) <= 1 |
matthiasm@8 | 11 % Hence |
matthiasm@8 | 12 % (P' 0n (pi = (pi |
matthiasm@8 | 13 % 1n 0) 1) 1) |
matthiasm@8 | 14 % or P2 pi2 = pi2. |
matthiasm@8 | 15 % Naively we can solve this using (P2 - I(n+1)) pi2 = 0(n+1) |
matthiasm@8 | 16 % or P3 pi2 = 0(n+1), i.e., pi2 = P3 \ zeros(n+1,1) |
matthiasm@8 | 17 % but this is singular (because of the sum-to-one constraint). |
matthiasm@8 | 18 % Hence we replace the last row of P' with 1s instead of appending ones to create P2, |
matthiasm@8 | 19 % and similarly for pi. |
matthiasm@8 | 20 |
matthiasm@8 | 21 n = length(P); |
matthiasm@8 | 22 P4 = P'-eye(n); |
matthiasm@8 | 23 P4(end,:) = 1; |
matthiasm@8 | 24 pi = P4 \ [zeros(n-1,1);1]; |
matthiasm@8 | 25 |
matthiasm@8 | 26 |