Strange results with asymptotics function
I have a Markov chain with K communication classes, each communication class forms an ergodic subchain. Then I use the asymptotics function to get the K limit distributions of this Markov chain. However, results from the asymptotics function seems very sensitive to the transition matrix. Then I impliment the asymptotics function by my self following the Spedicato’s algorithm (ref), i.e.:
And I found that Spedicato’s algorithm is not that sensitive to the transition matrix.
Here is the code:
%% load data
clear;clc;
load(‘data.mat’); % load the example transition matrix A and indicator matrix C
% note: A is a left stochastic matrix, i.e., all cols sum to 1
% A2, C2 is slightly changed from A1, C1
% disp(sum(sum(abs(A1-A2)))); %7.1385e-05
if true
A=A1;
C=C1;
else
A=A2;
C=C2;
end
%% matlab asymptotics function
K = 3;
S=size(A,1);
mc = dtmc(A’);
[pi, tmix] = asymptotics(mc);
pi=pi’;
prss = sum(C*pi);
fprintf(‘prss=%sn’, mat2str(prss));
fprintf(‘prss diff = %.3f%%n’, (max(prss)-min(prss))/min(prss)*100);
%% Spedicato’s algorithm
[bins,ClassStates,ClassRecurrence,ClassPeriod] = classify(mc);
num_classes = size(ClassStates, 2);
prss = zeros(1, num_classes);
for c =1 : num_classes
sc = subchain(mc,str2num(ClassStates{c}(1)));
Cs = zeros(size(C,1), sc.NumStates);
for i=1:sc.NumStates
idx = str2num(sc.StateNames(i));
Cs(:,i) = C(:, idx);
end
As = sc.P’;
% get stationary distribution
% the condition number is not very large, so I think using inv() is ok
pi = inv(ones(sc.NumStates)+eye(sc.NumStates)-As)*ones(sc.NumStates,1);
prss(c) = sum(Cs * pi);
end
fprintf(‘prss=%sn’, mat2str(prss));
fprintf(‘prss diff = %.3f%%n’, (max(prss)-min(prss))/min(prss)*100);I have a Markov chain with K communication classes, each communication class forms an ergodic subchain. Then I use the asymptotics function to get the K limit distributions of this Markov chain. However, results from the asymptotics function seems very sensitive to the transition matrix. Then I impliment the asymptotics function by my self following the Spedicato’s algorithm (ref), i.e.:
And I found that Spedicato’s algorithm is not that sensitive to the transition matrix.
Here is the code:
%% load data
clear;clc;
load(‘data.mat’); % load the example transition matrix A and indicator matrix C
% note: A is a left stochastic matrix, i.e., all cols sum to 1
% A2, C2 is slightly changed from A1, C1
% disp(sum(sum(abs(A1-A2)))); %7.1385e-05
if true
A=A1;
C=C1;
else
A=A2;
C=C2;
end
%% matlab asymptotics function
K = 3;
S=size(A,1);
mc = dtmc(A’);
[pi, tmix] = asymptotics(mc);
pi=pi’;
prss = sum(C*pi);
fprintf(‘prss=%sn’, mat2str(prss));
fprintf(‘prss diff = %.3f%%n’, (max(prss)-min(prss))/min(prss)*100);
%% Spedicato’s algorithm
[bins,ClassStates,ClassRecurrence,ClassPeriod] = classify(mc);
num_classes = size(ClassStates, 2);
prss = zeros(1, num_classes);
for c =1 : num_classes
sc = subchain(mc,str2num(ClassStates{c}(1)));
Cs = zeros(size(C,1), sc.NumStates);
for i=1:sc.NumStates
idx = str2num(sc.StateNames(i));
Cs(:,i) = C(:, idx);
end
As = sc.P’;
% get stationary distribution
% the condition number is not very large, so I think using inv() is ok
pi = inv(ones(sc.NumStates)+eye(sc.NumStates)-As)*ones(sc.NumStates,1);
prss(c) = sum(Cs * pi);
end
fprintf(‘prss=%sn’, mat2str(prss));
fprintf(‘prss diff = %.3f%%n’, (max(prss)-min(prss))/min(prss)*100); I have a Markov chain with K communication classes, each communication class forms an ergodic subchain. Then I use the asymptotics function to get the K limit distributions of this Markov chain. However, results from the asymptotics function seems very sensitive to the transition matrix. Then I impliment the asymptotics function by my self following the Spedicato’s algorithm (ref), i.e.:
And I found that Spedicato’s algorithm is not that sensitive to the transition matrix.
Here is the code:
%% load data
clear;clc;
load(‘data.mat’); % load the example transition matrix A and indicator matrix C
% note: A is a left stochastic matrix, i.e., all cols sum to 1
% A2, C2 is slightly changed from A1, C1
% disp(sum(sum(abs(A1-A2)))); %7.1385e-05
if true
A=A1;
C=C1;
else
A=A2;
C=C2;
end
%% matlab asymptotics function
K = 3;
S=size(A,1);
mc = dtmc(A’);
[pi, tmix] = asymptotics(mc);
pi=pi’;
prss = sum(C*pi);
fprintf(‘prss=%sn’, mat2str(prss));
fprintf(‘prss diff = %.3f%%n’, (max(prss)-min(prss))/min(prss)*100);
%% Spedicato’s algorithm
[bins,ClassStates,ClassRecurrence,ClassPeriod] = classify(mc);
num_classes = size(ClassStates, 2);
prss = zeros(1, num_classes);
for c =1 : num_classes
sc = subchain(mc,str2num(ClassStates{c}(1)));
Cs = zeros(size(C,1), sc.NumStates);
for i=1:sc.NumStates
idx = str2num(sc.StateNames(i));
Cs(:,i) = C(:, idx);
end
As = sc.P’;
% get stationary distribution
% the condition number is not very large, so I think using inv() is ok
pi = inv(ones(sc.NumStates)+eye(sc.NumStates)-As)*ones(sc.NumStates,1);
prss(c) = sum(Cs * pi);
end
fprintf(‘prss=%sn’, mat2str(prss));
fprintf(‘prss diff = %.3f%%n’, (max(prss)-min(prss))/min(prss)*100); markov chain, asymptotics MATLAB Answers — New Questions