% 
% 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; 

delta   = 0.2; 
epsilon = 0.7; 
rho     = 0.4; 
phi     = 0.6; 
gamma   = 0.8; 

P = [ 1-delta, delta*(1-epsilon), delta*epsilon; 
      0, 1-phi, phi; 
      gamma*rho, (1-gamma)*rho, 1-rho ]; 

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 

% the ratios requested: 
[ x(1)/x(2), x(1)/x(3) ] 

disp('hit enter to continue...'); 
pause 

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 

%M = ImP\UmD; 
%R = M + diag(mean_recurrence_times)

% some of the same calculations in double precision: 
%
format short; 
x = A\b 

% the ratios requested: 
[ x(1)/x(2), x(1)/x(3) ] 

% solve the equations by hand: 
% 

T_21 = (5/3)*(195/22); 
T_12 = 5/3; 
T_01 = 5*(1+(7/50)*T_21);
T_02 = 5*(1+(3/50)*T_12); 

A = [ -3/50, -7/50; 
        3/5, -3/5 ; ];
      %-2/25,  2/5 ];
b = [ -97/120; 1; ]; % 1 ];
x = A\b; 
T_10 = x(1); 
T_20 = x(2); 

M = [ 0,    T_01, T_02; 
      T_10,    0, T_12; 
      T_20, T_21,    0 ]; 
M

R = diag( mean_recurrence_times ) + M; 
R