clear all
clc

trans = [0.1 0.3 0.4 0.2;
        0.2 0.1 0.3 0.4;
        0.4 0.2 0.1 0.3;
        0.3 0.4 0.2 0.1];
 
 n=1000;
 state = zeros(4,n);
 path=zeros(n);
 
 for i=1:4;
     rand1 = rand/2;
     rand2 = rand/2;
     state(:,1)=[rand1 0.5-rand1 rand2 0.5-rand2];
     for j=1:n;
        x = rand;
        y = find(x<cumsum(state(:,j)));
        path(j,i) = min(y);
        state(:,j+1)=trans(path(j,i),:)';

         for k=1:4;
             count(j,i,k)=length(find(path(1:j,i)==k))/j;
         end;
     end;
 end;

%  for i=1:4;
%      figure(i);
%      stairs(path(:,i));
%      xlabel('time')
%      ylabel('X(i) state')
%      title(['Sample Path #' int2str(i)]);
%      saveas(figure(i),['Prob 2 -Sample path ' int2str(i) '.png'])
%  end;
% 
%  x1 = [1:1:100];
 for i=1:4;
     figure(i+4);
     hold on
     plot(x1,count(:,i,1));
     plot(x1,count(:,i,2));
     plot(x1,count(:,i,3));
     plot(x1,count(:,i,4));
     xlabel('i');
     ylabel('Relative frequency');
     title(['sample path relative freq' int2str(i)]);
     saveas(figure(i+4),['Prob 2 - sample path relative freq ' int2str(i) '.png'])
     hold off
 end;
 
 for i=1:4;
 [v,d] = eig(trans');
 v1_eigen = v/sum(v)
 v2_Rel_Freq = [count(n,i,1), count(n,i,2), count(n,i,3), count(n,i,4)]
 end;
 

f=[2 1 2.5 -1];
total2 = dot(f,v1_eigen)

n1=10000;

x(1) = ceil(4*rand);

for i=1:n1
[m,n]=sort(trans(x(i),:));
exitloop=0;
j=1;
u1=rand;
while exitloop==0;
    if u1<sum(m(1:j));
    x(i+1)=n(j);
    exitloop=1;
    else
    j=j+1;
    end
end
end


for i=1:n1;
    total1(i)=sum(f(x(i)));
end;

% for i=1:4
%     total2(i)=sum(f(i)*v1_eigen(i));
% end

avg1=sum(total1)/n1






 
 