markov_return.m

Markov分析的matlab工具包

%
%   随机数学(30250143, 07Fall) 大作业程序
%
%   作者： 王  音 2004011448 
%         葛元恺 2004011425 <[e-mail]>
%
%   All rights reserved. 2008.1

%   参看文献《马尔可夫基础及其应用》P151-P161

function result = Markov_Return(P)

[m n] = size(P);

if m~=n,
    error('P must be a square-matrix (m=n)');
end

N = n;
M = 0;
Z = calc_Z(P);
Z0= Z;

while(2^M < (N-1)*(N-1)),   %计算出符合要求的最小Z^(2^M)
    Z = calc_Z(Z*Z);
    M = M + 1;
end

disp(' ');              %空一行
Diff = 1:n;

while (sum(sum(Z))>0),  %判断Z是否已经为全零阵

    L=Z;
    
    C = find(sum(Z)==max(sum(Z)));
    C0 = C(1);          %找到一个常返态
    d = 1;              %记录周期信息
    
    G0 = calc_G(L,C0);  
    C  = G0;

    L = calc_Z(L*Z0);
    Gd = calc_G(L,C0);

    while Gd~=G0,   %计算各个Gd，参看P152中间
        C = union(C,Gd);
        L = calc_Z(L*Z0);
        Gd = calc_G(L,C0);
        d = d+1;
    end

    D = C;          %计算D，参看P153上侧
    for i=1:length(C),
        for j=1:n,
            if Z(j,C(i))==1,
                D = union(D,j);
            end
        end
    end

    for i=1:length(D),%把D中对于状态的行和列都置0，见P153中间
        Z(:,D(i))=0;
        Z(D(i),:)=0;
    end
    
    disp(sprintf('找到常返态闭集，周期为%d',d));    
    disp(C);
        
    Diff = setdiff(Diff,C); %剩下的状态

end

if Diff,                    %打印出剩下的非常返的状态
    disp('其余非常返闭集');
    disp(Diff);
end