% 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,:);

% statics-matlab.ppt , Page 3
% generate discrete r.v.
sampleN = 100;
P = [0.05 0.1 0.15 0.18 0.22 0.3];
Xvalue = [1 2 3 4 5 6];
n =length(P);
SumP = zeros(n, 1);
for i=1:n,
    SumP(i) = sum(P(1:i));
end
U = rand(sampleN,1);
Sample = zeros(sampleN,1);
for i=1:sampleN,
    for j=1:n,
        if U(i)< SumP(j),
            Sample(i) = Xvalue(j);
            break;
        end
    end
end

% draw CDF, Page 5
F = zeros(n,1);
for i = 1:n,
    for j=1:sampleN,
        if Sample(j)<= Xvalue(i),
            F(i) = F(i) + 1;
        end
    end
end
F = F./sampleN;
plot(Xvalue, F, '-ob');
xlabel('Random Variable X'); ylabel('Cumulative distr. function (CDF)');
hold on;
plot(SumP, '--^r');
legend('Simulation', 'Analytical');
title('CDF');

% draw pmf, Page 6
PMF = zeros(n,1);
for i=1:n,
    for j=1:sampleN,
        if Sample(j) == Xvalue(i),
            PMF(i)= PMF(i)+1;
        end
    end
end
PMF = PMF./sampleN;
figure;
plot(PMF, '-ob');
xlabel('Random Variable X'); ylabel('Cumulative distr. function (CDF)');
hold on;
plot(P, '--^r');
legend('Simulation', 'Analytical');
title('Probability Mass Function');

% The On-Off discrete-time Markov chain example (page 11)
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
State=zeros(stateNum, 1); 
State(init_state) = 1;  % count the initial state
currentState = init_state;
for i=1:sampleN,
    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');
legend('theory', 'simulation');

% 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
Edge = 0:0.05:5;
N = histc(Sample, Edge);
bar(Edge, N, 'histc');


