

% written by StudentDave
%for licensing and usage questions
%email scienceguy5000 at gmail. com

%% Welcome to the BAYESIAN DOJO! Here, we'll learn about markov chains

% our main examples will be of Ergodic regular markov chains
% these type of chains converge to a steady-state, and have some nice
% properties for rapid calculation of this steady state.


%% non-regular
%first, an example of a non-regular, but ergodic markov chain, which
%doesn't converge

P = [ 0 1 ;1 0]

t_all = [];
i_all = [];
figure(1)
clf
for i = 1:100
    t = P^i
    t_all = [t_all t(:)];
    i_all = [i_all  ones(4,1)*i];
    subplot(211)
    draw_states(t,i)
    subplot(212)
    plot(i_all',t_all','.-')
    xlabel('discrete time steps')
    ylabel('probability')
    title('evolution of transition probs. for each element')
    pause
end



%% Has zero elements in initial transition matrix (TM)), but not always, so
% thus still regular
P = [ 1/2 1/2 ;1 0]


t_all = [];
i_all = [];
figure(1)
clf
for i = 1:100
    t = P^i
    t_all = [t_all t(:)];
    i_all = [i_all  ones(4,1)*i];
    subplot(211)
    draw_states(t,i)
    subplot(212)
    plot(i_all',t_all','.-')
    xlabel('discrete time steps')
    ylabel('probability')
    title('evolution of transition probs. for each element')
    pause
end

%% Training: level one! Learn to read your opponent
% Our Bayesian Ninja has just learned about the sudden reappearance of an
% old enemy of the Bayesian Clan, The frequentisian Ninja Clan!
% The bayesian ninja must now train to fight!
% One critical skill is to know your opponent, know their tendencies, and
% learn their patterns!  The old master ninja is the last of the Bayesians who have fought the
% Frequentisian Ninja, and he describes the different close-combat fighting styles
% taught within their devious clan in terms of a markov process.

% In the first lesson, he describes a three state fighting style, comprised of
% 1)punch (red) and 2)kick (yellow), 3) and flying falcon punch (blue).
% here we look at how, overall, the Frequentisian Ninja's will fight, given
% the probabilities of how they mix up their punching, kicking, and
% "special attack"
%

% E honda style (Likes to punch)-------------------------------

a= .9
b = .3
c = .2
P = [ a (1-a)/2 (1-a)/2;  (1-b)/2 b (1-b)/2; (1-c)/2 (1-c)/2 c ;]


% chun Li Style (Likes to kick)-------------------------------

a= .1
b = .9
c = .3
P = [ a (1-a)/2 (1-a)/2;  (1-b)/2 b (1-b)/2; (1-c)/2 (1-c)/2 c ;]

% Captain Falcon Style (Obvious? :) -------------------------------

a= .1
b = .2
c = .9
P = [ a (1-a)/2 (1-a)/2;  (1-b)/2 b (1-b)/2; (1-c)/2 (1-c)/2 c ;]

% MASTER Frequentisian Ninja (Perfectly equal skill in each

a = .33333
b = .333333
c = .333333
P = [ a (1-a)/2 (1-a)/2;  (1-b)/2 b (1-b)/2; (1-c)/2 (1-c)/2 c ;]

t_all = [];
i_all = [];
figure(1)
clf
for i = 1:100
    t = P^i
    t_all = [t_all t(:)];
    i_all = [i_all  ones(size(t_all,1),1)*i];
    subplot(211)
    draw_states3(t(:),i)
    subplot(212)
    plot(i_all',t_all','.-')
    xlabel('discrete time steps')
    ylabel('probability')
    title('evolution of transition probs. for each element')
    axis([0 max(max(i_all)) min(min(t_all))-.5 max(max(t_all))+.5])
    pause
end



%% Training: level two! What's it take to beat the MASTER Frequentisian Ninja?
% here, we show how, even with just a 1st order markov chain, you can still
% simulation systems that depend on past events.
% If the Master ninja can land a 3 hit combo in the specific order of
% punch, kick, falcon punch...the bayesian ninja will get KO'd!
% Simulate this as a markov process and see how the Bayesian will perform
% given his ability, b, to interrupt the punches, and starting state u,
% after T number of attacks (time steps)


a= .5
b = .7 %interrupt probability
P = [1-a a 0 0; b 0 1-b 0; b 0 0 1-b; 1 0 0 0]

t_all = [];
i_all = [];
figure(1)
clf

for i = 1:100
    t = P^i
    t_all = [t_all t(:)];
    i_all = [i_all  ones(size(t_all,1),1)*i];
    subplot(211)
    draw_states4(t,i)
    subplot(212)
    plot(i_all',t_all','.-')
    xlabel(['Time steps = ', num2str(i)])
    ylabel('probability')
    title('evolution of transition probs. for each element')
    pause
end


%how to solve with eigen math :)--------------------

% get eigen decomposition 
[Evector,Evalue] = eig(P')
%get values out of matrix
values     = diag(Evalue);
%find the unitary Evalue, won't be exact, so use generic tool for finding
%closest element to N. There will only be one for regular transition matrix

N = 1
[min_v,coln] = min(abs(values-N))
%grab the corresponding vector, and normalize, this is your stationary
%distribution!
Evector = Evector(:,coln)

fixed_row_vector = (Evector/sum(Evector))'


