% set up parameters

Pw = .080527+0.01020;
Ps = .1473;
Pd = .0475;
Pt = .00501;
Ph = .02945; 
Po = 1-Pw-Ps-Pd-Pt-Ph;
count = 20;

% set up probability matrices

A = [Ph,Ph,Ph,Ph,Ph,Ph,Ph,Ph;
     Ps+Pw,0,Ps,Ps,0,0,Ps,0;
     Pd,0,Pd,Pd,0,0,Pd,0;
     Pt,Pt,Pt,Pt,Pt,Pt,Pt,Pt;
     0,Ps+Pw,Pw,0,Ps,Ps,0,Ps;
     0,0,0,Pw,0,0,0,0;
     0,Pd,0,0,Pd,Pd,0,Pd;
     0,0,0,0,Pw,Pw,Pw,Pw];
B = Po*eye(8);
Z = zeros(8,8);
Z1 = zeros(8,1);
Z2 = Z1';
O3 = Po*ones(1,8);

P = [A,Z,Z,Z1;
     B,A,Z,Z1;
     Z,B,A,Z1;
     Z2,Z2,O3,1];
  
%check  
check = ones(1,25);
check*P;

% eigenvalues as a curiosity
figure(1);
eigplot(P);

% set up initial probaliites
U0 = zeros(25,1);
U0(1,1) =1;

% show iterated distributions
U = U0;
V = U0;
for j=1:count
   U=P*U;
   V = [U,V];
end;
figure(2);
surf(V);
axis tight

% show iterated transition matrices
figure(3)
for j=1:count
   mesh(P^j);
   axis tight
   %pause(1)
   pause
end;


% set up run matrix
RA = [1,2,2,2,3,3,3,4;
      0,1,1,1,2,2,2,3; 
      0,1,1,1,2,2,2,3;
      0,1,1,1,2,2,2,3;
      0,0,0,0,1,1,1,2;
      0,0,0,0,1,1,1,2;      
      0,0,0,0,1,1,1,2;         
      0,0,0,0,0,0,0,1];   
R = [RA,Z, Z, Z1;
     Z, RA,Z, Z1;
     Z, Z, RA,Z1;
     Z2,Z2,Z2,0 ];
  
% compute the number of runs  
tnr = 0;
trh = [];
cnr = 0;
crh = [];
U = U0;
for j=1:count
   cnr =sum(sum(R.*(P*diag(U))));
   tnr = tnr+cnr;
   trh =[trh,tnr];
   crh= [crh,cnr];
   U = P*U;
end;

figure(4)
plot(1:count,trh,'o')
axis([1,count,0,1])
figure(5)
plot(1:count,crh,'o')
axis([1,count,0,1])