% 
% Written by:
% -- 
% John L. Weatherwax                2007-09-10
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

close all; drawnow; 
clc; clear; 

format rat; 

P = [   0, 1/2,  1/2; 
      1/3,   0,  2/3; 
        0,   1,    0; ]; 

nStates = size(P,1); 

fprintf('P equals \n'); P 

% compute the matrix and the right hand side to solve: 
A = [ P.' - eye(size(P)); 1, 1, 1, ]; 
b = [ 0; 0; 0; 1 ]; 

fprintf('the steady-state solution for a, b, and c is given by \n'); 
x = A\b 

mean_recurrence_times = 1./x; 

% Compute the matrix equations we need to solve to get the "R" matrix:
ImP = eye(size(P)) - P; 
UmD = repmat( ones(1,nStates), [nStates,1] ) - diag(mean_recurrence_times);
fprintf( 'The matrix I - P is \n' ); 
ImP
fprintf( 'The matrix U - D is \n' ); 
UmD 

%ImP\UmD

T_32 = 1; 
T_23 = 8/5; 
T_12 = -3*( -7/6 + (2/3)*T_23 ); 
T_13 = -3*( 1 - T_23 ); 

A = [ -1/2, -1/2 ; 1, -2/3 ]; 
b = [ -11/2; 1 ]; 
x = A\b; 
T_21 = x(1); T_31 = x(2); 

M = [ 0,    T_12, T_13; 
      T_21,    0, T_23; 
      T_31, T_32,    0 ]; 
M

R = diag( mean_recurrence_times ) + M; 
R