m=[1 0.1; 0.5 1]
day=[1 0]
for i=1:10000
    temp=day(i,:)*(m.*rand(2));
    day=[day; temp>min(temp)];
end
sum(day)

%%
D=[1,3,3,2,1,2,1,1,2,2,3,1,3,1,2,3,1,2,3,2,3,1,3,3,2,3,3,2,3,1,3,2,3,3,2,3,2,3,2,3,2,3,2,1,3,2,3,2,2,3,2,1,3,1,2,3,2,2,2,1,3,1,3,1,2,2,2,3,1,1,2,1,2,3,1,2,2,3,2,3,1,3,1,3,3,2,2,3,1,2,2,3,1,2,3,1,2,3,1,3,3,2,3,2,2,3,1,2,1,1,2,3,2,3,2,3,1,3,1,3,2,3,3,2,1,3,2,3,2,3,2,2,2,2,3,2,2,3,3,2,1,3,2,1,3,1,3,1,3,3,3,2,3,2,1,3,2,2,2,3,2,3,1,3,1,2,3,2,3,2,2,2,3,2,2,1,3,2,2,1,2,3,2,1,3,2,3,3,2,1,3,2,2,2,1,3,1,1,3,2,1,3,3,1,3,2,3,3,3,3,2,3,2,2,2,2,2,1,3,2,2,2,2,2,2,2,1,3,2,1,1,2,1,3,3,1,2,3,1,3,3,1,1,2,3,3,3,2,3,2,1,1,1,2,1,2,3,1,1,3,1,3,1,3,2,2,1,2,2,2,1,1,3,2,2,2,2,2,3,1,3,3,2,1,3,2,3,2,2,3,2,3,3,3,1,2,2,2,3,2,3,2,2,2,3,3,2,3,1,2,2,2,2,3,3,1,2,2,2,3,2,1,3,3,3,2,3,3,1,3,1,1,2,2,3,2,3,2,3,3,1,3,1,3,2,2,3,1,3,3,2,2,1,3,3,2,2,1,3,2,2,1,3,3,2,3,2,2,3,2,3,2,2,3,1,3,2,3,2,3,2,2,1,2,3,1,3,2,3,2,2,1,2,2,1,3,2,3,2,2,1,2,1,3,2,3,2,3,3,2,2,3,2,1,3,2,1,2,2,3,3,3,2,3,3,1,2,3,2,2,3,2,1,3,2,3,2,2,1,2,1,2,3,2,2,2,2,3,2,3,1,3,1,3,2,2,3,1,1,2,2,3,2,3,3,2,3,1,3,2,2,3,3,1,2,2,3,1,1,3,2,3,1,3,1,3,2,3,1,3,2,2,3,1,2,3,2,2,2,3,3,1,3,2,1,1,2,2,2,2,2,1,1,2,1,1,3,1,2,3,1,2,3,2,3,2,2,3,2,3,1,3,2,3,1,2,2,3,2,2,1,2,1,3,3,3,3,2,3,2,3,1,2,2,1,3,3,1,2,3,2,1,2,1,3,2,2,1,3,2,3,2,1,2,2,2,2,2,3,2,3,1,3,1,3,2,3,3,1,3,2,2,2,3,1,3,2,2,2,1,3,3,2,1,3,2,1,1,2,2,1,3,2,1,2,3,2,3,2,3,3,1,3,3,3,2,1,1,3,1,1,3,2,2,3,1,3,2,3,2,1,3,3,3,1,2,3,1,1,1,3,2,1,3,2,1,2,3,1,3,2,3,2,1,3,2,2,2,2,3,2,2,2,3,1,2,1,1,2,1,3,3,1,3,3,2,2,2,1,3,2,2,1,3,1,3,1,3,2,3,3,3,1,3,2,2,3,2,2,3,3,1,3,1,3,2,3,2,1,3,3,3,2,3,2,2,3,2,2,2,2,2,2,3,1,1,3,2,3,3,2,3,1,3,2,1,3,1,3,2,3,1,2,3,2,1,3,2,2,3,3,2,2,2,2,2,2,3,3,2,2,2,3,3,1,3,1,3,3,2,2,2,3,1,2,2,3,3,2,3,3,1,2,3,2,2,2,1,2,3,3,2,1,2,3,1,3,2,2,3,2,3,2,2,3,2,2,3,3,3,3,1,3,3,2,3,1,3,2,3,2,2,2,3,1,3,1,2,2,2,3,2,3,2,3,2,2,3,1,1,3,1,3,2,3,2,3,3,3,2,3,1,3,3,1,3,1,3,1,3,2,3,1,3,1,3,3,1,3,2,2,3,2,2,1,2,1,1,2,1,3,1,2,3,2,3,1,2,3,2,1,1,3,2,1,3,2,1,2,2,2,2,2,3,2,1,3,2,1,2,2,3,1,2,3,3,2,2,1,3,1,3,2,3,1,3,1,3,2,1,2,2,1,3,1,3,2,2,2,3,2,2,3,2,3,1,3,3,2,2,3,2,2,3,2,1,2,2,3,3,1,3,2,1,2,2,3,3,2,1,3,1,3,1,3,2,2,2,3,2,1,3,1,3,1,2,2,2,3,1,3];
S=(D==1);
R=(D==2);
C=(D==3);
%Total S
TS=sum(S)
%S after S
SS=sum(sum([S(1:end-1) ; S(2:end)])==2)
%R after S
SR=sum(sum([S(1:end-1) ; R(2:end)])==2)
%C after S
SC=sum(sum([S(1:end-1) ; C(2:end)])==2)
%percentage
SS/TS
SR/TS
SC/TS

