Home > wafo > markov > mc2stat.m

mc2stat

PURPOSE ^

Calculates the stationary distribution for a Markov chain.

SYNOPSIS ^

[ro,PP]=mc2stat(P)

DESCRIPTION ^

  MC2STAT  Calculates the stationary distribution for a Markov chain.
 
  CALL: ro = mc2stat(P);
 
  ro = Stationary distribution.   [1xr]
 
  P  = transition matrix.         [rxr]
 
  Example: 
    F = magic(5)
    P = mat2tmat(F)
    ro = mc2stat(P)
 
  See also  mat2tmat, mc2reverse, smc2stat, mctp2stat

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function [ro,PP]=mc2stat(P)
002 % MC2STAT  Calculates the stationary distribution for a Markov chain.
003 %
004 % CALL: ro = mc2stat(P);
005 %
006 % ro = Stationary distribution.   [1xr]
007 %
008 % P  = transition matrix.         [rxr]
009 %
010 % Example: 
011 %   F = magic(5)
012 %   P = mat2tmat(F)
013 %   ro = mc2stat(P)
014 %
015 % See also  mat2tmat, mc2reverse, smc2stat, mctp2stat
016   
017 % Tested  on Matlab  5.3
018 %
019 % History:
020 % Revised by PJ  23-Nov-1999
021 %   updated for WAFO
022 % Created by PJ (Pär Johannesson) 1998
023 %   from 'Toolbox: Rainflow Cycles for Switching Processes V.1.0'
024 %   previous name: 'statf_mc'
025 
026 % Check input arguments
027 
028 ni = nargin;
029 no = nargout;
030 error(nargchk(1,1,ni));
031 
032 n=length(P); % number of states
033 
034 PP=P;
035 
036 % Set negative elements to zero
037 [I,J] = find(PP<0);
038 if length(I) ~= 0
039   warning(['Warning: Negative elements in P. Setting to zero!']);
040   for k = 1:length(I)
041     PP(I(k),J(k)) = 0;
042   end
043 end
044 
045 J =[];
046 NoDeleted = 0; NoDeleted0 = -1; 
047 while NoDeleted ~= NoDeleted0
048 
049   % Check that the rowsums of PP are equal to 1
050 
051   sumPP = sum(PP');
052   if sum(abs(sumPP-1) <= n*eps) ~= length(PP)  % Number of rowsums equal to 1
053     warning(['Warning: Rowsums of P not equal to 1. Renormalizing!']);
054     for i = 1:length(PP)
055       if sumPP(i) == 0
056         if isempty( find(J==i) )  % Not allready deleted
057           warning(['Warning: Rowsum ' num2str(i) ' of P equals zero. Deleting state ' num2str(i) '!']);
058           J = [J i];
059           PP(:,i) = zeros(n,1);
060         end
061       else
062         PP(i,:) = PP(i,:)/sumPP(i);
063       end
064     end
065   end
066 
067   NoDeleted0 = NoDeleted;
068   NoDeleted = length(J);
069 
070 end
071 
072 I=1:n; I(J)=[]; 
073 P1 = PP(I,I);  % Delete states
074 
075 
076 % Compute stationary distribution
077 
078 r=size(P1,1);
079 P1=P1';
080 M=[eye(r-1) zeros(r-1,1)];
081 P1=P1(1:r-1,:);
082 P1=P1-M;
083 P1=[P1 ; ones(1,r)];
084 a=[zeros(r-1,1); 1];
085 roP1=P1\a;
086 %roP1=pinv(P1)*a;
087 roP1=roP1';
088 
089 % Set ro(i)=0 for 'singular states'
090 
091 ro = zeros(1,n);
092 ro(I)=roP1;
093 
094

Mathematical Statistics
Centre for Mathematical Sciences
Lund University with Lund Institute of Technology

Comments or corrections to the WAFO group


Generated on Thu 06-Oct-2005 02:21:16 for WAFO by m2html © 2003