I am trying to solve a set of equations to determine the stationary distribution of an ergodic Markov matrix.
Namely, the matrix is
P=[0 0 0 0.5 0 0.5;
0.1 0.1 0 0.4 0 0.4;
0 0.2 0.2 0.3 0 0.3;
0 0 0.3 0.5 0 0.2;
0 0 0 0.4 0.6 0;
0 0 0 0 0.4 0.6];
and the set of equations are the ones from the theorem below
How can I convert the equations above into valid Matlab syntax?
The stationary distribution is given by the left eigen vector with eigen-value 1.
>> [V D] = eig( P.' ); %// note the transpose .' - we are looking for the **left** EV
>> st = V(:,1).'; %//' the stationary distribution
st =
0.0051 0.0509 0.2291 0.6110 0.5346 0.5346
>> D(1)
1.0000
This is a different approach from @Shai's answer.
An alternate way to do it comes from solving the Pi*P = Pi equations for the for steady state and ignore the requirement that the sum of the pi_j needs to be one (for now).
A little matrix algebra...
Then we know Pi does not have a unique solution to this without the "sum to 1" requirement. Pi must be in the nullspace of (transpose(P) - I). MATLAB is good at this. Normalization after gives the desired result (as pointed out by @LuisMendo in the comments).
P=[0 0 0 0.5 0 0.5;
0.1 0.1 0 0.4 0 0.4;
0 0.2 0.2 0.3 0 0.3;
0 0 0.3 0.5 0 0.2;
0 0 0 0.4 0.6 0;
0 0 0 0 0.4 0.6];
I = eye(size(P));
Y = null(P'-I)
PI = Y./(sum(Y))
This is easy to check.
>> PI(:)' % Force into row vector
ans =
0.0026 0.0259 0.1166 0.3109 0.2720 0.2720
Compare with the 25-step transition matrix.
P5 = P*P*P*P*P;
P25 = P5*P5*P5*P5*P5;
>>P25 =
0.0026 0.0259 0.1166 0.3109 0.2720 0.2720
0.0026 0.0259 0.1166 0.3109 0.2720 0.2720
0.0026 0.0259 0.1166 0.3109 0.2720 0.2720
0.0026 0.0259 0.1166 0.3109 0.2720 0.2720
0.0026 0.0259 0.1166 0.3109 0.2720 0.2720
0.0026 0.0259 0.1166 0.3109 0.2720 0.2720
