%              Stationaer.M
%
%  Stationäre Wahrscheinlichkeitsverteilung einer Markovkette
%  
%
% J. Lunze
%
% für 2. Auflage: 16.6.2012
% für 3. Auflage: 11.1.2017
% für 4. Auflage: 1.8.2025
%
%  G = Matrix der Zustandsübergangswahrscheinlichkeiten der Markovkette
%
clear
close all
G=[0 0.5 1; 0.3 0.5 0; 0.7 0 0];
Markovprozess=ss(G, zeros(3,1), eye(3,3), zeros(3,1));
p0=[1; 0; 0];                   % Anfangswahrscheinlichkeitsverteilung
Tmax=30;
P1=dinitial(G, zeros(3,1), eye(3,3), zeros(3,1), p0, Tmax);
figure(1)
subplot(3,1,[1,2])
plot([0:Tmax-1], P1, 'b');
hold on
plot(0, P1(1,1), 'bo', 0, P1(1,2), 'bo', 0, P1(1,3), 'bo');
plot([0 0], [1.1, -0.1], 'k:');
latexxlabel('$$k$$');
xticks([0 5 10 15 20 25 30]);
latexylabel('$$p_i$$');
rotateY;
yticks([0 1]);
axis([-1 Tmax-1 -0.1 1.1]);
%
%  Berechnung der stationären Wahrscheinlichkeiten
[V, Lambda]=eig(G);
pbar=V(:,2);   % weil 2. Eigenwert den Wert 1 hat
pbar=pbar/sum(pbar)
% plot([30, Tmax], [pbar(1,1),pbar(1,1)], '--');
% plot([30, Tmax], [pbar(2,1),pbar(2,1)], '--');
% plot([30, Tmax], [pbar(3,1),pbar(3,1)], '--');
% hold off


%%  Realisierung der Markovkette
%
figure(2)
load('RandStationaer');
rng(RandStationaer);
subplot(6,1,[1,2])
kmax=40;
Rand=rand(1,kmax+1);    % Gleichverteilung im Intervall [0, 1]
%   Trajektorie im Zeitintervall [0 kend1]
Z=[];
Z(1)=1;    % Anfangszustand bekannt
k1=1;
for k1=2:kmax+1
    switch Z(k1-1)
        case 1           % System ist im Zustand 1
            if Rand(k1) >= 0.7    %  mit Wahrscheinlichkeit 0.7
               Z(k1)=2;           %  erhält man Zustand 2
            else
               Z(k1)=3;
            end
        case 2           % System ist im Zustand 2
            if Rand(k1) >= 0.5
               Z(k1)=2;
            else
               Z(k1)=1;
            end
        case 3               % System ist im Zustand 3
              Z(k1)=1;
    end 
end
Time1=[0:length(Z)-1];
fill([1.5 2.5 2.5 1.5 1.5],[0.5 0.5 3.5 3.5 0.5],'y');
hold on
fill([35.5 36.5 36.5 35.5 35.5],[0.5 0.5 3.5 3.5 0.5],'y');
plot(Time1, Z, 'b');
axis([0 kmax 0.5 3.5]);
xleer;
latexylabel('$$Z(k)$$');
rotateY;
yticks([1 2 3]);
hold off
%   2. Realisierung
subplot(6,1,[3,4])
Rand=rand(1,kmax+1);    % Gleichverteilung im Intervall [0, 1]
%   Trajektorie im Zeitintervall [0 kend1]
Z=[];
Z(1)=1;    % Anfangszustand bekannt
k1=1;
for k1=2:kmax+1
    switch Z(k1-1)
        case 1           % System ist im Zustand 1
            if Rand(k1) >= 0.7    %  mit Wahrscheinlichkeit 0.7
               Z(k1)=2;           %  erhält man Zustand 2
            else
               Z(k1)=3;
            end
        case 2           % System ist im Zustand 2
            if Rand(k1) >= 0.5
               Z(k1)=2;
            else
               Z(k1)=1;
            end
        case 3               % System ist im Zustand 3
              Z(k1)=1;
    end 
end
Time1=[0:length(Z)-1];
fill([1.5 2.5 2.5 1.5 1.5],[0.5 0.5 3.5 3.5 0.5],'y');
hold on
fill([35.5 36.5 36.5 35.5 35.5],[0.5 0.5 3.5 3.5 0.5],'y');
plot(Time1, Z, 'b');
axis([0 kmax 0.5 3.5]);
xleer;
latexylabel('$$Z(k)$$');
rotateY;
yticks([1 2 3]);
hold off
%   3. Realisierung
subplot(6,1,[5,6])
Rand=rand(1,kmax+1);    % Gleichverteilung im Intervall [0, 1]
%   Trajektorie im Zeitintervall [0 kend1]
Z=[];
Z(1)=1;    % Anfangszustand bekannt
k1=1;
for k1=2:kmax+1
    switch Z(k1-1)
        case 1           % System ist im Zustand 1
            if Rand(k1) >= 0.7    %  mit Wahrscheinlichkeit 0.7
               Z(k1)=2;           %  erhält man Zustand 2
            else
               Z(k1)=3;
            end
        case 2           % System ist im Zustand 2
            if Rand(k1) >= 0.5
               Z(k1)=2;
            else
               Z(k1)=1;
            end
        case 3               % System ist im Zustand 3
              Z(k1)=1;
    end 
end
Time1=[0:length(Z)-1];
fill([1.5 2.5 2.5 1.5 1.5],[0.5 0.5 3.5 3.5 0.5],'y');
hold on
fill([35.5 36.5 36.5 35.5 35.5],[0.5 0.5 3.5 3.5 0.5],'y');
plot(Time1, Z, 'b');
axis([0 kmax 0.5 3.5]);
%latexxlabel('$$k$$');
xticks([0 10 20 30 40]);
latexylabel('$$Z(k)$$');
rotateY;
yticks([1 2 3]);
hold off


%%
epsfigc('Stationaer');
