%      ConvergenceDemo.M
%   
%  Figures to demonstrate the convergence of random sequences
%  
% J. Lunze
% 
% 29.11.2017
% Version v. 22.11.2018
% for 2nd edition: 28.3.2021
%
echo off
clear
close all

%%   Sequence not converging everywhere
%
figure(1)
subplot(4,1,[1 2])
kend=16;
K=[0:kend-1];
omega=1.3;
X1=[];
for k1=K
    X1(k1+1)=(-1/omega)^(k1-1);
end
plot(K, X1,'b--')
hold on
plot(K, X1, 'bo');
xleer;
latexylabel('$$X(k,1.3)\;$$');
rotateY;
axis([0 kend-1 -1.35 1.35]);
latextitle('$${\tt ConvergenceDemo:}$$ Two random sequences');
yticks([-1 0 1]);
%
subplot(4,1,[3 4])
omega=1;
X2=[];
for k1=K
    X2(k1+1)=(-1/omega)^(k1-1);
end
plot(K, X2,'b--')
hold on
plot(K, X2, 'bo');
latexxlabel('$$k$$');
latexylabel('$$X(k,1)\;$$');
rotateY;
axis([0 kend-1 -1.35 1.35]);
yticks([-1 0 1]);
xticks([0 5 10 15]);
hold off;

%%  Convergence with probability
%
%   X(k) = 1 with probability 1/k
%        = 0 otherwise
figure(2)
subplot(4,1,[1 2])
kend=90;
%RandConvergence=rng;                         % save s to restore the generator setting
%save('RandConvergence', 'RandConvergence');
load('RandConvergence');       % Seed (initialisation) of the number generator
rng(RandConvergence);
P=rand(1,kend);
X=[];               %  stochastic sequence
for k1=1:kend
    if  P(k1) <= 1/k1
        X(k1)=1;
    else
        X(k1)=0;
    end
end
plot([1:kend], X,'bo');
hold on
for k1=1:kend
   plot([k1 k1],[0 0.96*X(k1)],'b'); 
end
latexxlabel('$$k$$');
latexylabel('$$X\;$$');
rotateY;
latextitle('$${\tt ConvergenceDemo:}$$ Random sequence that converges in probability');
axis([-1 kend -0.3 1.3]);
yticks([0 1]);
xticks([1 20 40 60 80]);
hold off


%%  Weak law of large numbers
%   relative frequency of getting a "6" in dice throwing
figure(3)
N1=201;
P=[];
S=[];     % overall success: sum X(k)
R=[];    % Difference S(k) - pk
%RandConvergence2=rng;                         % save s to restore the generator setting
%save('RandConvergence2', 'RandConvergence2');
load('RandConvergence2');       % Seed (initialisation) of the number generator
rng(RandConvergence2);
p=1/6;
P=rand(1,N1);
S=[];
Sigma=[];
k2=0;
if P(1)<1/6
    S(1)=1;
    k2=k2+1;
else
    S(1)=0;
end
for k1=2:length(P)
    if P(k1) < 1/6
       S(k1)=S(k1-1)+1;
       k2=k2+1;
    else
        S(k1)=S(k1-1);
    end
end
for k1=1:length(P)
    Sigma(k1)=sqrt((k1+1)*p*(1-p));
end
Time=[0:length(P)-1];
Exp=6*p*Time;        % 6*p =1
PolygonX=[Time, flip(Time), Time(1,1)];
PolygonY=[Exp+2*6*Sigma, flip(Exp-2*6*Sigma), Exp(1,1)+2*6*Sigma(1,1)];
subplot(4,1,[1 2])
fill(PolygonX, PolygonY, [0.9 0.9 0.9]);
hold on
plot(Time, 6*S, 'b', [1, length(P)], [0 6*p*length(P)],'k:');
latexylabel('$$S\;$$');
rotateY;
latexxlabel('$$k$$');
%xleer;
xticks([0 50 100 150 200]);
yticks([0 100 200]);
latextitle('$${\tt ConvergenceDemo:}$$ Dice throwing');
axis([-1 200 -20 270]);
hold off

%%   Estimation of the event probability
%
figure(4)
for k1=1:length(P)
    R(k1)=k1-6*S(k1);
end
for k1=1:length(P)
    H(k1)=S(k1)/k1;