%%
D=[1,3,3,2,1,2,1,1,2,2,3,1,3,1,2,3,1,2,3,2,3,1,3,3,2,3,3,2,3,1,3,2,3,3,2,3,2,3,2,3,2,3,2,1,3,2,3,2,2,3,2,1,3,1,2,3,2,2,2,1,3,1,3,1,2,2,2,3,1,1,2,1,2,3,1,2,2,3,2,3,1,3,1,3,3,2,2,3,1,2,2,3,1,2,3,1,2,3,1,3,3,2,3,2,2,3,1,2,1,1,2,3,2,3,2,3,1,3,1,3,2,3,3,2,1,3,2,3,2,3,2,2,2,2,3,2,2,3,3,2,1,3,2,1,3,1,3,1,3,3,3,2,3,2,1,3,2,2,2,3,2,3,1,3,1,2,3,2,3,2,2,2,3,2,2,1,3,2,2,1,2,3,2,1,3,2,3,3,2,1,3,2,2,2,1,3,1,1,3,2,1,3,3,1,3,2,3,3,3,3,2,3,2,2,2,2,2,1,3,2,2,2,2,2,2,2,1,3,2,1,1,2,1,3,3,1,2,3,1,3,3,1,1,2,3,3,3,2,3,2,1,1,1,2,1,2,3,1,1,3,1,3,1,3,2,2,1,2,2,2,1,1,3,2,2,2,2,2,3,1,3,3,2,1,3,2,3,2,2,3,2,3,3,3,1,2,2,2,3,2,3,2,2,2,3,3,2,3,1,2,2,2,2,3,3,1,2,2,2,3,2,1,3,3,3,2,3,3,1,3,1,1,2,2,3,2,3,2,3,3,1,3,1,3,2,2,3,1,3,3,2,2,1,3,3,2,2,1,3,2,2,1,3,3,2,3,2,2,3,2,3,2,2,3,1,3,2,3,2,3,2,2,1,2,3,1,3,2,3,2,2,1,2,2,1,3,2,3,2,2,1,2,1,3,2,3,2,3,3,2,2,3,2,1,3,2,1,2,2,3,3,3,2,3,3,1,2,3,2,2,3,2,1,3,2,3,2,2,1,2,1,2,3,2,2,2,2,3,2,3,1,3,1,3,2,2,3,1,1,2,2,3,2,3,3,2,3,1,3,2,2,3,3,1,2,2,3,1,1,3,2,3,1,3,1,3,2,3,1,3,2,2,3,1,2,3,2,2,2,3,3,1,3,2,1,1,2,2,2,2,2,1,1,2,1,1,3,1,2,3,1,2,3,2,3,2,2,3,2,3,1,3,2,3,1,2,2,3,2,2,1,2,1,3,3,3,3,2,3,2,3,1,2,2,1,3,3,1,2,3,2,1,2,1,3,2,2,1,3,2,3,2,1,2,2,2,2,2,3,2,3,1,3,1,3,2,3,3,1,3,2,2,2,3,1,3,2,2,2,1,3,3,2,1,3,2,1,1,2,2,1,3,2,1,2,3,2,3,2,3,3,1,3,3,3,2,1,1,3,1,1,3,2,2,3,1,3,2,3,2,1,3,3,3,1,2,3,1,1,1,3,2,1,3,2,1,2,3,1,3,2,3,2,1,3,2,2,2,2,3,2,2,2,3,1,2,1,1,2,1,3,3,1,3,3,2,2,2,1,3,2,2,1,3,1,3,1,3,2,3,3,3,1,3,2,2,3,2,2,3,3,1,3,1,3,2,3,2,1,3,3,3,2,3,2,2,3,2,2,2,2,2,2,3,1,1,3,2,3,3,2,3,1,3,2,1,3,1,3,2,3,1,2,3,2,1,3,2,2,3,3,2,2,2,2,2,2,3,3,2,2,2,3,3,1,3,1,3,3,2,2,2,3,1,2,2,3,3,2,3,3,1,2,3,2,2,2,1,2,3,3,2,1,2,3,1,3,2,2,3,2,3,2,2,3,2,2,3,3,3,3,1,3,3,2,3,1,3,2,3,2,2,2,3,1,3,1,2,2,2,3,2,3,2,3,2,2,3,1,1,3,1,3,2,3,2,3,3,3,2,3,1,3,3,1,3,1,3,1,3,2,3,1,3,1,3,3,1,3,2,2,3,2,2,1,2,1,1,2,1,3,1,2,3,2,3,1,2,3,2,1,1,3,2,1,3,2,1,2,2,2,2,2,3,2,1,3,2,1,2,2,3,1,2,3,3,2,2,1,3,1,3,2,3,1,3,1,3,2,1,2,2,1,3,1,3,2,2,2,3,2,2,3,2,3,1,3,3,2,2,3,2,2,3,2,1,2,2,3,3,1,3,2,1,2,2,3,3,2,1,3,1,3,1,3,2,2,2,3,2,1,3,1,3,1,2,2,2,3,1,3];
for i=1:3
    for j=1:3
        model(i,j)=sum(sum([0 D==i ; D==j 0])==2) / sum(D==i);
    end
end

model'