%Aca demostramos como funciona el proceso de Markov de las fichas
%rojas y blancas, simulando varias jugadas regidas por la
%probabilidad de transición

%la matriz de transicion es
M=[ 1/2 1/2 0
    3/8 1/2 1/8
    0   1   0];

%defino la suma acumulada de probabilidades para asignar numeros
%aleatorios que cumplan con las probabilidades. Esto es necesario
%solo en la simulacion, no en la solucion matematica del problema
%de Markov
MAc=cumsum(M,2);

%elijo una condicion inicial cualquiera para el sistema
%nc veces voy a jugar (ie repito el experimento)
nc=1e6;
%voy a hacer ns pasos en la cadena de Markov
ns=10;

%en blancas voy a guardar cuantas blancas hay en la primer urna
%para cada realización/jugada y en cada
%paso del proceso.

%dimensiono el vector
blancas=zeros(ns,nc);

%elijo como condicion inicial ninguna ficha blanca.
blancas(1,:)=0;
%o puedo elegir 1 
%blancas(1,:)=1;
%o puedo elegir 2
%blancas(1,:)=2;


%probas iniciales, 

p0(1)=length(find(blancas(1,:)==0))/nc;
p1(1)=length(find(blancas(1,:)==1))/nc;
p2(1)=length(find(blancas(1,:)==2))/nc;

%evolucion de los pasos de Markov, 
for i=2:ns;

%voy a hacer la transicion con probabilidad dada por M
%tiro nc numeros al azar, una para cada jugada/realizacion

rands=rand(nc,1);

% me fijo a que transicion corresponde, aca uso la proba
% acumulada que defini arriba. 
ip1=find(rands < MAc(blancas(i-1,:)+1,1));
rands(ip1)=NaN;%'borro'
ip2=find(rands < MAc(blancas(i-1,:)+1,2));
rands(ip2)=NaN; % 'borro'
ip3=find(rands < MAc(blancas(i-1,:)+1,3));

blancas(i,ip1)=0;
blancas(i,ip2)=1;
blancas(i,ip3)=2;

%ahora calculo las probas, a cada paso, contando 
%en cuantas jugadas tengo 0, 1 o 2 fichas blancas.

p0(i)=length(find(blancas(i,:)==0))/nc;
p1(i)=length(find(blancas(i,:)==1))/nc;
p2(i)=length(find(blancas(i,:)==2))/nc;

end %termino mi ciclo de n pasos.


%ahora empiezo a analizar un poco los resultados

figure(1) % mirar en pantalla completa
plot(0:ns-1,[p0;p1;p2]','o-');legend('0','1','2')
xlabel('paso n de la cadena');
ylabel('probabilidad');
title({['Problemas de 2 fichas blancas y 4 rojas - ' num2str(nc) ' jugadas'], ...
['p_0=' num2str(p0(end)) ' p_1=' num2str(p1(end)) 'p_2=' ... 
num2str(p2(end))]});

figure(2)
%otro plot para animar
for i=1:floor(log(nc)/log(2));
cou=histc(blancas(:,1:2^i)',0:2);bar(1:10,cou'/2^i);
xlabel('paso n de la cadena');
ylabel('histograma/probabilidad');
title([num2str(2^i) ' jugadas']);
legend('0 fichas blancas','1 ficha blanca','2 fichas blancas');
drawnow;
end


figure(3)

isel=[3:5 30:10:50 300:100:500];

for i=1:3
for j=1:3
subplot('position',[0.06+0.32*(i-1),0.08+0.32*(j-1),0.25,0.25]);
ij = j + 3*(i-1);
cou=histc(blancas(:,1:isel(ij))',0:2);bar(1:10,cou'/isel(ij));
xlim([0 11]);
title([num2str(isel(ij)) ' jugadas']);
legend('0 fichas blancas','1 ficha blanca','2 fichas blancas');
end
end