end
H4=H(1,4);
H100=H(1,101);
Z24=R(4);
Z2100=R(101);
subplot(4,1,[1 2])
plot(Time, H, 'b', [1, length(P)], [1/6 1/6],'k:');
hold on;
plot(4,H4,'bo');
plot(100, H100, 'bo');
hold on
latextext(1.02*N1, 1/6, '$$\frac{1}{6}$$');
latexylabel('$$R\;$$');
rotateY;
%latexxlabel('$$k$$');
xleer;
xticksempty([0 50 100 150 200 250 300]);yticks([0 0.2 0.4 0.6]);
latextitle('$${\tt ConvergenceDemo:}$$ Dice throwing');
axis([-1 200 0 0.6]);
hold off
%
subplot(4,1,[3 4])
plot(Time, R, 'b');
hold on;
plot(4,Z24,'bo');
plot(100, Z2100, 'bo');
latexxlabel('$$k$$');
latexylabel('$$D\;$$');
rotateY;
axis([-1 200 -60 30]);
yticks([-50 -25 0 25 50]);
xticks([0 50 100 150 200 250 300]);
hold off

%%    Variance
%
figure(5)
Sigma=[];
for k1=1:length(P)
    Sigma(k1)=sqrt(5/36/k1);
end
subplot(4,1,[1 2])
PolygonX=[Time, flip(Time), Time(1,1)];
PolygonY=[1/6+2*Sigma, flip(1/6-2*Sigma), 1/6+2*Sigma(1,1)];
fill(PolygonX, PolygonY, [0.9 0.9 0.9]);
hold on
%plot([1:length(P)], 1/6+2*Sigma, 'k', [1:length(P)], 1/6-2*Sigma, 'k');
plot([0, length(P)], [1/6 1/6], 'k:');
%   Realisation from figure(3)
plot(Time, H, 'b');
latextext(1.02*N1, 1/6, '$$\frac{1}{6}$$');
latexylabel('$$R\;$$');
rotateY;
latexxlabel('$$k$$');
xticks([0 50 100 150 200 250 300]);
yticks([-0.2 0 0.2 0.4 0.6]);
latextitle('$${\tt ConvergenceDemo:}$$ Dice throwing');
axis([0 200 -0.2 0.6]);
hold off



%%  Weak law of large numbers
%   arithmetic average converges i.P.
%
figure(6)
N1=901;
P=[];
S=[];     % overall success: sum X(k)
%RandConvergence2=rng;                         % save s to restore the generator setting
%save('RandConvergence2', 'RandConvergence2');
load('RandConvergence2');       % Seed (initialisation) of the number generator
rng(RandConvergence2);
p=0.5;
P=rand(N1,1);
S=0;
for k1=1:length(P)
    if P(k1) < p
       X(k1)=1;
    else
        X(k1)=-1;
    end
    S=S+X(k1);
    Xbar(k1)=S/k1;
end
Time=[0:length(P)-1];
subplot(3,1,1)
plot(Time, Xbar, 'b');
hold on
latexylabel('$$\overline{X}\;\;\;$$');
rotateY;
%latexxlabel('$$k$$');
%xleer;
xticks([0 50 100 150 200]);
yticks([-1 0 1]);
latextitle('$${\tt ConvergenceDemo:}$$ Convergence i.P.');
axis([-1 100 -1 1]);
hold off
%
subplot(3,1,2)
plot(Time, Xbar, 'b');
hold on
latexylabel('$$\overline{X}\;\;$$');
rotateY;
%latexxlabel('$$k$$');
%xleer;
xticks([300 350 400]);
yticks([-0.1 0 0.1]);
axis([300 400 -0.1 0.1]);
hold off
%
subplot(3,1,3)
plot(Time, Xbar, 'b');
hold on
latexylabel('$$\overline{X}\;$$');
rotateY;
latexxlabel('$$k$$');
%xleer;
xticks([800 850 900]);
yticks([-0.01 0 0.01]);
axis([800 900 -0.01 0.01]);
hold off



%%  Stochastic stability of a re-set system
%   = convergence w.p.1
%
%   X(k+1) = A x(k) if omega=0
%   X(k+1) = 1      if omega=1
%        = 0 otherwise
figure(7)
subplot(4,1,1)
kend=100;
A=0.8;
%   Stochastic experiment
ap=0.9;
P=[];
for k1=0:kend
    P(k1+1,1)=ap^k1;
