When I am working with Markov Chains I'm normally looking at stationary distributions which, by definition, don't depend on the initial state. Alternatively, if I'm looking at a set of emissions, I typically do so after the burn in.
I'm attaching some code to show a simple use case.
%%Estimating the convergence probabilities for a Markov Chain
clear all
clc
% Assume a triple modular redundant system (three independent modules)
m = 3;
% At any given point in time, a module is either working or broken
% The probability that a module that is working at time T is broken at
% time T+1 is given by
pFail = .05; % Probability of a failure
% The probability that a module that is broken at time T is repaired at
% time T+1 is given by
pFix = .40; % Probability of a fix
% The system has 8 possible states
% (0,0,0) All 3 modules are working
% (0,0,1) Modules 2 and 3 are working, module 1 is broken
% (0,1,0) Modules 1 and 3 are working, module 2 is broken
% (0,1,1) Modules 1 and 2 are broken, module 3 is working
% (1,0,0) Modules 1 and 2 are working, module 3 is broken
% (1,0,1) Modules 1 and 3 are broken, module 2 is working
% (1,1,0) Modules 2 and 3 are broken, module 1 is working
% (1,1,1) All three modules are broken
% The transition matrix for the system is
n = 2^m; % Number of states
[x, y] = meshgrid(0:n-1, 0:n-1); % All possible combos of state1 and state2
state1 = dec2bin(x)-48; % Convert to binary
state2 = dec2bin(y)-48; % Convert to binary
P = pFail.^(state1 == 0 & state2 == 1) .* ... Subsystem fails
(1-pFail).^(state1 == 0 & state2 == 0) .* ... Subsystem does not fail
pFix.^(state1 == 1 & state2 == 0) .* ... Subsystem is fixed
(1-pFix).^(state1 == 1 & state2 == 1); % Subsystem is not fixed
% P contains individual probabilities for each subsystem. Multiply across
% the row to calculate total probability assuming subsystem events are
% independent. Then reshape and transpose to get right symmetric
P = reshape(prod(P,2), n, n);
P = P'
Simple illustrative example
% Let Pi = a vector of probabilities that describe the initial state of the
% system.
% Chose a random starting position for Pi at T = 0
Pi = 100 * rand(1, length(P));
Pi = Pi/sum(Pi)
% Multiple Pi * the transition matrix to get a vector of describing the
% expected state of the system at T = 1
Pi = Pi * P
% Repeat 19 more times to get the expected state of the system at time T =
% 20
for i = 1: 19
Pi = Pi * P;
end
% After 20 iterations, Pi will have converged very close to the steady
% state of the system
Pi
% Alternatively (and easier)
Pi_other = P^20;
Pi_other = Pi_other(1,:)
%%Use Eigenvectors to calculate the steady state
% The purpose of the original example was to illustrate that the steady
% state probability vector can be described using the equation
% Pi*P = Pi We can use this information to create a system of linear
% equations and solve for Pi
[V D] = eigs(P');
Steady_State = V(:,1)/sum(V(:,1))
