% The On-Off discrete-time Markov chain example 
clear all;
alpha= 0.6; beta = 0.9;
Pi_theory = [beta/(alpha+beta); alpha/(alpha+beta)];
sampleN= 100;  stateNum = 2;  % 2 states
U = rand(sampleN,1);
init_state = randi(2);  %either 1 or 2
StateVector = zeros(1, sampleN); % record system's state according to time
State=zeros(stateNum,1); 
State(init_state) = 1;  % count the initial state
currentState = init_state;
for i=1:sampleN,
    StateVector(i) = currentState;
    if currentState ==1, 
        if U(i) > alpha,
            nextState = 1;
        else
            nextState = 2;
        end
    else  % currentState = 2
        if U(i) > beta,
            nextState = 2;
        else
            nextState = 1;
        end
    end
    State(nextState) =State(nextState)+1; % count the current state
    currentState = nextState;
end
Pi_simulation = State./sampleN;
figure;
%plot(Pi_theory, 'r-^'); hold on; plot(Pi_simulation, 'k:o');
bar([Pi_theory Pi_simulation]); 
xlabel('system state'); ylabel('probability');
legend('theory', 'simulation');
% draw the simulation run scenario according to time
figure;
plot(1:20, StateVector(1:20), 'k-o'); xlabel('discrete time'); ylabel('system state');
axis([1 20 0 4]);
return;

% draw pdf for exponential distribution.
lambda = 2; sampleN = 1000;
Sample = zeros(1, sampleN);
U = rand(1, sampleN);
for i=1:sampleN,
    Sample(i) = -log(1-U(i))/lambda;
end
figure; Edge = 0:0.5:5;
N = histc(Sample, Edge);
bar(Edge, N, 'histc');

% Matlab.ppt,  compute discrete-time MC steady state
P = [0.4 0.6 0 0; 0.12 0.46 0.42 0; 0.036 0.222 0.448 0.294; 0.027 0.189 0.441 0.343];
I = eye(4);
Pi = zeros(1, 4);   Unit = ones(4, 1);
A = P - I;
A = [Unit A(:, 2:4)];
Pi = [ 1 0 0 0]*inv(A);
%  second way to calculate pi.
PP = P;
for i=1:100,
    PP = P*PP;
end
Pi = PP(1,:);


