% PRACTICA de CADENAS de MARKOV con MatLab.

% Sea la matriz de transicion de una cadena de Markov
P = [.1 .7 .2
    .4 0 .6
    0 .5 .5];

% La matriz de transicion despues de 20 etapas es:
P^20

% Supongamos la distribucion inicial:
p0 = [.3 .4 .3];

% La distribucion de los estados despues de 20 etapas es
p20 = p0*(P^20)

%.....................................................................

% Calculas autovalores y autovectores:
[c,d] = eig(P);

% c es la matriz de autovectores colocados en columna
c

% d es una matriz diagonal con los autovalores
d

% Reescalas los autovectores, para que uno de los componentes sea
% igual a 1
for (i=1:length(P))
	c(:,i) = c(:,i)./min(c(:,i));
end

%.....................................................................

% Calculas la distribucion limite

% Forma chapuza de calcular la distribucion limite:
Plim = c*(d^(1.E10))*inv(c)

% Si esta instalada la toolbox Simbolic:
n = sym('n')
Q = d;
for (i=1:length(P))
pre = limit(d(i,i)^n,n,inf);
Q(i,i) = double(pre);
end
Plim = c*Q*inv(c)

%.....................................................................

% Calculas la distribucion estacionaria:
% pi=pi*P

% Primero guardas la matriz original
P1 = P;

% Tal como sale, pi=pi*P es un sistema linealmente dependiente
% (sobra una ecuacion): Sustituyes la ultima columna de P por una
% columna de 1's para tener en cuenta que pi(1)+ ... + pi(n) = 1
P1(:,3) = [1 1 1]';

% Construyes la matriz diagonal con un 0 en el ultimo elemento
J = diag([1 1 0], 0);

% pi=pi*P y pi*1=1 ->
% pi*(P1-J)=[0 0 1]

pi = [0 0 1]*inv(P1-J)

% Alternativa:
%
% Como pi*P=pi implica que pi es un autovector por la izquierda de la
% matriz Q asociado al autovalor 1:
% pi=pi*P -> pi-pi*P=0 -> pi*(I-P)=0

% Si pones
% [v,d] = eig(P)
% d(i,i) es el autovalor asociado al autovector v(:,i)
% Pero un autovector por la izquierda de P equivale al
% autovector por la derecha de P' para el mismo autovalor.

[v,d] = eig(P');

% Para que las probabilidades sumen uno reescalo:

pi = v(:,3)/sum(v(:,3))

% Ejemplo: probad con

Q = [0 2/3 1/3 ; 3/8 1/8 1/2 ; 1/2 1/2 0];

%.....................................................................

% Proceso de POISSON

% Se simulan las llegadas de un proceso de Poisson

% Tasa de llegadas
lambda=10;

% Tiempo total
tmax=3;

T = [];

T(1) = exprnd(1/lambda);
i=1;

while T(i) < tmax,
  pre = T(i)+ exprnd(1/lambda);

  	if(pre>tmax)
  		break
  	else
  		T(i+1)=pre;
  	end
  i=i+1;
end

Nllegadas = length(T);

% Dibujas un grafico en escalon que indica las llegadas
% frente a los tiempos en que aparecen
stairs(T(1:Nllegadas), 0:(Nllegadas-1));
xlabel('tiempos de llegadas');
ylabel('Numero dellegadas');

pause
clf

% Dibujas un grafico de los tiempos de llegadas
Y = zeros(1,Nllegadas);
plot(T(1:Nllegadas),Y,'.')
xlabel('tiempos de llegadas');