%              MARKOV.M
%
%  Verhalten eines Markovprozesses
%  
%
% J. Lunze
%
% 19.5.2006
% für 2. Auflage: 21.6.2012
% für 3. Auflage: 26.1.2017
% für 4. Auflage: 12.1.2025
%
%
%   Berechnung des Verhaltens eines kontinuierlichen Markovprozesses
%   und Vergleich mit einer eingebetteten Markovkette
%   und einer zeitdiskreten Darstellung des kont. Markovprozesses
%
%  Q = Systemmatrix des Markovprozesses
%  G = Matrix der Zustandsübergangswahrscheinlichkeiten der Markovkette
%
%
clear
close all

%%  
lambda1=1;
lambda2=2;
lambda3=0.2;
g21=0.95;
g31=1-g21;
Q=[-lambda1        0       0
   g21*lambda1 -lambda2  lambda3
   g31*lambda1  lambda2 -lambda3];
Markovprozess=ss(Q, zeros(3,1), eye(3,3), zeros(3,1));
p0=[1; 0; 0];                   % Anfangswahrscheinlichkeitsverteilung
Tmax=4;
T1=[0:0.01:Tmax];
[P1, T1]=initial(Markovprozess, p0, T1);
figure(1)
%
subplot(4,1,[1,2])
plot(T1, P1,'b');
hold on
%latexxlabel('$$t$$ in s');
xleer;
latexylabel('$$p_i$$');
rotateY;
yticks([0 0.5 1]);
axis([0 Tmax 0 1]);
latextext(0.5, 0.85, '$$p_1$$');
latextext(1.3, 0.85, '$$p_3$$');
latextext(3.5, 0.3, '$$p_2$$');
%   Berechnung der Aufenthaltswahrscheinlichkeiten zur Zeit Tquer
Tquer=0.5;
PTquer=expm(Q*Tquer)*p0;
plot(0.5, PTquer, 'bo');

%  Berechnung der stationären Verteilung
[V, D]=eig(Q);
k=0;
abbruch=-4;
while abbruch<0
    k=k+1;
    abbruch=D(k,k);
end
v=V(1:3,k);
pquer=v/(v(1)+v(2)+v(3))
%
%  Veränderte Parameter
lambda1=1;
lambda2=0.2;
lambda3=2;
Q=[-lambda1        0       0
   g21*lambda1 -lambda2  lambda3
   g31*lambda1  lambda2 -lambda3];
Markovprozess=ss(Q, zeros(3,1), eye(3,3), zeros(3,1));
p0=[1; 0; 0];                   % Anfangswahrscheinlichkeitsverteilung
T=[0:0.01:Tmax];
[P, T]=initial(Markovprozess, p0, T);
subplot(4,1,[3,4])
plot(T, P,'b');
latexxlabel('$$t$$ in s');
xticks([0 1 2 3 4]);
latexylabel('$$p_i$$');
rotateY;
yticks([0 0.5 1]);
axis([0 Tmax 0 1]);
latextext(0.5, 0.85, '$$p_1$$');
latextext(0.9, 0.85, '$$p_2$$');
latextext(3.5, 0.3, '$$p_3$$');

%%  Berechnung der stationären Verteilung
[V, D]=eig(Q);
k=0;
abbruch=-4;
while abbruch<0
    k=k+1;
    abbruch=D(k,k);
end
v=V(1:3,k);
pquer=v/(v(1)+v(2)+v(3))


%%  Vergleich mit der eingebetteten Markovkette
%
figure(2)
subplot(2,1,1)
G=[0   0 0
   g21 0 1
   g31 1 0];
%                           zeitdiskrete Markovkette
eingebMarkovkette=ss(G, zeros(3,1), eye(3,3), zeros(3,1), 1);
p0=[1; 0; 0];                   % Anfangswahrscheinlichkeitsverteilung
kmax=6;
P=dinitial(G, zeros(3,1), eye(3,3), zeros(3,1), p0, kmax);
stairs([0:kmax-1], P,'b', 'LineWidth',1.5);
latexxlabel('$$k$$');
xticks([0 1 2 3 4 5 6]);
latexylabel('$$p_i$$');
rotateY;
yticks([0 0.5 1]);
latextitle('Verhalten der eingebetteten Markovkette');


%%  Zeitdiskrete Betrachtung des kontinuierlichen Markovprozesses
%
figure(3)
lambda1=1; %   erster Parametersatz
lambda2=2;
lambda3=0.2;
g21=0.95;
g31=1-g21;
Q=[-lambda1        0       0
   g21*lambda1 -lambda2  lambda3
   g31*lambda1  lambda2 -lambda3];
%
subplot(4,1,[1,2])
T=0.25;     %   Abtastzeit
G=expm(Q*T);
zeitdiskrMarkovkette=ss(G, zeros(3,1), eye(3,3), zeros(3,1), 1);
p0=[1; 0; 0];                   % Anfagswahrscheinlichkeitsverteilung
kmax=ceil(Tmax/T);
P=dinitial(G, zeros(3,1), eye(3,3), zeros(3,1), p0, kmax);
plot([0:kmax-1]*T, P,'bo');
hold on
%  Zum Vergleich Eintragen der genauen Werte
plot(T1, P1, 'b--');
latextext(0.5, 0.8, '$$p_1$$');
latextext(1.3, 0.8, '$$p_3$$');
latextext(3.5, 0.3, '$$p_2$$');
latexylabel('$$p_i$$');
rotateY;
yticks([0 0.5 1]);
xleer;
axis([0 Tmax -0.1 1.1]);
latextitle('Zeitdiskrete Approximation der kontinuierlichen Markovkette');
hold off
%
%  zweite Darstellung
subplot(4,1,[3,4])
stairs([0:kmax-1]*T, P,'b', 'LineWidth',1.5);
hold on
%  Zum Vergleich Eintragen der genauen Werte
plot(T1, P1, 'b--');
latextext(0.6, 0.9, '$$p_1$$');
latextext(1.3, 0.8, '$$p_3$$');
latextext(3.5, 0.3, '$$p_2$$');
latexylabel('$$p_i$$');
rotateY;
yticks([0 0.5 1]);
latexxlabel('$$t$$ in s');
xticks([0 1 2 3 4]);
axis([0 Tmax -0.1 1.1]);
hold off


%%
epsfigc('Markov');