end
plot(0:kend, P, 'bo');
xleer;
latexylabel('$$p_1\;$$');
rotateY;
yticks([0 0.5 1]);
axis([0 kend 0 1.1]);
latextitle('$${\tt ConvergenceDemo:}$$ Behaviour of a re-set system');
%RandConvergence=rng;             % save s to restore the generator setting
%save('RandConvergence2', 'RandConvergence2');
load('RandConvergence2');       % Seed (initialisation) of the number generator
rng(RandConvergence2);
Prob=rand(kend+1,1);
for k1=0:kend
     Omega(k1+1,1)=ceil(Prob(k1+1,1)-(1-P(k1+1,1)));
end
subplot(4,1,2)
plot(0:kend, Omega, 'bo')
latexylabel('$$\omega(k)\;$$');
rotateY;
xleer;
yticks([0 1]);
axis([0 kend -0.2 1.2]);
%
%   System behaviour
subplot(4,1,[3,4]);
X=[];               %  stochastic sequence
X(1)=1;
for k1=1:kend
    if  Omega(k1,1)==1
        X(k1+1)=1;
    else
        X(k1+1)=A*X(k1);
    end
end
plot(0:kend, X,'bo');
hold on
plot(0:kend, X, 'b--');
plot([0 kend],[0 0],'k:');
latexxlabel('$$k$$');
latexylabel('$$\quad\quad X\;\;$$');
rotateY;
axis([0 kend  -0.05 1.1]);
yticks([0 0.5 1]);
xticks([0 20 40 60 80 100]);
hold off


%%  Martingale and supermartingale
%
figure(8)
N1=31;
X=[];     % random sequence
V=[];     % supermartingale sequence
%RandConvergence3=rng;                         % save s to restore the generator setting
%save('RandConvergence3', 'RandConvergence3');
load('RandConvergence3');       % Seed (initialisation) of the number generator
rng(RandConvergence3);
%    Martingale
X=rand(N1,1)-0.5;
V(1,1)=X(1,1);
for k1=2:length(X)
    V(k1,1)=V(k1-1,1)+X(k1,1);
end
Time=[0:length(X)-1];
subplot(5,1,[1,2])
plot(Time, X, 'bo');
hold on
plot([0 N1],[0 0],'k--');
latexylabel('$$X\;$$');
rotateY;
xleer;
xticksempty([0 10 20 30]);
yticks([-1 0 1]);
latextitle('$${\tt ConvergenceDemo:}$$ Martingale');
axis([0 N1-1 -1 1]);
hold off
%
subplot(5,1,[3,4])
plot(Time, V, 'b');
hold on
plot(Time, V, 'bo');
latexylabel('$$V\;$$');
rotateY;
latexxlabel('$$k$$');
xticks([0 10 20 30]);
yticks([-2 0 2]);
axis([0 N1-1 -2 2]);
hold off



%%  Supermartingale
%
figure(9)
N1=31;
X=[];     % random sequence
V=[];     % supermartingale sequence
%RandConvergence4=rng;                         % save s to restore the generator setting
%save('RandConvergence4', 'RandConvergence4');
load('RandConvergence4');       % Seed (initialisation) of the number generator
rng(RandConvergence4);
mu=-0.3;
X=rand(N1,1)-0.5+mu;
V(1,1)=X(1,1);
for k1=2:length(X)
    V(k1,1)=V(k1-1,1)+X(k1,1);
end
Time=[0:length(X)-1];
subplot(5,1,[1,2])
plot(Time, X, 'bo');
hold on
plot([0 N1],[mu mu],'k--');
latexylabel('$$X\;$$');
rotateY;
xleer;
xticksempty([0 10 20 30]);
yticks([-1 0 1]);
latextitle('$${\tt ConvergenceDemo:}$$ Supermartingale');
axis([0 N1-1 -1 1]);
latextext(20,0.7, ['$$\mu=\,$$',num2str(mu)]);
hold off
%
subplot(5,1,[3,4])
plot(Time, V, 'b');
hold on
plot(Time, V, 'bo');
latexylabel('$$V\;$$');
rotateY;
latexxlabel('$$k$$');
xticks([0 10 20 30]);
yticks([-10 -5 0]);
axis([0 N1-1 -10 0]);
hold off





%%   Figures
epsfigc15('ConvergenceDemo');
